怎么利用r语言做em算法估计混合双参数指数分布的数值模拟

建议你先看一下这本书:

Modeling Survival Data Using Frailty Models

chap 2. Some Parametric Methods

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.2 Exponential Distribution . . . . . . . . . . . . . . . . . . . 20

2.3 Weibull Distribution . . . . . . . . . . . . . . . . . . . . . 21

2.4 Extreme Value Distributions . . . . . . . . . . . . . . . . 23

2.5 Lognormal . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.6 Gamma . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.7 Loglogistic . . . . . . . . . . . . . . . . . . . . . . . . . 29

2.8 Maximum Likelihood Estimation . . . . . . . . . . . . . 30

2.9 Parametric Regression Models

chap 6. Estimation Methods for Shared Frailty Models

6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . .105

6.2 Inference for the Shared Frailty Model . . . . . . . . . . 106

6.3 The EM Algorithm . . . . . . . . . . . . . . . . . . . . . . .108

6.4 The Gamma Frailty Model . . . . . . . . . . . . . . . . . . . 110

6.5 The Positive Stable Frailty Model . . . . . . . . . . . . . . 111

6.6 The Lognormal Frailty Model . . . . . . . . . . . . . . . . . 113

6.6.1 Application to Seizure Data . . . . . . . . . . . . . . . 113

6.7 Modified EM (MEM) Algorithm for Gamma Frailty Models 114

6.8 Application

然后用最基本的package "survival"

并参考你的模型可能用到的一些functions:

survreg(formula, data, weights, subset,na.action, dist="weibull",....)

survreg.distributions include "weibull", "exponential", "gaussian",

"logistic","lognormal" and "loglogistic"

frailty(x, distribution="gamma", ...)

distribution: either the gamma, gaussian or t distribution may be specified.

frailty.gamma(x, sparse = (nclass >5), theta, df, eps = 1e-05,

method = c("em","aic", "df", "fixed"),...)

  姓名：崔升    学号：14020120005

【嵌牛导读】：

  EM作为一种经典的处理大数据的算法，是我们在学习互联网大数据时不得不去了解的一种常用算法

【嵌牛鼻子】：经典大数据算法之EM简单介绍

【嵌牛提问】：EM是一种怎么的算法，其如何去观测其中隐变量的？

【嵌牛正文】：

1. 极大似然

极大似然（Maximum Likelihood）估计为用于已知模型的参数估计的统计学方法。比如，我们想了解抛硬币是正面（head）的概率分布θθ；那么可以通过最大似然估计方法求得。假如我们抛硬币1010次，其中88次正面、22次反面；极大似然估计参数θθ值：

θ^=argmaxθl(θ)=argmaxθθ8(1−θ)2θ^=arg⁡maxθl(θ)=arg⁡maxθθ8(1−θ)2

其中，l(θ)l(θ)为观测变量序列的似然函数（likelihood function of the observation sequence）。对l(θ)l(θ)求偏导

∂l(θ)∂θ=θ7(1−θ)(8−10θ)⇒θ^=0.8∂l(θ)∂θ=θ7(1−θ)(8−10θ)⇒θ^=0.8

因为似然函数l(θ)l(θ)不是凹函数（concave），求解极大值困难。一般地，使用与之具有相同单调性的log-likelihood，如图所示

凹函数（concave）与凸函数（convex）的定义如图所示：

从图中可以看出，凹函数“容易”求解极大值，凸函数“容易”求解极小值。

2. EM算法

EM算法（Expectation Maximization）是在含有隐变量（latent variable）的模型下计算最大似然的一种算法。所谓隐变量，是指我们没有办法观测到的变量。比如，有两枚硬币A、B，每一次随机取一枚进行抛掷，我们只能观测到硬币的正面与反面，而不能观测到每一次取的硬币是否为A；则称每一次的选择抛掷硬币为隐变量。

用Y表示观测数据，Z表示隐变量；Y和Z连在一起称为完全数据( complete-data )，观测数据Y又称为不完全数据(incomplete-data)。观测数据的似然函数：

P(Y|θ)=∑ZP(Z|θ)P(Y|Z,θ)P(Y|θ)=∑ZP(Z|θ)P(Y|Z,θ)

求模型参数的极大似然估计：

θ^=argmaxθlogP(Y|θ)θ^=arg⁡maxθlog⁡P(Y|θ)

因为含有隐变量，此问题无法求解。因此，Dempster等人提出EM算法用于迭代求解近似解。EM算法比较简单，分为两个步骤：

E步（E-step），以当前参数θ(i)θ(i)计算ZZ的期望值

Q(θ,θ(i))=EZ[logP(Y,X|θ)|Y,θ(i)]Q(θ,θ(i))=EZ[log⁡P(Y,X|θ)|Y,θ(i)]

M步（M-step），求使Q(θ,θ(i))Q(θ,θ(i))极大化的θθ，确定第i+1i+1次迭代的参数的估计值θ(i+1)θ(i+1)

θ(i+1)=argmaxθQ(θ,θ(i))θ(i+1)=arg⁡maxθQ(θ,θ(i))

如此迭代直至算法收敛。关于算法的推导及收敛性证明，可参看李航的《统计学习方法》及Andrew Ng的《CS229 Lecture notes》。 这里 有一些极大似然以及EM算法的生动例子。

3. 实例

[2]中给出极大似然与EM算法的实例。如图所示，有两枚硬币A、B，每一个实验随机取一枚抛掷10次，共5个实验，我们 可以 观测到每一次所取的硬币，估计参数A、B为正面的概率θ=(θA,θB)θ=(θA,θB)，根据极大似然估计求解

如果我们 不能 观测到每一次所取的硬币，只能用EM算法估计模型参数，算法流程如图所示：

隐变量ZZ为每次实验中选择A或B的概率，则第一个实验选择A的概率为

P(z1=A|y1,θ(0))=P(z1=A|y1,θ(0))P(z1=A|y1,θ(0))+P(z1=B|y1,θ(0))=0.65∗0.450.65∗0.45+0.510=0.45P(z1=A|y1,θ(0))=P(z1=A|y1,θ(0))P(z1=A|y1,θ(0))+P(z1=B|y1,θ(0))=0.65∗0.450.65∗0.45+0.510=0.45

按照上面的计算方法可依次求出隐变量ZZ，然后计算极大化的θ(i)θ(i)。经过10次迭代，最终收敛。

4. 参考资料

[1] 李航，《统计学习方法》.

[2] Chuong B Do &Serafim Batzoglou, What is the expectation maximization algorithm?

[3] Pieter Abbeel, Maximum Likelihood (ML), Expectation Maximization (EM) .

[4] Rudan Chen, 【机器学习算法系列之一】EM算法实例分析 .

---