原文鏈接：http://tecdat.cn/?p=18794

原文出處：拓端數(shù)據(jù)部落公眾號

為了在統(tǒng)計過程中發(fā)現(xiàn)更多有趣的結(jié)果，我們將解決極大似然估計沒有簡單分析表達(dá)式的情況。舉例來說，如果我們混合了各種分布，

作為說明，我們可以使用樣例數(shù)據(jù)

1. 


2. > X=height

第一步是編寫混合分布的對數(shù)似然函數(shù)

1. > logL=function(theta){

2. + p=theta[1]

3. + m1=theta[2]

4. + s1=theta[3]

5. + m2=theta[4]

6. + s2=theta[5]

7. + logL=-sum(log(p*dnorm(X,m1,s1)+(1-p)*dnorm(X,m2,s2)))

8. + return(logL)

9. + }

極大似然性的最簡單函數(shù)如下（從一組初始參數(shù)開始，只是為了獲得梯度下降的起點）

1. > optim(c(.5,160,1,180,1 ,logL ?> ?theta=opt$par)

2. [1] 0.5987635 165.2547700 5.9410993 178.4856961 6.3547038

因為我們可以通過使用約束優(yōu)化算法來做到“更好”，例如，概率一定在0到1之間。

為了可視化估計的密度，我們使用

1. > hist(X,col="light green probability=TRUE)

2. > lines(density(X )

另一個解決方案是使用EM算法。我們將從參數(shù)的初始值開始，并比較屬于每個類的機(jī)會

1. 


2. > p=p1/(p1+p2)

從屬于每個類別的這些概率中，我們將估算兩個正態(tài)分布的參數(shù)。使用極大似然

1. > m1=sum(p*X)/sum(p)

2. 


3. + logL=-sum(log(p*dnorm(X,m1,s1)+(1-p)*dnorm(X,m2,s2)))

4. + return(logL)

這個想法實際上是有一個循環(huán)的：我們估計屬于這些類的概率（考慮到正態(tài)分布的參數(shù)），一旦有了這些概率，就可以重新估計參數(shù)。然后我們再次開始

1. 


2. > for(s in 1:100){

3. 


4. + p=p1/(p1+p2)

5. 


6. + s1=sqrt(sum(p*(X-m1)^2)/sum(p))

7. + s2=sqrt(sum((1-p)*(X-m2)^2)/sum(1-p))

8. 


9. + }

然后，我們恢復(fù)混合分布的“最佳”參數(shù)

1. > hist(X,col="light green",probability=TRUE)

2. > lines(density(X))

這相對接近我們的估計。

最受歡迎的見解

1.R語言泊松Poisson回歸模型分析案例

2.R語言進(jìn)行數(shù)值模擬：模擬泊松回歸模型

3.r語言泊松回歸分析

4.R語言對布豐投針（蒲豐投針）實驗進(jìn)行模擬和動態(tài)可視化

5.用R語言模擬混合制排隊隨機(jī)服務(wù)排隊系統(tǒng)

6.GARCH（1,1），MA以及歷史模擬法的VaR比較

7.R語言做復(fù)雜金融產(chǎn)品的幾何布朗運動的模擬

8.R語言進(jìn)行數(shù)值模擬：模擬泊松回歸模型

9.R語言對巨災(zāi)風(fēng)險下的再保險合同定價研究案例：廣義線性模型和帕累托分布Pareto distributions