1

Rのoptim関数を使用してモデルの3つのパラメーターを最適化しようとしていますが、「optimize」関数を使用して可能なように、値の範囲を検索する方法を理解できません。forループを使用して試してみましたが、これは私の試みの中で最も成功しましたが、何らかの理由で355の値で停止しているようです。理想的には、これよりも高い組み合わせを試してみたいと思います。これに加えて、optimを何度も呼び出す関数を記述し、ベクトル化を試み、optim内の「par」引数にリスト値を入れることを試みましたが、これらすべての試みでエラーメッセージが生成されました。

"unable to evaluate at initial parameters".

簡単に言うと、「最適化」関数のように、optim関数を使用してパラメータの値の範囲を検索する方法を知っている人はいますか？

どんな助けやポインタも大歓迎です!!!

私のコードは次のようになります。対応するスケールの3つの最尤関数と、optimの使用を3回試みます。

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))

}

どうもありがとう！！

4

1 に答える 1

11

unable to evaluate at initial parameters いくつかのポイントで関数を最適化できないため、エラーが発生します（ここに多くあります）。ご了承ください：

  • sqrtデータが正でなければならないので使用します。それ以外の場合は、負の観測を削除する必要があります。dat <- dat[dat>0]
  • ログ機能と同じ問題、dat <- dat[dat > 1]
  • log(ptp1)-log(ptp)n のベクトルを n-1 のベクトルで減算するため、ここでリサイクルが発生します。私ppt1はc(1,ppt1)
  • functionのlinear場合、指数関数に大きな r を使用するため発散します (たとえば、exp(365) を参照)。

データを簡単にプロットし、関数で何が起こるかを確認できるので、R は素晴らしいと思います。たとえば、ここでwireframeは、関数の 1 つの 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)

ここに画像の説明を入力

于 2013-03-02T23:43:17.150 に答える