使用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

回答

相關問題