1

我正在编写一些代码来对某些参数值进行最大似然估计，并且我正在尝试创建从 optim 函数中获取的参数值的曲面图，并且需要创建一个网格来执行此操作。这是我需要创建一个让我感到困惑的网格的部分，我的 MLE 函数如下所示：

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

我已经从 optim 中创建了一些数据，为我提供了相应的参数值

我尝试从 optim 函数中获取输出并将其放入 expand.grid 函数中，例如：

gridlog<-expand.grid(logs[,"r"],logs[,"sigma"],logs[,"k"])

但这所做的只是创建一个填充了所有相同值的大型矩阵。

进入 expand.grid 函数的数据从以下位置填充：

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

}
logs<-do.call(rbind,logs)

这将创建一个 300 长的对应 sigma 的 r 和 k 的矩阵

我的数据是：

c(100, 128.675595618645, 75.436115414503, 146.398449792328, 102.419994706974, 
207.397726741841, 23.4579309898438, 42.4085746569567, 119.498216389673, 
59.7845591706614, 119.37201616882, 252.047672957539, 28.3165331949818, 
57.4918213065119, 311.615538092141, 8.53779749227741, 31.5382580618134, 
115.617013730077, 43.6907812963781, 70.9139870053552, 123.004040266686, 
132.575148404208, 114.813947981006, 115.950032495637, 120.891472762661, 
97.0207348527786, 235.618894638631, 17.0936655960759, 49.4419128844531, 
112.476950569973, 58.3241789008329, 80.0300102105128, 103.248819284132, 
99.1968765946717, 113.905769052605, 143.181386861766, 62.962989192695, 
174.054591300157, 39.9156352770331, 81.8344415290292, 176.631480374326, 
51.5564038694108, 131.542259464434, 72.5981749979889, 38.9733086158719, 
126.808054274927, 73.6960412245896, 62.5484608101147, 55.539355637003, 
137.888502803112, 106.921926717155, 140.000738390606, 162.512046122238, 
26.2949484171288, 80.4110888678422, 74.0481779531392, 33.9890286552257, 
142.477859644323, 55.1820570626643, 107.242498924143, 56.8497685792794, 
143.676120209843, 84.2334844367379, 67.0330079913484, 109.96246704725, 
157.216290273118, 59.4585552091703, 67.2986524284706, 55.2529503291083, 
38.932960005221, 62.7454169122216, 210.687014199037, 38.7348882392115, 
75.6645116341029, 115.924283193145, 117.772958122253, 45.5313134644358, 
112.306998515583, 38.7001172906923, 66.1308507048062, 122.516808638813, 
38.8283932430479, 168.014298040365, 38.0902373313928, 117.414876109978, 
168.615976661456, 66.5037228223079, 94.4482610053865, 505.254990783834, 
1.05181785078369, 1.77594058056118, 4.36034444400473, 12.1485473106491, 
82.2373017835424, 58.9775202042162, 132.907299665772, 51.2346939236555, 
123.251093218535, 143.077217943039, 96.1524852870813)

任何人都可以提供的任何帮助将不胜感激！

4

1 回答 1

3 
#find optimum:
fit<-optim(par=c(r=1,k=1,sigma=1),fn=loglike,dat=dat,scale='log',
           method='Nelder-Mead',control=list(fnscale=-1))

fit$par
        r           k       sigma 
0.3911590 254.4989317   0.5159761 

# make grid around optimum with few selected sigma values:

rs<-seq(0.01,1,length=30)
ks<-seq(230,280,length=30)
sigmas<-c(0.25,0.5159761,0.75)

# this will contains all parameter combinations 
# and the corresponding likelihood values

mlegrid<-cbind(as.matrix(expand.grid(rs,ks,sigmas)),0)  #Matrix
colnames(mlegrid)<-c('r','k','sigma','likelihood')

for(i in 1:nrow(mlegrid)){ #go through all combinations
  mlegrid[i,4]<- loglike(par=mlegrid[i,1:3],dat=dat,scale='log') 
}
mlegrid[which.max(mlegrid[,4]),] 
      r           k       sigma  likelihood 
0.3855172 257.5862069   0.5159761 -74.9940496 
# almost the same as from optim 
# (differences due to sparse grid, more dense gives more accurate results)

#for interactive plots, static versions with `persp` function
library(rgl)
persp3d(x=rs,y=ks,
        z=matrix(mlegrid[mlegrid[,3]==sigmas[1],4],nrow=length(rs)),col=2)
#with sigma from optim
persp3d(x=rs,y=ks,
        z=matrix(mlegrid[mlegrid[,3]==sigmas[2],4],nrow=length(rs)),col=2) 
persp3d(x=rs,y=ks,
        z=matrix(mlegrid[mlegrid[,3]==sigmas[3],4],nrow=length(rs)),col=2)
于 2013-03-12T10:36:59.390 回答