使用R中的optim（）实现高斯混合MLE

我试图在R中使用optim（）使用R的本地数据集（来自MASS的Geyser）来实现R中高斯混合的MLE。我的代码如下。问题是，优化工作正常，但是，返回我传递给它的原始参数，并且还说它已经收敛。如果您能指出我要脱轨的地方，我将不胜感激。我的期望是它至少会产生不同的结果如果不是完全不同的话。使用R中的optim（）实现高斯混合MLE

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

来源

2015-08-19 user2007598

为什么你没有使用你传递给MLE.x的参数？（MLE.x的第一行） –

@ antoine-sac没有先生，只是业余编码 – user2007598

另外，使用'='进行赋值被认为是不好的形式，虽然[正反两方都不是那么清楚]（http： //stackoverflow.com/questions/1741820/assignment-operators-in-r-and）。国际海事组织使用'='很好，但我只是确保它是您的选择。 –

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

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

但是，当您在评论中指出代码仍然无法正常运行。我看不到任何错误，所以我说下面的代码片段在MLE.x结束：）

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

它是什么，如果有一些NA（或NaN），我们调用浏览器（这是一个非常方便的调试工具：它停止代码并调出控制台，以便我们可以探索环境。否则我们返回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

这里的问题是：-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。

但是你实际上不应该这样做，而是想想你想做什么，因为我不太清楚你为什么要适应单变量高斯。

来源

2015-08-19 09:28:31

按照您的建议更改代码。它现在提出：'警告信息： 1：在dnorm（x.vector，mu.sigma [i，1]，mu.sigma [i，2]）：产生的NaNs $ par [1] 16146.894787 10919.923359 81.029617 54.062756 6.818465 5.615605 $值 [1] -1888.043 $计数功能梯度 $收敛 [1] 1个 $消息 NULL' – user2007598

概率超过1和MLE现在是否定的。任何线索？ – user2007598

我编辑了我的答案。 –

使用R中的optim（）实现高斯混合MLE

回答

相关问题