为了在统计过程中发现更多有趣的结果,我们将解决极大似然估计没有简单分析表达式的情况。举例来说,如果我们混合了各种分布,
作为说明,我们可以使用样例数据
> X=height
第一步是编写混合分布的对数似然函数
> 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算法。我们将从参数的初始值开始,并比较属于每个类的机会
> p=p1/(p1+p2)
从属于每个类别的这些概率中,我们将估算两个正态分布的参数。使用极大似然
> 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))
这相对接近我们的估计。
最受欢迎的见解
1.R语言泊松Poisson回归模型分析案例
2.R语言进行数值模拟:模拟泊松回归模型
3.r语言泊松回归分析
4.R语言对布丰投针(蒲丰投针)实验进行模拟和动态可视化
5.用R语言模拟混合制排队随机服务排队系统
6.GARCH(1,1),MA以及历史模拟法的VaR比较
7.R语言做复杂金融产品的几何布朗运动的模拟
8.R语言进行数值模拟:模拟泊松回归模型
9.R语言对巨灾风险下的再保险合同定价研究案例:广义线性模型和帕累托分布Pareto distributions