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

[Shengqiao Li]

Thanks.
 log.arg worked to make sure nonnegative value. But the other side
is overflow. Beyond the upper limit, the  exp(arg) is large and bessel
overflows even if exponentially scaled one is used. This will happen at
1/500 chance. Due to I'm doing 10,000 simulations, this problem blocked
me.

Regards,

Li

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