R help - Aug 2003 - Method of L-BFGS-B of optim evaluate function outside of box constraints

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

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

On Wednesday 20 August 2003 15:44, 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?
>
No idea... I can only say that a  similar behaviour of optim
with  method="L-BFGS-B" (i.e. evaluation of the function 
outside the lower/upper limits) has occurred to me too. 
It is a rare but possible behaviour.

Last June, I have sent a full description of the problem to 
r-bugs at r-project.org

Regards,

Adelchi Azzalini


-- 
Adelchi Azzalini  <azzalini at stat.unipd.it>
Dipart.Scienze Statistiche, Universit? di Padova, Italia
http://azzalini.stat.unipd.it/
(please, no ms-word/ms-excel/alike attachments)