Using optim in R

Question

I'm trying to use the optim function in R to optimise three parameters in a model, but can't figure out how to get it to search over a range of values, as is possible using the "optimize" function. I've tried doing it using a for loop, and this was the most successful of my attempts but it seems to stop at values of 355 for some reason, Ideally I'd like to try higher combinations than this. As well as this I have tried writing functions that call optim many times, tried vectorising and tried just putting list values into the "par" arguments within optim however all these attempts produced the error message

"unable to evaluate at initial parameters".

Long short does anyone know how I can use the optim function to search over a range of values for parameters as the "optimize" function will???

Any help or pointers would be GREATLY appreciated!!!

My code looks like: It's three maximum likelihood functions for corresponding scale and then three attempts at using optim!

rm(list=ls())

load('Dat.RData')

mean(dat)
var(dat)


loglike<-function(par,dat,scale)
{ ptp<-dat[1:length(dat)-1]
  ptp1<-dat[2:length(dat)]

  r<-par['r']
  k<-par['k']
  sigma<-par['sigma']

  if(scale=='log')
  {
    return(sum(dnorm(log(ptp1)-log(ptp)*exp(r-(ptp/k)),mean=0,sd=sigma,log=T)))
  }

  if (scale=='sqrt')
  {
    return(sum(dnorm(sqrt(ptp1)-sqrt(ptp)*exp(r-(ptp/k)),mean=0,sd=sigma,log=T)))
  }

  if (scale=='linear')
  {
    return(sum(dnorm(ptp1-ptp*exp(r-(ptp/k)),mean=0,sd=sigma,log=T)))
  }
}

sqrts<-c()
for(i in 1:4000){
  sqrts[i]<-optim(par=c(r=i,k=i,sigma=i),fn=loglike,dat=dat,scale='sqrt',method='Nelder-Mead',control=list(fnscale=-1))

}

logs<-c()
for(i in 1:4000){
  logs[i]<-optim(par=c(r=i,k=i,sigma=i),fn=loglike,dat=dat,scale='log',method='Nelder-Mead',control=list(fnscale=-1))

}

lins<-c()
for(i in 1:4000){
  lins[i]<-optim(par=c(r=i,k=i,sigma=i),fn=loglike,dat=dat,scale='linear',method='Nelder-Mead',control=list(fnscale=-1))

}

Many thanks!!

post some of your data? It will help. Try posting the output of head(dput(dat)) it will help people here reconstruct a portion of your data so they can run the code more easily — Simon O'Hanlon

agstudy agstudy · Accepted Answer · 2013-03-02T23:43:17

The error unable to evaluate at initial parameters is due because you optim can't evaluate your functions at some points.(many here). Note that:

You use sqrt so your data must be positive otherwise you need to remove negative observations. dat <- dat[dat>0]
Same problem with log function ,dat <- dat[dat > 1]
Recycling happen here log(ptp1)-log(ptp) because you substract a vector of n by a vector of n-1. I would replace ppt1 by c(1,ppt1)
for linear function , it diverge since you use give a big r to exponential function ( see exp(365) for example).

I think, R is great because you can easily plot your data and see what happens with your function. For example Here I am using wireframe to plot 3-dimensional surfaces of one of your function.

dat <- seq(1,100)
ptp <- head(dat,-1)
ptp1 <- c(tail(dat,-2),1)

g <- expand.grid( k = seq(0.1,2,length.out=100),     ## k between [0.1,2]
                  sigma = seq(0.1,1,length.out=100), ## sigma [0.1,1]
                  r= c(0.1,0.5,0.8,1))               ## some r points forgrouping

z <- rep(0,nrow(g))
for(i in seq_along(z))
  z[i] <- sum(dnorm(log(ptp1)-log(ptp)*exp(g[i,'r']-(ptp/g[i,'k'])),
                 mean=0,
                 sd=g[i,'sigma'],
                 log=T))
g$z <- z
any(is.infinite(g$z))     ## you can test if you have infinite value       
FALSE

wireframe(z ~ k * sigma, data = g, groups = r,
          scales = list(arrows = FALSE),
          drape = TRUE, colorkey = TRUE)

enter image description here

Using optim in R

1 Answers