期望最大化（EM）是一种对数据进行分类的概率方法。如果我错了，如果它不是分类器，请纠正我。
这种 EM 技术的直观解释是什么？这里是expectation什么，是什么maximized？
问问题
44629 次
8 回答
140
注意：这个答案背后的代码可以在这里找到。
假设我们从两个不同的组（红色和蓝色）中采样了一些数据：
在这里，我们可以看到哪个数据点属于红色或蓝色组。这使得很容易找到表征每个组的参数。例如，红色组的平均值在 3 左右，蓝色组的平均值在 7 左右（如果需要，我们可以找到确切的平均值）。
一般而言，这称为最大似然估计。给定一些数据，我们计算最能解释该数据的参数（或多个参数）的值。
现在想象一下，我们看不到哪个值是从哪个组中采样的。一切对我们来说都是紫色的：
这里我们知道有两组值，但我们不知道任何特定值属于哪个组。
我们还能估计最适合这个数据的红色组和蓝色组的均值吗？
是的，我们通常可以！期望最大化为我们提供了一种方法。该算法背后的一般思想是这样的：
- 从每个参数可能的初始估计开始。
- 计算每个参数产生数据点的可能性。
- 根据参数产生的可能性，计算每个数据点的权重，指示它是更红还是更蓝。将权重与数据（期望）结合起来。
- 使用权重调整后的数据（最大化）计算更好的参数估计值。
- 重复步骤 2 到 4，直到参数估计收敛（该过程停止产生不同的估计）。
这些步骤需要进一步解释，所以我将逐步解决上述问题。
示例：估计均值和标准差
我将在此示例中使用 Python，但如果您不熟悉该语言，代码应该很容易理解。
假设我们有两个组，红色和蓝色，其值分布如上图所示。具体来说，每个组都包含一个从具有以下参数的正态分布中提取的值：
import numpy as np
from scipy import stats
np.random.seed(110) # for reproducible results
# set parameters
red_mean = 3
red_std = 0.8
blue_mean = 7
blue_std = 2
# draw 20 samples from normal distributions with red/blue parameters
red = np.random.normal(red_mean, red_std, size=20)
blue = np.random.normal(blue_mean, blue_std, size=20)
both_colours = np.sort(np.concatenate((red, blue))) # for later use...
这是这些红色和蓝色组的图像（以免您不得不向上滚动）：
当我们可以看到每个点的颜色（即它属于哪个组）时，就很容易估计每个组的均值和标准差。我们只需将红色和蓝色值传递给 NumPy 中的内置函数。例如：
>>> np.mean(red)
2.802
>>> np.std(red)
0.871
>>> np.mean(blue)
6.932
>>> np.std(blue)
2.195
但是如果我们看不到点的颜色怎么办？也就是说，每个点都被染成紫色，而不是红色或蓝色。
为了尝试恢复红色组和蓝色组的均值和标准差参数，我们可以使用期望最大化。
我们的第一步（上面的第1步）是猜测每个组的均值和标准差的参数值。我们不必聪明地猜测；我们可以选择我们喜欢的任何数字：
# estimates for the mean
red_mean_guess = 1.1
blue_mean_guess = 9
# estimates for the standard deviation
red_std_guess = 2
blue_std_guess = 1.7
这些参数估计产生如下所示的钟形曲线：
这些都是不好的估计。例如，对于合理的点组，这两种方法（垂直虚线）看起来都远离任何类型的“中间”。我们希望改进这些估计。
下一步（步骤 2）是计算每个数据点出现在当前参数猜测下的可能性：
likelihood_of_red = stats.norm(red_mean_guess, red_std_guess).pdf(both_colours)
likelihood_of_blue = stats.norm(blue_mean_guess, blue_std_guess).pdf(both_colours)
在这里，我们使用我们当前对红色和蓝色均值和标准差的猜测，简单地将每个数据点放入正态分布的概率密度函数中。例如，这告诉我们，根据我们目前的猜测，1.761 处的数据点更可能是红色（0.189）而不是蓝色（0.00003）。
对于每个数据点，我们可以将这两个似然值转换为权重（步骤 3），以便它们总和为 1，如下所示：
likelihood_total = likelihood_of_red + likelihood_of_blue
red_weight = likelihood_of_red / likelihood_total
blue_weight = likelihood_of_blue / likelihood_total
使用我们当前的估计和新计算的权重，我们现在可以计算红色和蓝色组的均值和标准差的新估计（步骤 4）。
我们使用所有数据点两次计算平均值和标准差，但权重不同：一次用于红色权重，一次用于蓝色权重。
直觉的关键是，颜色在数据点上的权重越大，数据点对该颜色参数的下一次估计的影响就越大。这具有将参数“拉”到正确方向的效果。
def estimate_mean(data, weight):
"""
For each data point, multiply the point by the probability it
was drawn from the colour's distribution (its "weight").
Divide by the total weight: essentially, we're finding where
the weight is centred among our data points.
"""
return np.sum(data * weight) / np.sum(weight)
def estimate_std(data, weight, mean):
"""
For each data point, multiply the point's squared difference
from a mean value by the probability it was drawn from
that distribution (its "weight").
Divide by the total weight: essentially, we're finding where
the weight is centred among the values for the difference of
each data point from the mean.
This is the estimate of the variance, take the positive square
root to find the standard deviation.
"""
variance = np.sum(weight * (data - mean)**2) / np.sum(weight)
return np.sqrt(variance)
# new estimates for standard deviation
blue_std_guess = estimate_std(both_colours, blue_weight, blue_mean_guess)
red_std_guess = estimate_std(both_colours, red_weight, red_mean_guess)
# new estimates for mean
red_mean_guess = estimate_mean(both_colours, red_weight)
blue_mean_guess = estimate_mean(both_colours, blue_weight)
我们对参数有新的估计。要再次改进它们，我们可以跳回第 2 步并重复该过程。我们这样做直到估计收敛，或者在执行了一些迭代之后（步骤 5）。
对于我们的数据，这个过程的前五次迭代看起来像这样（最近的迭代有更强的外观）：
我们看到均值已经在某些值上收敛，曲线的形状（由标准差控制）也变得更加稳定。
如果我们继续进行 20 次迭代，我们将得到以下结果：
EM 过程已经收敛到以下值，结果非常接近实际值（我们可以看到颜色 - 没有隐藏变量）：
| EM guess | Actual | Delta
----------+----------+--------+-------
Red mean | 2.910 | 2.802 | 0.108
Red std | 0.854 | 0.871 | -0.017
Blue mean | 6.838 | 6.932 | -0.094
Blue std | 2.227 | 2.195 | 0.032
在上面的代码中，您可能已经注意到新的标准差估计是使用先前迭代的均值估计来计算的。最终，如果我们首先计算平均值的新值并不重要，因为我们只是在某个中心点周围找到值的（加权）方差。我们仍然会看到参数的估计收敛。
于 2017-04-22T15:59:08.440 回答
37
EM 是一种在模型中的某些变量未被观察到时（即当您有潜在变量时）最大化似然函数的算法。
你可能会问，如果我们只是想最大化一个函数，为什么我们不直接使用现有的机制来最大化一个函数。好吧，如果您尝试通过获取导数并将它们设置为零来最大化它，您会发现在许多情况下，一阶条件没有解决方案。有一个先有鸡还是先有蛋的问题，要解决模型参数，您需要知道未观察到的数据的分布；但是您未观察到的数据的分布是您的模型参数的函数。
EM 试图通过迭代猜测未观察到的数据的分布来解决这个问题，然后通过最大化实际似然函数的下限来估计模型参数，并重复直到收敛：
EM算法
从猜测模型参数的值开始
E-step：对于每个具有缺失值的数据点，在给定模型参数的当前猜测和观察到的数据的情况下，使用您的模型方程求解缺失数据的分布（请注意，您正在求解每个缺失值的分布值，而不是期望值）。现在我们有了每个缺失值的分布，我们可以计算似然函数对未观察到的变量的期望值。如果我们对模型参数的猜测是正确的，那么这个预期的可能性就是我们观察到的数据的实际可能性；如果参数不正确，它将只是一个下限。
M-step：现在我们已经有了一个没有未观察到的变量的预期似然函数，像在完全观察到的情况下一样最大化该函数，以获得模型参数的新估计。
重复直到收敛。
于 2012-08-04T20:09:38.927 回答
28
这是理解期望最大化算法的简单方法：
1-阅读 Do 和 Batzoglou 的这篇EM 教程论文。
2-您的脑海中可能有问号，请查看此数学堆栈交换页面上的解释。
3-看看我用 Python 编写的这段代码，它解释了第 1 项的 EM 教程文件中的示例：
警告：代码可能混乱/次优，因为我不是 Python 开发人员。但它完成了这项工作。
import numpy as np
import math
#### E-M Coin Toss Example as given in the EM tutorial paper by Do and Batzoglou* ####
def get_mn_log_likelihood(obs,probs):
""" Return the (log)likelihood of obs, given the probs"""
# Multinomial Distribution Log PMF
# ln (pdf) = multinomial coeff * product of probabilities
# ln[f(x|n, p)] = [ln(n!) - (ln(x1!)+ln(x2!)+...+ln(xk!))] + [x1*ln(p1)+x2*ln(p2)+...+xk*ln(pk)]
multinomial_coeff_denom= 0
prod_probs = 0
for x in range(0,len(obs)): # loop through state counts in each observation
multinomial_coeff_denom = multinomial_coeff_denom + math.log(math.factorial(obs[x]))
prod_probs = prod_probs + obs[x]*math.log(probs[x])
multinomial_coeff = math.log(math.factorial(sum(obs))) - multinomial_coeff_denom
likelihood = multinomial_coeff + prod_probs
return likelihood
# 1st: Coin B, {HTTTHHTHTH}, 5H,5T
# 2nd: Coin A, {HHHHTHHHHH}, 9H,1T
# 3rd: Coin A, {HTHHHHHTHH}, 8H,2T
# 4th: Coin B, {HTHTTTHHTT}, 4H,6T
# 5th: Coin A, {THHHTHHHTH}, 7H,3T
# so, from MLE: pA(heads) = 0.80 and pB(heads)=0.45
# represent the experiments
head_counts = np.array([5,9,8,4,7])
tail_counts = 10-head_counts
experiments = zip(head_counts,tail_counts)
# initialise the pA(heads) and pB(heads)
pA_heads = np.zeros(100); pA_heads[0] = 0.60
pB_heads = np.zeros(100); pB_heads[0] = 0.50
# E-M begins!
delta = 0.001
j = 0 # iteration counter
improvement = float('inf')
while (improvement>delta):
expectation_A = np.zeros((5,2), dtype=float)
expectation_B = np.zeros((5,2), dtype=float)
for i in range(0,len(experiments)):
e = experiments[i] # i'th experiment
ll_A = get_mn_log_likelihood(e,np.array([pA_heads[j],1-pA_heads[j]])) # loglikelihood of e given coin A
ll_B = get_mn_log_likelihood(e,np.array([pB_heads[j],1-pB_heads[j]])) # loglikelihood of e given coin B
weightA = math.exp(ll_A) / ( math.exp(ll_A) + math.exp(ll_B) ) # corresponding weight of A proportional to likelihood of A
weightB = math.exp(ll_B) / ( math.exp(ll_A) + math.exp(ll_B) ) # corresponding weight of B proportional to likelihood of B
expectation_A[i] = np.dot(weightA, e)
expectation_B[i] = np.dot(weightB, e)
pA_heads[j+1] = sum(expectation_A)[0] / sum(sum(expectation_A));
pB_heads[j+1] = sum(expectation_B)[0] / sum(sum(expectation_B));
improvement = max( abs(np.array([pA_heads[j+1],pB_heads[j+1]]) - np.array([pA_heads[j],pB_heads[j]]) ))
j = j+1
于 2013-07-11T15:11:02.877 回答
18
从技术上讲，“EM”这个术语有点未指定，但我假设您指的是高斯混合建模聚类分析技术，这是一般 EM 原理的一个实例。
实际上，EM 聚类分析不是分类器。我知道有些人认为聚类是“无监督分类”，但实际上聚类分析是完全不同的。
人们对聚类分析的主要区别，以及人们对分类的最大误解是：在聚类分析中，没有“正确的解决方案”。它是一种知识发现方法，它实际上是为了发现新的东西！这使得评估非常棘手。它通常使用已知的分类作为参考进行评估，但这并不总是合适的：您拥有的分类可能会或可能不会反映数据中的内容。
让我举个例子：你有一个庞大的客户数据集，包括性别数据。将此数据集拆分为“男性”和“女性”的方法在与现有类进行比较时是最佳的。以“预测”的方式思考这很好，因为对于新用户，您现在可以预测他们的性别。在“知识发现”的思维方式中，这实际上很糟糕，因为您想在数据中发现一些新结构。然而，将数据拆分为老年人和儿童的方法在男性/女性类别中得分会尽可能低。但是，这将是一个很好的聚类结果（如果没有给出年龄）。
现在回到 EM。本质上，它假设您的数据由多个多元正态分布组成（请注意，这是一个非常强的假设，特别是当您固定集群数量时！）。然后它试图通过交替改进模型和模型的对象分配来为此找到一个局部最优模型。
为了在分类上下文中获得最佳结果，请选择大于类数的聚类数，或者甚至仅将聚类应用于单个类（以找出类中是否存在某种结构！）。
假设你想训练一个分类器来区分“汽车”、“自行车”和“卡车”。假设数据恰好包含 3 个正态分布几乎没有用处。但是，您可能会假设有不止一种类型的汽车（以及卡车和自行车）。因此，不是为这三个类别训练分类器，而是将汽车、卡车和自行车分成 10 个集群（或者可能是 10 辆汽车、3 辆卡车和 3 辆自行车，等等），然后训练一个分类器来区分这 30 个类别，然后将类结果合并回原始类。您可能还会发现有一个集群特别难以分类，例如 Trikes。它们有点像汽车，有点像自行车。或送货卡车，这更像是超大型汽车而不是卡车。
于 2012-08-04T21:45:11.330 回答
2
其他答案很好，我将尝试提供另一种观点并解决问题的直观部分。
摘录（强调我的）：
在数学中，一般来说，对偶性通常（但不总是）通过对合运算以一对一的方式将概念、定理或数学结构转换为其他概念、定理或结构：如果A 是 B，那么 B 的对偶就是 A。这样的对合有时有不动点，所以 A 的对偶就是 A 本身
通常，对象A 的对偶B以某种方式与 A 相关，以保持某种对称性或兼容性。例如 AB =常量
使用对偶性（在前面的意义上）的迭代算法的例子是：
- 最大公约数的欧几里得算法及其变体
- Gram–Schmidt 向量基算法和变体
- 算术平均 - 几何平均不等式及其变体
- 期望最大化算法及其变体（另请参阅此处了解信息几何视图）
- （..其他类似的算法..）
以类似的方式，EM 算法也可以看作是两个对偶最大化步骤：
..[EM] 被视为最大化参数和未观测变量分布的联合函数。关于参数的 M 步..
在使用对偶的迭代算法中，有一个平衡（或固定）收敛点的显式（或隐式）假设（对于 EM，这是使用 Jensen 不等式证明的）
所以这类算法的概要是：
- 类 E 步骤：找到关于给定y保持不变的最佳解决方案x 。
- M-like step (dual)：找到关于x的最佳解y（如上一步中计算的）保持不变。
- 终止/收敛步骤的标准：使用更新的x、y值重复步骤 1、2直到收敛（或达到指定的迭代次数）
请注意，当这种算法收敛到（全局）最优时，它已经找到了在两种意义上（即在x域/参数和y域/参数中）都最佳的配置。然而，该算法只能找到局部最优而不是全局最优。
我会说这是对算法大纲的直观描述
对于统计论点和应用，其他答案已经给出了很好的解释（也请查看此答案中的参考资料）
于 2015-03-03T01:22:51.580 回答
2
接受的答案引用了Chuong EM Paper，它在解释 EM 方面做得不错。还有一个youtube 视频更详细地解释了这篇论文。
回顾一下，这里是场景：
1st: {H,T,T,T,H,H,T,H,T,H} 5 Heads, 5 Tails; Did coin A or B generate me?
2nd: {H,H,H,H,T,H,H,H,H,H} 9 Heads, 1 Tails
3rd: {H,T,H,H,H,H,H,T,H,H} 8 Heads, 2 Tails
4th: {H,T,H,T,T,T,H,H,T,T} 4 Heads, 6 Tails
5th: {T,H,H,H,T,H,H,H,T,H} 7 Heads, 3 Tails
Two possible coins, A & B are used to generate these distributions.
A & B have an unknown parameter: their bias towards heads.
We don't know the biases, but we can simply start with a guess: A=60% heads, B=50% heads.
在第一次试验的问题的情况下，直觉上我们会认为 B 生成了它，因为正面的比例非常符合 B 的偏见......但这个值只是一个猜测，所以我们不能确定。
考虑到这一点，我喜欢这样考虑 EM 解决方案：
- 每次翻转试验都会“投票”选出它最喜欢的硬币
- 这是基于每枚硬币与其分布的匹配程度
- 或者，从硬币的角度来看，相对于另一枚硬币（基于对数似然） ，人们对看到这个试验的期望很高。
- 根据每个试验对每个硬币的喜欢程度，它可以更新对该硬币参数（偏差）的猜测。
- 试验越喜欢硬币，它就越能更新硬币的偏见以反映它自己！
- 本质上，硬币的偏差是通过在所有试验中组合这些加权更新来更新的，这个过程称为（最大化），它指的是在给定一组试验的情况下尝试对每个硬币的偏差进行最佳猜测。
这可能过于简单化了（甚至在某些层面上是根本错误的），但我希望这在直观层面上有所帮助！
于 2016-09-30T02:20:17.750 回答
1
EM 用于最大化具有潜在变量 Z 的模型 Q 的似然性。
这是一个迭代优化。
theta <- initial guess for hidden parameters
while not converged:
#e-step
Q(theta'|theta) = E[log L(theta|Z)]
#m-step
theta <- argmax_theta' Q(theta'|theta)
e-step：给定 Z 的当前估计，计算预期的对数似然函数
m-step：找到使这个 Q 最大化的 theta
GMM 示例：
e-step：在给定当前 gmm 参数估计的情况下估计每个数据点的标签分配
m-step：在给定新标签分配的情况下最大化新的 theta
K-means 也是一种 EM 算法，在 K-means 上有很多解释动画。
于 2012-08-04T14:56:43.243 回答
1
使用在 Zhubarb 的回答中引用的 Do 和 Batzoglou 的同一篇文章，我在Java中为该问题实现了 EM 。对他的回答的评论表明，该算法陷入了局部最优，如果参数 thetaA 和 thetaB 相同，我的实现也会发生这种情况。
下面是我的代码的标准输出，显示了参数的收敛。
thetaA = 0.71301, thetaB = 0.58134
thetaA = 0.74529, thetaB = 0.56926
thetaA = 0.76810, thetaB = 0.54954
thetaA = 0.78316, thetaB = 0.53462
thetaA = 0.79106, thetaB = 0.52628
thetaA = 0.79453, thetaB = 0.52239
thetaA = 0.79593, thetaB = 0.52073
thetaA = 0.79647, thetaB = 0.52005
thetaA = 0.79667, thetaB = 0.51977
thetaA = 0.79674, thetaB = 0.51966
thetaA = 0.79677, thetaB = 0.51961
thetaA = 0.79678, thetaB = 0.51960
thetaA = 0.79679, thetaB = 0.51959
Final result:
thetaA = 0.79678, thetaB = 0.51960
下面是我在（Do and Batzoglou，2008）中解决问题的 EM Java 实现。实现的核心部分是运行 EM 直到参数收敛的循环。
private Parameters _parameters;
public Parameters run()
{
while (true)
{
expectation();
Parameters estimatedParameters = maximization();
if (_parameters.converged(estimatedParameters)) {
break;
}
_parameters = estimatedParameters;
}
return _parameters;
}
下面是整个代码。
import java.util.*;
/*****************************************************************************
This class encapsulates the parameters of the problem. For this problem posed
in the article by (Do and Batzoglou, 2008), the parameters are thetaA and
thetaB, the probability of a coin coming up heads for the two coins A and B,
respectively.
*****************************************************************************/
class Parameters
{
double _thetaA = 0.0; // Probability of heads for coin A.
double _thetaB = 0.0; // Probability of heads for coin B.
double _delta = 0.00001;
public Parameters(double thetaA, double thetaB)
{
_thetaA = thetaA;
_thetaB = thetaB;
}
/*************************************************************************
Returns true if this parameter is close enough to another parameter
(typically the estimated parameter coming from the maximization step).
*************************************************************************/
public boolean converged(Parameters other)
{
if (Math.abs(_thetaA - other._thetaA) < _delta &&
Math.abs(_thetaB - other._thetaB) < _delta)
{
return true;
}
return false;
}
public double getThetaA()
{
return _thetaA;
}
public double getThetaB()
{
return _thetaB;
}
public String toString()
{
return String.format("thetaA = %.5f, thetaB = %.5f", _thetaA, _thetaB);
}
}
/*****************************************************************************
This class encapsulates an observation, that is the number of heads
and tails in a trial. The observation can be either (1) one of the
experimental observations, or (2) an estimated observation resulting from
the expectation step.
*****************************************************************************/
class Observation
{
double _numHeads = 0;
double _numTails = 0;
public Observation(String s)
{
for (int i = 0; i < s.length(); i++)
{
char c = s.charAt(i);
if (c == 'H')
{
_numHeads++;
}
else if (c == 'T')
{
_numTails++;
}
else
{
throw new RuntimeException("Unknown character: " + c);
}
}
}
public Observation(double numHeads, double numTails)
{
_numHeads = numHeads;
_numTails = numTails;
}
public double getNumHeads()
{
return _numHeads;
}
public double getNumTails()
{
return _numTails;
}
public String toString()
{
return String.format("heads: %.1f, tails: %.1f", _numHeads, _numTails);
}
}
/*****************************************************************************
This class runs expectation-maximization for the problem posed by the article
from (Do and Batzoglou, 2008).
*****************************************************************************/
public class EM
{
// Current estimated parameters.
private Parameters _parameters;
// Observations from the trials. These observations are set once.
private final List<Observation> _observations;
// Estimated observations per coin. These observations are the output
// of the expectation step.
private List<Observation> _expectedObservationsForCoinA;
private List<Observation> _expectedObservationsForCoinB;
private static java.io.PrintStream o = System.out;
/*************************************************************************
Principal constructor.
@param observations The observations from the trial.
@param parameters The initial guessed parameters.
*************************************************************************/
public EM(List<Observation> observations, Parameters parameters)
{
_observations = observations;
_parameters = parameters;
}
/*************************************************************************
Run EM until parameters converge.
*************************************************************************/
public Parameters run()
{
while (true)
{
expectation();
Parameters estimatedParameters = maximization();
o.printf("%s\n", estimatedParameters);
if (_parameters.converged(estimatedParameters)) {
break;
}
_parameters = estimatedParameters;
}
return _parameters;
}
/*************************************************************************
Given the observations and current estimated parameters, compute new
estimated completions (distribution over the classes) and observations.
*************************************************************************/
private void expectation()
{
_expectedObservationsForCoinA = new ArrayList<Observation>();
_expectedObservationsForCoinB = new ArrayList<Observation>();
for (Observation observation : _observations)
{
int numHeads = (int)observation.getNumHeads();
int numTails = (int)observation.getNumTails();
double probabilityOfObservationForCoinA=
binomialProbability(10, numHeads, _parameters.getThetaA());
double probabilityOfObservationForCoinB=
binomialProbability(10, numHeads, _parameters.getThetaB());
double normalizer = probabilityOfObservationForCoinA +
probabilityOfObservationForCoinB;
// Compute the completions for coin A and B (i.e. the probability
// distribution of the two classes, summed to 1.0).
double completionCoinA = probabilityOfObservationForCoinA /
normalizer;
double completionCoinB = probabilityOfObservationForCoinB /
normalizer;
// Compute new expected observations for the two coins.
Observation expectedObservationForCoinA =
new Observation(numHeads * completionCoinA,
numTails * completionCoinA);
Observation expectedObservationForCoinB =
new Observation(numHeads * completionCoinB,
numTails * completionCoinB);
_expectedObservationsForCoinA.add(expectedObservationForCoinA);
_expectedObservationsForCoinB.add(expectedObservationForCoinB);
}
}
/*************************************************************************
Given new estimated observations, compute new estimated parameters.
*************************************************************************/
private Parameters maximization()
{
double sumCoinAHeads = 0.0;
double sumCoinATails = 0.0;
double sumCoinBHeads = 0.0;
double sumCoinBTails = 0.0;
for (Observation observation : _expectedObservationsForCoinA)
{
sumCoinAHeads += observation.getNumHeads();
sumCoinATails += observation.getNumTails();
}
for (Observation observation : _expectedObservationsForCoinB)
{
sumCoinBHeads += observation.getNumHeads();
sumCoinBTails += observation.getNumTails();
}
return new Parameters(sumCoinAHeads / (sumCoinAHeads + sumCoinATails),
sumCoinBHeads / (sumCoinBHeads + sumCoinBTails));
//o.printf("parameters: %s\n", _parameters);
}
/*************************************************************************
Since the coin-toss experiment posed in this article is a Bernoulli trial,
use a binomial probability Pr(X=k; n,p) = (n choose k) * p^k * (1-p)^(n-k).
*************************************************************************/
private static double binomialProbability(int n, int k, double p)
{
double q = 1.0 - p;
return nChooseK(n, k) * Math.pow(p, k) * Math.pow(q, n-k);
}
private static long nChooseK(int n, int k)
{
long numerator = 1;
for (int i = 0; i < k; i++)
{
numerator = numerator * n;
n--;
}
long denominator = factorial(k);
return (long)(numerator / denominator);
}
private static long factorial(int n)
{
long result = 1;
for (; n >0; n--)
{
result = result * n;
}
return result;
}
/*************************************************************************
Entry point into the program.
*************************************************************************/
public static void main(String argv[])
{
// Create the observations and initial parameter guess
// from the (Do and Batzoglou, 2008) article.
List<Observation> observations = new ArrayList<Observation>();
observations.add(new Observation("HTTTHHTHTH"));
observations.add(new Observation("HHHHTHHHHH"));
observations.add(new Observation("HTHHHHHTHH"));
observations.add(new Observation("HTHTTTHHTT"));
observations.add(new Observation("THHHTHHHTH"));
Parameters initialParameters = new Parameters(0.6, 0.5);
EM em = new EM(observations, initialParameters);
Parameters finalParameters = em.run();
o.printf("Final result:\n%s\n", finalParameters);
}
}
于 2014-11-22T00:03:23.333 回答