期望最大化技术的直观解释是什么?
期望最大化如果一种概率的方法来分类数据。 如果不是分类,请纠正我的错误。
什么是这种EM技术的直观解释? 这里的期待是什么,什么是最大化?
假设我们有一些来自两个不同组的红色和蓝色的数据:
在这里,我们可以看到哪个数据点属于红色或蓝色组。 这可以很容易地find表征每个组的参数。 例如,红色组的平均值在3左右,蓝色组的平均值在7左右(如果需要,我们可以find确切的方法)。
一般来说,这被称为最大似然估计 。 给定一些数据,我们计算最好解释数据的一个或多个参数的值。
现在想象一下,我们看不到从哪个组采样哪个值。 一切看起来都是紫色的:
在这里,我们知道有两组值,但不是每个值来自的组。 我们仍然可以估计出最适合我们所看到的数据的红色和蓝色组的方法吗?
是的,我们经常可以: 预期最大化让我们有办法做到这一点。 algorithm背后的一般概念是这样的:
- 从每个参数可能的初始估计开始。
- 计算每个参数产生数据点( 期望 )的可能性 。
- 根据参数产生的可能性计算每个数据点的权重。
- 将这些权重与数据组合在一起以计算对参数( 最大化 )的更好估计。
- 重复步骤2到4,直到参数估计收敛(过程停止产生不同的估计值)。
其中一些步骤需要进一步解释。 这是最好的例子和图片。
例如:估算平均值和标准偏差
在这个例子中,我将使用Python,但是如果您不熟悉这种语言,代码应该相当容易理解。
假设我们有两个组 – 红色和蓝色 – 数值分布在上面的图像中。 具体而言,每个组包含一个从正态分布中绘制的值,其中包含以下参数:
import numpy as np from scipy import stats np.random.seed(110) # for reproducible random 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
但是如果我们看不到点的颜色呢? 也就是说,不是红色或蓝色,每一点对我们来说都是紫色的。
为了尝试恢复红色和蓝色组的平均值和标准偏差参数,我们可以使用期望最大化。
我们的第一步(上面的第一步 )是猜测每组的平均值和标准偏差的参数值。 我们不需要聪明地猜测,我们可以select任何我们喜欢的数字:
# 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): return np.sum(data * weight) / np.sum(weight) def estimate_std(data, weight, mean): 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过程已经收敛到以下值,这些值非常接近实际值(在这里我们可以看到颜色 – 没有隐藏variables):
| 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
当模型中的某些variables是不可观测的(即当有潜在variables时),EM是一种用于最大化似然函数的algorithm。
你可能会问,如果我们只是试图最大化一个函数,为什么我们不使用现有的机器来最大化函数。 那么,如果你试图通过把衍生物设置为零来最大化这一点,你会发现在许多情况下,一阶条件没有解决scheme。 有一个鸡与鸡蛋的问题,为了解决你的模型参数,你需要知道你的未观测数据的分布; 但是您的未观测数据的分布是您的模型参数的函数。
EM试图通过迭代猜测未观测数据的分布来估计模型参数,然后通过最大化某个对实际似然函数下界的东西来估计模型参数,并重复直到收敛:
EMalgorithm
从猜测你的模型参数的值开始
E步骤:对于每个缺失值的数据点,使用您的模型方程来求解缺失数据的分布,给出您当前对模型参数的猜测,并给出观察到的数据(注意,您正在解决每个缺失的分布价值,而不是预期的价值)。 现在我们对每个缺失值都有一个分布,我们可以计算似然函数相对于未观测variables的期望值 。 如果我们对模型参数的猜测是正确的,那么这个预期的可能性将是我们观察到的数据的实际可能性; 如果参数不正确,它只是一个下界。
M-step:现在我们已经得到了一个没有未被观察到的variables的预期似然函数,就像在完全观察的情况下那样最大化函数,得到模型参数的一个新的估计。
重复,直到收敛。
这是一个直观的方法来理解期望最大化algorithm:
1-请阅读Do和Batzoglou的EM教程 。
2-你可能会有问题,看看这个math堆栈交换页面上的解释。
3-看看我在Python中编写的代码,该代码解释了项目1的EM教程中的示例:
警告:由于我不是Python开发人员,因此代码可能是凌乱/不理想的。 但它做的工作。
import numpy as np import math #### EM 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 # EM 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
从技术上来说,“EM”这个术语有点低估,但是我假设你参考了高斯混合build模聚类分析技术,这是一般EM原理的一个实例 。
实际上, EM聚类分析不是一个分类器 。 我知道有些人认为聚类是“无监督分类”,但实际上聚类分析是完全不同的。
关键的区别,人们对聚类分析总是有大的误解分类,即聚类分析中没有“正确的解决办法” 。 这是一种知识发现的方法,实际上是为了find新的东西! 这使评估非常棘手。 它通常使用已知的分类作为参考进行评估,但这并不总是适当的:您所拥有的分类可能会或可能不会反映数据中的内容。
让我举一个例子:你有一个庞大的客户数据集,包括性别数据。 将数据集分为“男”和“女”的方法在与现有类进行比较时是最佳的。 以“预测”的思维方式来说,这是好的,至于新用户,现在可以预测他们的性别。 在“知识发现”的思维方式中,这实际上是不好的,因为你想在数据中发现一些新的结构 。 例如,将数据分成老年人和小孩的方法会比男性和女性的得分更差 。 但是,这将是一个很好的聚类结果(如果没有给出年龄)。
现在回到EM。 从本质上讲,它假定你的数据是由多个多元正态分布组成的(注意这是一个非常强的假设,尤其是当你修正了簇的数目!)。 然后,通过交替地改进模型和模型的对象分配,试图find一个局部最优的模型 。
为了在分类上下文中获得最佳结果,请select大于类数的簇数,甚至只将聚类应用于单个类(以查明类中是否存在某种结构!)。
假设你想训练一个分类器来区分“汽车”,“自行车”和“卡车”。 假设数据完全由3个正态分布构成几乎没有用处。 但是,您可能会认为有多种types的汽车 (卡车和自行车)。 因此,不是为这三个class级分类,而是将汽车,卡车和自行车分成10个集群(或者10辆汽车,3辆卡车和3辆自行车,不pipe),然后训练一个分类器来区分这30个class级,然后将类结果合并回原始类。 你也可能会发现有一个特别难以分类的聚类,例如Trikes。 他们有些汽车,有些自行车。 还是送货卡车,比卡车更像超大的车。
其他答案是好的,我会尝试提供另一个angular度和解决问题的直观部分。
EM(Expectation-Maximization)algorithm是一类使用二元性的迭代algorithm的变体
摘录(强调我的):
在math中,一般而言,二元性通常以(但不总是)通过对合操作的方式将概念,定理或math结构转换为其他概念,定理或结构:一对一方式: A是B,那么B的对偶就是A.这样的对合有时有固定点 ,所以A的对偶就是A本身
对象 A的双 B通常以某种方式与A保持一定的对称性或兼容性 。 例如AB = const
使用对偶(在前面的意义上)的迭代algorithm的例子是:
- 最大公约数的欧几里德algorithm及其变体
- Gram-Schmidt Vector Basisalgorithm和变体
- 算术平均值 – 几何平均不等式及其变体
- 期望最大化algorithm及其变体 ( 这里也可以参见信息几何视图 )
- (..其他类似的algorithm..)
以类似的方式, EMalgorithm也可以被看作两个双重最大化步骤 :
.. [EM]被看作是最大化参数和分布在未观测variables上的联合函数.E步骤使关于未观测variables的分布最大化该函数; 关于参数的M步骤
在使用对偶性的迭代algorithm中,存在收敛点(或固定点)的显式(或隐式)假设(对于EM,这是使用Jensen不等式certificate的)
所以这些algorithm的概要是:
- 类似于E的步骤:针对给定y保持不变,find最佳解x 。
- 类似M的步骤(双重):find关于x的最佳解(如前一步所计算的)保持不变。
- 终止/收敛步骤的准则:重复步骤1,2,更新x , y的值直到收敛(或达到指定的迭代次数)
注意 ,当这样的algorithm收敛到(全局)最优时,它已经find了在两个方面都是最好的configuration(即,在x域/参数和y域/参数两者中)。 然而,algorithm只能find局部最优而不是全局最优。
我会说这是对algorithm轮廓的直观描述
对于统计参数和应用,其他答案已经给出了很好的解释(在这个答案中也检查参考)
被接受的答案引用了Chuong EM Paper ,他在解释EM方面做了不错的工作。 还有一个YouTubevideo ,更详细地解释了纸张。
回顾一下,这里是这种情况:
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解决scheme:
- 每次翻转试验都可以“投票”选出最喜欢的硬币
- 这是基于每枚硬币如何适合其分配
- 或者,从硬币的angular度来看,相对于另一枚硬币(根据对数可能性 )看到这个试验的期望值很高。
- 根据每个试验对每个硬币的喜好程度,它可以更新硬币参数(偏差)的猜测。
- 试验越喜欢硬币,越能更新硬币的偏见来反映自己的想法!
- 基本上硬币的偏见是通过在所有试验中结合这些加权更新来更新的,这个过程被称为( 最大化 ),它指的是通过一系列的试验来试图对每个硬币的偏差进行最好的猜测。
这可能是过于简单化(甚至在某些层面上是根本错误的),但我希望这在直观的层面上有所帮助!
EM用于最大化潜在variablesZ的模型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:find最大化这个Q的θ
GMM例子:
e步骤:在给定当前gmm参数估计的情况下估计每个数据点的标签分配
m-step:在给定新标签分配的情况下最大化新的theta
K-means也是一种EMalgorithm,在K-means上有很多解释animation。
使用在Zhubarb的答案中引用的Do和Batzoglou的同一篇文章,我在Java中使用 EM来解决这个问题。 对他的回答的评论显示,algorithm陷入局部最优,这在我的实现中也会出现,如果参数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
下面是我的Java的Java实现来解决(Do和Batzoglou,2008)中的问题。 实现的核心部分是运行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 (ie 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)^(nk). *************************************************************************/ 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, nk); } 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); } }