[R] optim error

[Dieter Menne]

 <keunhchoi <at> gmail.com> writes:

> 
> I saw a similar question but I still don't fully understand how to implement
> optim.
> 
> > vol<-rep(c(0.03, 0.5, 2, 4, 8, 16, 32), 3)
> > time<-rep(c(2,4,8),each=7)
> > p.mated<-c(0.47, 0.48, 0.43, 0.43, 0.26, 0.23, "null", 0.68, 0.62, 0.64,
> 0.58, 0.53, 0.47,
> + 0.24, 0.8, 0.79, 0.71, 0.56, 0.74, 0.8, 0.47)
> > eury<-as.data.frame(cbind(vol, time, p.mated))
> > eury<-eury[-7,]; eury
>     vol time p.mated
> 1  0.03    2    0.47
> 2   0.5    2    0.48
> 3     2    2    0.43
> >
> > p0<- c(f=0.87, b=0.1, c=150)
> > eury.fit <- function (f, time, vol)
> + {
> +             f<-p[1]; b<-p[2]; c<-p[3]
> +             p.mated = p[1] * ( (1 -
> exp(-p[2]*time))-(p[2]/(p[2]-(p[3]/vol)))
> +                            * (exp(-p[3]/vol*time)-exp(-p[2]*time)))
> + }
> > eury.opt<- optim(p0, fn=eury.fit, NULL, method = "BFGS", hessian = TRUE)
> Error in fn(par, ...) : argument "time" is missing, with no default


Thanks for providing a full data set. When I cleaned up the >, I got a totally
difference error message, saying that p was unknown. This happens frequently
when there are eggshell of first attempts in the environment, i.e. you had some
valid p floating around. So when you are lost or before posting, restart R, and
paste your example in again. And post without the >...

Use nls, not optim for this type of fits. It is possible with optim, but much
more complicated, because you have to do something like (predicted-data)^2
manually in your target function.

Be absolutely sure that your starting values are reasonable. I always plot the
data together with the curve at the start values, nls runs into bad gradient
very easily. I tried below, but I gave up, because the function is
incomprehensible, so there might be other errors in my example.

However, the main problem with you data set is the complex function you want to
fit. I did not try to get this to work, because it is way to ambitious for 20
data points, and will certainly run into some "step size decreased beyond
recognition" problem. I suggest to use a very simple function with 2 parameter
first. Only if that works, try the more complex model. You might also look at
Gabor Grothendiek's nls2 to help you find good starting values, but don't bother
if you cannot get it to work with nls first.

Other software packages are more friendly when you try the fit, but I learned
the hard way that nls-brutality is more realistic in most cases, rudely telling
me about my overparameterized models. 

Dieter

vol<-rep(c(0.03, 0.5, 2, 4, 8, 16, 32), 3)
time<-rep(c(2,4,8),each=7)
# You used "null" with quoted, converting everything to a factor
# Use NA instead 
p.mated<-c(0.47, 0.48, 0.43, 0.43, 0.26, 0.23, NA, 0.68, 0.62, 0.64,
   0.58, 0.53, 0.47, 0.24, 0.8, 0.79, 0.71, 0.56, 0.74, 0.8, 0.47)
eury<-data.frame(vol=vol, time=time, p.mated=p.mated)
eury<-eury[-7,] 
# as a more flexible alternative, you could have used na.omit
#eury<- na.omit(eury)

# do not use c as a parameters, this is a much used function
# replaced by d
pfit <- function(f,b,d,time,vol) {
  f * ( (1 - exp(-b*time))-(b/(b-(d/vol)))
     * (exp(-d/vol*time)-exp(-d*time)))    
}
# Always plot your data with the start value! nls is very picky
# when you give bad starters
start = c( f=0.87,b=0.1,d=150)
vol1 = sort(unique(eury$vol))
time1= 4
# 
plot(vol1,pfit(start[1],start[2],start[3],time1,vol1))
# This does not work yet
#nls(p.mated~pfit(f,b,d,time,vol),data = eury,  start = start, trace=TRUE)