?p=18794

为了在统计过程中发现更多有趣的结果，我们将解决极大似然估计没有简单分析表达式的情况。举例来说，如果我们混合了各种分布，

作为说明，我们可以使用样例数据

第一步是编写混合分布的对数似然函数


  
   
    
     
    
    
     
      > logL=function(theta){
     
    
   
    
     
    
    
     
      + p=theta[1]
     
    
   
    
     
    
    
     
      + m1=theta[2]
     
    
   
    
     
    
    
     
      + s1=theta[3]
     
    
   
    
     
    
    
     
      + m2=theta[4]
     
    
   
    
     
    
    
     
      + s2=theta[5]
     
    
   
    
     
    
    
     
      + logL=-sum(log(p
      *dnorm(X,m1,s1)+(1-p)*dnorm(X,m2,s2)))
     
    
   
    
     
    
    
     
      + return(logL)
     
    
   
    
     
    
    
     
      + }

极大似然性的最简单函数如下（从一组初始参数开始，只是为了获得梯度下降的起点）


  
   
    
     
    
    
     
      > 
      optim(c(.5,160,1,180,1 
      ,logL  
      >  
      theta=opt$par)
     
    
   
    
     
    
    
     
      [
      1] 
      0.5987635 
      165.2547700 
      5.9410993 
      178.4856961 
      6.3547038

因为我们可以通过使用约束优化算法来做到“更好”，例如，概率一定在0到1之间。

为了可视化估计的密度，我们使用


  
   
    
     
    
    
     
      > hist(X,col="light green probability=TRUE)
     
    
   
    
     
    
    
     
      > lines(density(X )

另一个解决方案是使用EM算法。我们将从参数的初始值开始，并比较属于每个类的机会

从属于每个类别的这些概率中，我们将估算两个正态分布的参数。使用极大似然


  
   
    
     
    
    
     
      > m1=sum(
      p*X)/sum(
      p)
     
    
   
    
     
    
    
      
     
    
   
    
     
    
    
     
      + logL=-sum(
      log(
      p*dnorm(
      X,m1,s1)+(
      1-p)*dnorm(X,m2,s2)))
     
    
   
    
     
    
    
     
      + return(logL)

这个想法实际上是有一个循环的：我们估计属于这些类的概率（考虑到正态分布的参数），一旦有了这些概率，就可以重新估计参数。然后我们再次开始


  
   
    
     
    
    
      
     
    
   
    
     
    
    
     
      > for(s in 1:100){
     
    
   
    
     
    
    
      
     
    
   
    
     
    
    
     
      + p=p1/(p1+p2)
     
    
   
    
     
    
    
      
     
    
   
    
     
    
    
     
      + s1=sqrt(sum(p
      *(X-m1)^2)/sum(p))
     
    
   
    
     
    
    
     
      + s2=sqrt(sum((1-p)*(X-m2)^2)/sum(1-p))
     
    
   
    
     
    
    
      
     
    
   
    
     
    
    
     
      + }

然后，我们恢复混合分布的“最佳”参数


  
   
    
     
    
    
     
      > hist(
      X,col=
      "light green",probability=TRUE)
     
    
   
    
     
    
    
     
      > lines(
      density(
      X))