search for: besseli

...here: library(numDeriv)   fprime <- function(z) { alp  <- 2.0165;   rho  <- 0.868;   # simplified expressions   a      <- alp-0.5   c1     <- sqrt(pi)/(gamma(alp)*(1-rho)^alp)   c2     <- sqrt(rho)/(1-rho)   t1     <- exp(-z/(1-rho))   t2     <- (z/(2*c2))^a   bes1   <- besselI(z*c2,a)   t1bes1 <- t1*bes1   c1*t1bes1*t2 }   ## Newton iteration newton_gam <- function(z) { n   <- length(z)   r   <- runif(n)   tol <- 1E-6   cdf <- vector(length=n, mode="numeric")     for (i in 1:1000)   { # numerical intergration to find the cdf     cdf  <- int...

Hi, R guys: I'm using L-BFGS-B method of optim for minimization problem. My function called besselI function which need non-negative parameter and the besselI will overflow if the parameter is too large. So I set the constraint box which is reasonable for my problem. But the point outside the box was test, and I got error. My program and the error follows. This program depends on CircStats packag...

#bug 1: besselI() for nu<0 and expon.scaled=TRUE #tested with R-devel (2007-06-17 r41981) x <- 2.3 nu <- -0.4 print(paste(besselI(x, nu, TRUE), "=", exp(-x)*besselI(x, nu, FALSE))) #fix: #$ diff bessel_i_old.c bessel_i_new.c #57c57 #< bessel_k(x, -alpha, expo) * ((ize == 1)? 2. : 2.*exp(-x)...

...) z<-as.vector(z) x<-rvonmises(100,0,2+cos(z)) x<-as.vector(x) f<-function(x){ k1<-x[1] k2<-x[2] k12<-x[3] g<-x[4] -2*sum(log( exp(k1*cos(x))*integrate(function(y)exp((k2+k12*cos(x))*cos(y)),g,2*pi,subdivisions=1000000)$value /(2*pi*integrate(function(y)besselI(k1+k12*cos(y),0)*exp(k2*cos(y)),g,2*pi,subdivisions=1000000)$value) )) } gr<-function(u){ k1<-u[1] k2<-u[2] k12<-u[3] g<-u[4] grad1<-sum(integrate(function(y)integrate(function(a)(cos(x)-cos(a))*exp(k2*cos(y)+(k1+k12*cos(y))*cos(a)),subdivisions=10000000,-pi,pi)$va...

...z  <- pars[1]       ## Constant value        alp  <- 2.0165 ; rho <- 0.868;    c    <- sqrt(pi)/(gamma(alp)*(1-rho)^alp)       for (i in 1:n)    {  t1   <- exp(-pars[1]/(1-rho))                             t2   <- (pars[1]*(1-rho)/(2*sqrt(rho)))^(alp-0.5)         bes1 <- besselI(pars[1]*sqrt(rho)/(1-rho),alp-0.5)        bes2 <- besselI(pars[1]*sqrt(rho)/(1-rho),alp-1.5)       bes3 <- besselI(pars[1]*sqrt(rho)/(1-rho),alp+0.5)      ## Equation       f   <- c*t1*t2*bes1 - runi       ## derivative       fprime   <- c*t1*t2*( -bes1/(1-rho) + (alp-0.5)*bes1/pars[1...

Here is a copy of the last few lines in base-Ex.Rout.fail: > x <- seq(3,500);yl <- c(-.3, .2) > plot(x,x, ylim = yl, ylab="",type='n', main = "Bessel Functions Y_nu(x)") > for(nu in nus){xx <- x[x > .6*nu]; lines(xx,besselY(xx,nu=nu), col = nu+2)} > legend(300,-.08, leg=paste("nu=",nus), col = nus+2, lwd=1) > > x <-

Hi R-users,   I wrote a code to evaluate double sum as follows:   ff2 <- function(bb,eta,z,k) { r <- length(z) for (i in 1:r) { sm1 <- sum((z[i]*bb/2)*(psigamma((0:k)+eta+1,deriv=0)/(factorial(0:k)*gamma((0:k)+eta+1))))  sm2 <- sum((besselI(z[i]*bb,eta)*log(z[i]*bb/2) - sm1)/besselI(z[i]*bb,eta))  sm2 } ff2(bb,eta,z,10)     but it gave me the following message:   > source(.trPaths[5], echo=TRUE, max.deparse.length=10000) Error in source(.trPaths[5], echo = TRUE, max.deparse.length = 10000) :   C:\Documents and Settings\zakry001\Ap...

Dear R-Group: Recently, I ran into a problem. I was using a function called "I.1", which evaluates the first-order modified Bessel function of the first kind, in the package "CircStats". This function is not vectorized, since it uses a couple of "if" conditions. However, when I called this function with a vector argument, I got no error/warning messages in

...has one parameter (vector[2]). So in total I want to estimate two parameters: the AR(1) paramter (vector[1]) and the distribution parameter. My function is the negative loglikelihood derived from a mixing operator. f=function(vector) -sum(log((vector[1]*data$I)+((1-vector[1])*((exp(-2*vector[2]))*besselI(2*vector[2],data$PC))))) nlm(f,p=c(-0.2354,0.00003269)) Error in nlm(f, p = c(-0.2354, 3.269e-05)) : non-finite value supplied by 'nlm' In addition: There were 50 or more warnings (use warnings() to see the first 50) warnings() Warning messages: 1: In besselI(2 * vector[2], data$PC...

Hi, I have written out the log-likelihood function to fit some data I have (called ONES20) to the non-central chi-squared distribution. >library(stats4) >ll<-function(lambda,k){x<-ONES20; 25573*0.5*lambda-25573*log(2)-sum(-x/2)-log((x/lambda)^(0.25*k-0.5))-log(besselI(sqrt(lambda*x),0.5*k-1,expon.scaled=FALSE))} > est<-mle(minuslog=ll,start=list(lambda=0.05,k=0.006)) R accepts the function definition without a problem, but gives this error when I ask for the mle of the parameters. Error in besselI(x, nu, 1 + as.logical(expon.scaled)) :...

In order to fit a probability distribution proposed by Sichel [Journal of the Royal Statistical Society. Series A (General), Vol. 137, No. 1. (1974), pp. 25-34], I need a modified Bessel function of the 2nd kind. I notice that the base package of "R" only has modified Bessel functions of the 1st and 3rd kind. Does a modified Bessel function of the 2nd kind exist anywhere? Many

Hi, I've got an ugly but fairly simple function: mdevstdev <- function(a){ l <- dnorm(a)/(1-pnorm(a)) integrand <- function(z)(abs(z-l)*dnorm(z)) inted <- integrate(integrand, a, Inf) inted[[1]]/((1- pnorm(a))*sqrt((1 + a*l - l^2))) } I wanted to quickly produce a graph of this over the range [-3,3] so I used: plotit <-function(x=seq(-3,3,0.01),...){

Hi R-users, Does R has a topic on newton's method? Thank you for the info.

...z  <- pars[1]       ## Constant value        alp  <- 2.0165 ; rho <- 0.868;    c    <- sqrt(pi)/(gamma(alp)*(1-rho)^alp)       for (i in 1:n)    {  t1   <- exp(-pars[1]/(1-rho))                             t2   <- (pars[1]*(1-rho)/(2*sqrt(rho)))^(alp-0.5)         bes1 <- besselI(pars[1]*sqrt(rho)/(1-rho),alp-0.5)        bes2 <- besselI(pars[1]*sqrt(rho)/(1-rho),alp-1.5)       bes3 <- besselI(pars[1]*sqrt(rho)/(1-rho),alp+0.5)      ## Equation       f   <- c*t1*t2*bes1 - runi       ## derivative       fprime   <- c*t1*t2*( -bes1/(1-rho) + (alp-0.5)*bes1/pars[1...

Hi, is there support (i.e. packages or similar) in R for functional data analysis as suggested by Ramsay (1982) and Besse and Ramsay (1986)? Cheers Roland -- E-Mail : Roland.Goecke at ieee.org -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or

Hello Everybody in order to get some needed results out of my function i need to get my besselI function evaluated at some values which normally gave Inf or 0 (expon.scaled NAN) back. So I would like to increase the range in R from approxamittly 1e+320 to aabout 1e+500 or something like that. Is there any possibility or pacckage to do this easily? Thank You Sebastian Kaiser Institut for Stati...

...one parameter (vector[2]). So in total I want to estimate two parameters: the AR(1) paramter (vector[1]) and the distribution parameter. My function is the negative loglikelihood derived from a mixing operator. f=function(vector) -sum(log((vector[1]*data$I)+((1-vector[1])*((exp(-2*vector[2]))*besselI(2*vector[2],data$PC))))) nlm(f,p=c(-0.2354,0.00003269)) Error in nlm(f, p = c(-0.2354, 3.269e-05)) : non-finite value supplied by 'nlm' In addition: There were 50 or more warnings (use warnings() to see the first 50) warnings() Warning messages: 1: In besselI(2 * vector[2], d...

...y out the calculation.   integrand <- function(z) { alp  <- 2.0165   rho  <- 0.868   # simplified expressions   a      <- alp-0.5   c1     <- sqrt(pi)/(gamma(alp)*(1-rho)^alp)   c2     <- sqrt(rho)/(1-rho)   t1     <- exp(-z/(1-rho))   t2     <- (z/(2*c2))^a   bes1   <- besselI(z*c2,a)   t1bes1 <- t1*bes1    c1*t1bes1*t2 }   z1  <- 20 cdf <- integrate(integrand, lower = 0, upper = z1,abs.tol = FALSE) r   <- runif(1)   ## Newton iteration z2  <- z1 - (cdf - r)/integrand(z1)   Output   > z1  <- 20   > cdf <- integrate(integrand, lower = 0, upper...

..._K in src/nmath/bessel.h (and referred to only in bessel_k.c), equal to 705.342, which is checked if expon.scaled is FALSE. (The equivalent number for bessel_i.c is 709, defined as exparg_BESS in bessel.h.) However, besselK(x) silently returns +Inf if x>705.342. This behavior is reasonable for besselI (where the problem is an overflow), but seems wrong for besselK (since the problem here is an *underflow*). The answer may be as simple as returning 0 rather than +Inf in this case. If this is appropriate, the following patch should do the job: *** bessel_k.c Thu Oct 17 13:55:55 2002 --- bess...

...The best I have come up with so far just tests for z being longer than 1 and if so, uses sapply() recursively. This is fine, except that it doesn't preserve the dimensions correctly if z is a matrix or an array. And the function doesn't work properly if z is a scalar but A is a vector. besselI() does The Right Thing (tm), but it is internal ; what is the best way to vectorize this type of function? hypergeo <- function(A,B,C,z,tol=1e-6){ if(length(z) > 1) { return(sapply(z,hypergeo,A=A,B=B,C=C,tol=tol)) } else { term <- tmp <- 1 for(n in 1:100){...