在R中使用優化

Question

我試圖在R中使用優化函數來優化模型中的三個參數，但無法弄清楚如何讓它在一系列值上進行搜索，因爲使用「優化「功能。 我嘗試過使用for循環來做這件事，這是我嘗試中最成功的一次，但由於某種原因它似乎停止在355的值，理想情況下我想嘗試比此更高的組合。 除了這個我試過，調用的Optim多次寫入功能，試圖向量化，並試圖只是把列表的值到「相提並論」論點的Optim內但所有這些嘗試所產生的錯誤信息

"unable to evaluate at initial parameters". 

長較短時，任何人都知道我可以如何使用優化功能來搜索參數值的範圍，因爲「優化」功能將？

任何幫助或指針將非常感謝！

我的代碼看起來像這樣： 這是三個最大似然函數的對應比例尺，然後三次嘗試使用優化！

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所以您的數據必須是正數，否則您需要移除負數觀察值。因爲你用n-1的載體。減去n的矢量dat <- dat[dat>0]
• 同樣的問題具有日誌功能，dat <- dat[dat > 1]
• 回收發生在這裏log(ptp1)-log(ptp)。我會替換ppt1c(1,ppt1)
• 對於linear函數，它會因使用給出一個大的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)