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

Question

我正在尝试使用 R 的本地数据集（来自 MASS 的间歇泉）使用 optim() 在 R 中为高斯混合实现 MLE。我的代码如下。问题是 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 的末尾添加了以下 sn-p：

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...所以 dnorm 返回 NaN...

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

这就是问题所在：SD 为 -20，这不太好......

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。

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