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

[Spencer Graves]

I see two problems, which I handle differently:

	  First, if arg must be nonnegative, I program fn in terms of log.arg, 
so it can never become nonpositive.

	  Second, is besselI always positive?  If yes, then program fn to 
compute log(besselI) and use that.  Standard references such as 
Abramowitz and Stegun (1970) Handbook of Mathematical Functions (US 
Gov't printing office) give asymptotic expansions that give good answers 
for values of arguments near very large or very small.  If you test the 
values of the arguments, you can develop good numbers to use.  Then 
optim should work fine.

hope this helps.  spencer graves

Shengqiao Li wrote:
> 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))
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help