[R] Method of L-BFGS-B of optim evaluate function outside of box constraints

[Shengqiao Li]

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 package.


Anyone has any idea about this?

Thanks in advance.

Li

#################### source code ###################################

Dk2<- function(pars,theta)
{
	kappa<- pars[1]; mu<- pars[2];
 	IoK<- besselI(kappa, nu=0);
	res<- besselI(2*kappa, nu=0)/2/IoK^2 -
mean(exp(kappa*cos(theta-mu)))/IoK;
	if(is.na(res)||is.infinite(res)){
		 print(pars);
		# assign("Theta", theta, env=.GlobalEnv);
	}
    return(res);
}


mse.Dk2<- function(pars, s, n)
{
	sum.est <- SSE <- numeric(2);
	j<- 0;
    while(j<=n){
      	theta<- rvm(s, pi, k=pars[1]) - pi;
      	est<- optim(par=pars, fn=Dk2, lower=c(0.001, -pi), upper=c(10,
pi), method="L-BFGS-B", theta=theta);
      	i<- 0;
        while(est$convergence!=0 && i< 30){
	          est<- optim(par=est$par, fn=Dk2, lower=c(0.001, -pi),
upper=c(10, pi), method="L-BFGS-B", theta=theta);
	          i<- i+1;
        }
        if(est$convergence!=0) {
	        #print(j);
	        next;
	     }
        else { j<- j+1; }

        #est<- nlm(p=pars, f=Dk2, theta=theta);
        mu.hat<- est$par[2];
        while(mu.hat< -pi) mu.hat<- mu.hat + 2*pi;
        while(mu.hat > pi) mu.hat<- mu.hat  -2*pi;
        est<- c(est$par[1], mu.hat);
        sum.est <- sum.est +  est;
	  	SSE <- SSE + (est - pars)^2;


	}
	Est <-  sum.est/n;
	Bias<- Est - pars;
	MSE<- SSE/n;

	res<- c(Kappa=pars[1], Kappa.hat= Est[1], Kappa.Bias=Bias[1],
Kappa.MSE=MSE[1], Mu.hat=Est[2], Mu.MSE=MSE[2])

   	return(res);
}
kappas <- c(0.01, 0.05, 0.1, 0.20, 0.5, 1, 2, 5);
N<- 10000;
for ( s in c(5, 10, 20, 30, 50)){
	cat("\nSample size = ", s);
	cat("\n=====================================\n");
	res<- NULL;
	for(i in 1:8){
		res<- rbind(res, mse.Dk2(c(kappas[i], 0), s, N));

	}
	print(round(res,4));
}

#Error message. -32.7 is far lower then the lower limit 0.001.
Sample size =  5
=====================================
[1] -32.736857  -3.141593
Error in optim(par = pars, fn = Dk2, lower = c(0.001, -pi), upper = c(10,
:
        L-BFGS-B needs finite values of fn
In addition: Warning messages:
1: NaNs produced in: besselI(x, nu, 1 + as.logical(expon.scaled))
2: NaNs produced in: besselI(x, nu, 1 + as.logical(expon.scaled))