在R中使用optim
在R中使用optim
提问于 2013-03-03 03:45:06
我正在尝试使用R中的optim函数来优化模型中的三个参数，但是我不知道如何让它在一个值范围内进行搜索，因为使用"optimize“函数是可能的。我尝试过使用for循环，这是我最成功的一次尝试，但由于某些原因，它似乎在值355处停止了，理想情况下，我想尝试比这个更高的组合。除此之外，我还多次尝试编写调用optim的函数，尝试矢量化，并尝试将列表值放入optim内的"par“参数中，但是所有这些尝试都会产生错误消息
"unable to evaluate at initial parameters".长短有人知道如何使用optim函数来搜索参数的值范围，就像"optimize“函数所做的那样？
任何帮助或指示都将不胜感激！
我的代码看起来是这样的:它是对应尺度的三个最大似然函数，然后使用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))
}非常感谢！
回答 1
Stack Overflow用户
回答已采纳
发布于 2013-03-03 07:43:17
出现错误unable to evaluate at initial parameters是因为optim在某些点上无法计算函数(这里有很多点)。请注意：
，
- ，你使用
sqrt，所以你的数据必须是正数，否则你需要删除负数。对数函数的问题，在这里发生，因为你用n的向量减去n-1的向量。，dat <- dat[dat > 1] - Recycling，
dat <- dat[dat>0] - Same，
dat <- dat[dat > 1] - Recycling，
log(ptp1)-log(ptp)。我会用c(1,ppt1) - for
linear函数来代替ppt1，因为你使用了一个大的r来表示指数函数(例如参见exp(365) )，所以它是不同的。
我认为，R很棒，因为您可以轻松地绘制数据，并查看函数发生了什么。例如，我在这里使用wireframe绘制您的一个函数的三维曲面。
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)
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/15178647
复制相关文章
相似问题
腾讯云开发者
