[R] Modified Bessel Function - 2nd kind

[Douglas Bates]

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)