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
Modified Bessel Function - 2nd kind
4 messages · Dr Andrew Wilson, David Pearson, Douglas Bates
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
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
Dr Andrew Wilson <eia018 at comp.lancs.ac.uk> writes:
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:
It appears that the value of q is being driven toward negative values. The numerical derivative routine is trying to evaluate the predictions at a negative value of a and q. This suggests that your model may not fit the data or that you need better starting values for the parameters.
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
One way of enforcing constraints on the parameters in a nonlinear
regression is to use transformed parameters. In this case you can use
a logistic transformation for q and the logarithm of a. I suggest
that you write a function for the model rather than an expression
modl <- function(x, loga, logisq) {
xm5 <- x - 0.5
a <- exp(loga)
q <- 1/(1 + exp(-logisq))
sqrt((2*a)/pi)*exp(a*sqrt(1-q))*((((a*q)/2)^x)/factorial(x))* ((pi/2) *
(besselI(a,-xm5) - besselI(a,xm5))/sin(xm5 * pi))
}
and try to fit
y2 <- nls(y ~ modl(x, loga, logisq), start=list(loga=log(0.1),
logisq=log(0.1/0.9), trace = TRUE)