[R] question on "optim"

[Hey Sky]

thanks. Ben

after read your email, I realized the initial value of w[5]=0 is a stupid 
mistake. and I have changed it. but I am sorry I cannot reproduce the result, 
convergence, as you get. the error message is <non-finite finite difference 
value [12]>. any suggestion about it?

and could you plz recommend some R books on optimization, such as tips for setup 
gradient and others, or common mistakes?  thanks

Nan







----- Original Message ----
From: Ben Bolker <bbolker at gmail.com>
To: r-help at stat.math.ethz.ch
Sent: Tue, September 7, 2010 11:15:43 AM
Subject: Re: [R] question on "optim"

Hey Sky <heyskywalker <at> yahoo.com> writes:

> I do not know how to describe my question. I am a new user for R and
> write the 
> following code for a dynamic labor economics model and use OPTIM to get 
> optimizations and parameter values. the following code does not work due to 
> the equation:
> 
>    wden[,i]<-dnorm((1-regw[,i])/w[5])/w[5]
> 
> where w[5] is one of the parameters (together with vector a, b and other 
> elements in vector w) need to be estimated. if I
>  delete the w[5] from the upper 
> equation. that is:
> 
>  wden[,i]<-dnorm(1-regw[,i])
> 
> optim will give me the estimated parameters.

  Thank you for the reproducible example!

  The problem is that you are setting the initial value of w[5]
to zero, and then trying to divide by it ...

  I find that


guess<-rep(0,times=npar)
guess[16] <- 1

system.time(r1<-optim(guess,myfunc1,data=mydata, method="BFGS",hessian=TRUE,
                      control=list(trace=TRUE)))

seems to work OK (I have no idea if the answers are sensible are not ...)

If you're going to be doing a lot of this it might be wise to see
if you can specify the gradient of your objective function for R --
it will speed up and stabilize the fitting considerably.

  By the way, you should be careful with this function: if we try
this with Nelder-Mead instead, it appears to head to a set of
parameters that lead to some sort of singularity in the objective
function:

system.time(r2<-optim(guess,myfunc1,data=mydata, 
   method="Nelder-Mead",hessian=FALSE,
                      control=list(trace=TRUE,maxit=5000)))

## still thinks it hasn't converged, but objective function is
##   much smaller

## plot 'slice' through objective space where 0 corresponds to
##  fit-1 parameters and 1 corresponds to fit-2 parameters;
## adapted from emdbook::calcslice
range <- seq(-0.1,1.1,length=400)
slicep <- seq(range[1], range[2], length = 400)
slicepars <- t(sapply(slicep, function(x) (1 - x) * r1$par +  x * r2$par))
v <- apply(slicepars, 1, myfunc1)
plot(range,v,type="l")


  Ideally, you should be able to look at the parameters of fit #2
and figure out (a) what the result means in terms of labor economics
and (b) how to keep the objective function from going there, or at
least identifying when it does.

  Ben Bolker

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.