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

Spencer Graves spencer.graves at pdf.com
Wed Aug 20 18:25:57 CEST 2003

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