在 R 中使用优化

Question

我正在尝试使用 R 中的 optim 函数来优化模型中的三个参数，但无法弄清楚如何让它搜索一系列值，就像使用“优化”函数一样。 我尝试过使用 for 循环，这是我尝试中最成功的一次，但由于某种原因它似乎停止在 355 的值，理想情况下我想尝试比这更高的组合。 除此之外，我还尝试编写多次调用 optim 的函数，尝试矢量化并尝试将列表值放入 optim 中的“par”参数中，但是所有这些尝试都产生了错误消息

"unable to evaluate at initial parameters".

有谁知道我如何使用 optim 函数来搜索参数值的范围，因为“优化”函数会？？？

任何帮助或指点将不胜感激！！！

我的代码如下所示： 它是对应尺度的三个最大似然函数，然后是使用 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))

}

非常感谢！！

Answer 1

错误unable to evaluate at initial parameters 是由于您优化无法在某些时候评估您的函数。（这里有很多）。请注意：

• 您使用sqrt，因此您的数据必须是正面的，否则您需要删除负面观察。 dat <- dat[dat>0]
• 日志功能同样的问题，dat <- dat[dat > 1]
• 这里发生了回收log(ptp1)-log(ptp)，因为您用 n-1 的向量减去了 n 的向量。我会将ppt1 替换为c(1,ppt1)
• 对于 linear 函数，它会出现分歧，因为您使用给指数函数一个很大的 r（例如，参见 exp(365)）。

我认为，R 很棒，因为您可以轻松地绘制数据并查看函数会发生什么。例如，在这里我使用wireframe 来绘制您的一个函数的 3 维表面。

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)