用optim()在R中实现高斯混合MLE

Question

我正在尝试用optim()实现混合高斯混合的MLE，使用R的本地数据集(来自质量的Geyser)。我的密码在下面。问题是optim工作得很好，但是，返回我传递给它的原始参数，并且还说它已经收敛。如果你能指出我偏离轨道的地方，我将不胜感激。我的期望是，它至少会产生不同的结果，即使不是大不相同。

library(ggplot2)
library(MASS)
data("geyser")
externaldata=geyser$waiting
x.vector=externaldata


MLE.x= function(combined.vector)
{ combined.vector=bigvec
  x.vector = externaldata
  K = k #capital K inside this MLE function, small K defined in the global environment
  prob.vector = combined.vector[1:K] 
  mu.vector =combined.vector[(K+1):((2*K))]
  sd.vector=combined.vector[(2*K+1):(3*K)]
  prob=matrix(rep(prob.vector,length(x.vector)),byrow=TRUE,nrow = length(x.vector))
  mu.sigma=cbind(mu.vector,sd.vector)
  x.by.K=matrix(nrow = length(x.vector), ncol = k)
  for (i in 1:K){
    x.by.K[,i]=dnorm(x.vector,mu.sigma[i,1],mu.sigma[i,2])
  }
  prob.mat=x.by.K*prob
  density=apply(prob.mat,1,sum)
  log.density=sum(-log(density))
  return(log.density)
}



## k=2 set ##
meanvec=c(50,80)
sigmavec=c(5,5)
k=2
probvec=c(1/3,2/3)
bigvec=c(probvec,meanvec,sigmavec)
est.k2.MLE=MLE.x(bigvec)
z=optim(bigvec,
        fn=MLE.x,
        method = "L-BFGS-B")
z


#### k=3 set #####
meanvec=c(50,70,80)
sigmavec=c(5,5,5)
k=3
probvec=rep(1/3,3)
bigvec=c(probvec,meanvec,sigmavec)
est.k3.MLE=MLE.x(bigvec)
z=optim(bigvec,
        fn=MLE.x,
        method = "BFGS")
z

Answer 1

删除MLE.x函数的第一行。

它总是返回相同的东西，因为它的参数被全局变量"bigvec“所取代。所以MLE不能收敛，我认为你应该达到最大迭代。您可以通过访问z$convergence来检查这一点，其中z是optim返回的值。这将是一个整数代码。0表示一切正常，1表示已达到最大迭代次数。其他值是不同的错误代码。

但是，正如您在注释中指出的那样，代码仍然不能正常运行。我看不到任何错误，所以在MLE.x末尾添加了以下代码片段：

if(any(is.na(density))) {
    browser()
  } else {
    log.density
  }

它所做的是，如果有NA (或NaN)，我们调用browser()，这是一个非常方便的调试工具:它停止代码并打开控制台，以便我们可以探索环境。否则我们会返回log.density。

然后我运行了代码，当密度中存在NA时，它现在打开控制台，而不是失败：

你可以看到：

Browse[1]> head(x.by.K)
     [,1]       [,2]
[1,]  NaN 0.01032407
[2,]  NaN 0.01152576
[3,]  NaN 0.01183521
[4,]  NaN 0.01032407
[5,]  NaN 0.01107446
[6,]  NaN 0.01079706

第一列x.by.K是NaN..。所以德诺姆返回NaN..。

Browse[1]> mu.sigma
     mu.vector sd.vector
[1,]  64.70180 -20.13726
[2,]  61.89559  33.34679

问题是:20岁的SD，这不太好.

Browse[1]> combined.vector
[1] 1267.90677 1663.42604   64.70180   61.89559  -20.13726   33.34679

但这是MLE.x的输入。

在这里，我刚刚向您展示了我如何调试代码:)

所以，在优化过程中，参数5和6取负值，这会导致dnorm失效。为什么他们不会是阴性的？Optim不知道这些应该是积极的！

因此，您必须找到一种使用SD>0约束进行约束优化的方法。

但实际上你不应该这么做，而是想想你想做什么，因为我不太确定为什么你想要拟合单变量高斯。