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 thanks, Andrew Wilson
Dr Andrew Wilson wrote:> > 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?G'day, According to Mathworld, you can make it from modified Bessel functions of the 1st kind: http://mathworld.wolfram.com/ModifiedBesselFunctionoftheSecondKind.html Regards, David. -- David Pearson david.pearson at mail.nerc-essc.ac.uk Environmental Systems Science Centre Tel: (0118) 9318741 Harry Pitt Building Fax: (0118) 9316413 University of Reading WWW: http://www.nerc-essc.ac.uk/~dwcp/Home.html RG6 6AL, United Kingdom
Maybe this hint is of help to you: rather to enter the parameters (a,q) directly into your equation, you could use some transformations to circumvent the inequalities, such as: 0 < q < 1 q = exp(x)/(1+exp(x)) and insert the rhs for q into your argument. Likewise instead, of a > 0 use: a = y^2 HTH, Bernhard -----Original Message----- From: Dr Andrew Wilson [mailto:eia018 at comp.lancs.ac.uk] Sent: 11 December 2002 13:58 To: r-help at stat.math.ethz.ch Subject: Re: [R] Modified Bessel Function - 2nd kind Many thanks for this pointer. Using the formula from the page you referenced, I now have the formula with the modified Bessel function of the second kind:> x <- c(1,2,3,4,5,6,7,8) > y <- c(1,4,5,7,5,4,1,1) > library(nls) > library(gregmisc) > y2 <- nls(y ~sqrt((2*a)/pi)*exp(a*sqrt(1-q))*((((a*q)/2)^x)/factorial(x))* ((pi/2) * (besselI(a,-(x-0.5)) - besselI(a,(x-0.5)))/sin((x-0.5) * pi)), start=list(a=0.1,q=0.1),trace=TRUE) However, for some reason, I get the following error message: 133.9999 : 0.100 0.001 Error in numericDeriv(form[[3]], names(ind), env) : Missing value or an Infinity produced when evaluating the model In addition: Warning messages: 1: NaNs produced in: sqrt((2 * a)/pi) 2: NaNs produced in: sqrt(1 - q) 3: NaNs produced in: besselI(x, nu, 1 + as.logical(expon.scaled)) 4: NaNs produced in: besselI(x, nu, 1 + as.logical(expon.scaled)) Could anyone tell me what I'm doing wrong (and how to fix it)? The following constraints should apply to the parameters, but I'm not aware of a "constrain" option in nls that allows me to set these minima and maxima: 0 < q < 1 a > 0 Many thanks, Andrew Wilson ______________________________________________ R-help at stat.math.ethz.ch mailing list http://www.stat.math.ethz.ch/mailman/listinfo/r-help ---------------------------------------------------------------------- If you have received this e-mail in error or wish to read our e-mail disclaimer statement and monitoring policy, please refer to http://www.drkw.com/disc/email/ or contact the sender. ----------------------------------------------------------------------
Hi Andrew: There is the "besselK()" function in the base distribution. This is exactly the Modified Bessel's function of the Second Kind. You don't need to write this in terms of I.0 and I.1. Ravi. ----- Original Message ----- From: Dr Andrew Wilson <eia018 at comp.lancs.ac.uk> Date: Wednesday, December 11, 2002 5:55 am Subject: [R] Modified Bessel Function - 2nd kind> 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 thanks, > Andrew Wilson > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > http://www.stat.math.ethz.ch/mailman/listinfo/r-help >