期望最大化(EM)是一种概率方法进行分类的数据。 请纠正我,如果我错了,如果它不是一个分类。
这是什么EM技术的一个直观的解释? 什么是expectation
这里的是什么maximized
?
期望最大化(EM)是一种概率方法进行分类的数据。 请纠正我,如果我错了,如果它不是一个分类。
这是什么EM技术的一个直观的解释? 什么是expectation
这里的是什么maximized
?
注意:这个答案背后的代码可以找到这里 。
假设我们有两个不同的群体,红色和蓝色采样一些数据:
在这里,我们可以看到哪些数据点属于红色或蓝色组。 这可以很容易地发现,表征每个组的参数。 例如,红组的平均大约是3,蓝组的平均大约是7(如果我们希望我们能找到确切的手段)。
这是,一般来说,被称为最大似然估计 。 鉴于一些数据,我们计算一个参数(或参数)的值最好的解释数据。
现在想象一下,我们无法看到其值从组的抽样。 一切看起来紫色给我们:
在这里,我们有知识,有值的两个群体,但我们不知道任何特定值属于哪个组。
我们还可以估算红组和蓝组最适合这个数据的手段?
是的,我们经常可以! 期望最大化为我们提供了一种方式来做到这一点。 该算法背后的很一般的思路是这样的:
这些步骤需要一些进一步的解释,所以我会通过上述问题走。
我将在这个例子中使用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
在你上面的代码可能已经注意到了,标准偏差的新的估计是使用前一次迭代的估计平均计算。 最终,如果我们的平均首次计算新的价值并不重要,因为我们刚刚发现周围的一些中心点值(加权)变化。 我们仍将看到这些参数的估计值收敛。
EM是最大化似然函数的算法,当一些模型中的变量是不可观测(即当你有潜变量)。
你可能会问公平,如果我们只是想最大限度地发挥功能,为什么我们不利用现有的机械最大化的功能。 好吧,如果你试图通过求导,并将其设置为零,以最大限度地发挥这一点,你会发现,在许多情况下,第一阶条件不具备的解决方案。 有一个在一个鸡和蛋的问题解决了,你需要知道你的不可观测数据的分布模型参数; 但你不可观测数据的分布模型参数的函数。
EM试图通过反复猜测分布的未观测到的数据,然后估算通过最大化东西是实际似然函数的下界,并重复直到收敛模型参数来解决这个问题:
EM算法
有猜测开始为你的模型参数的值
E-步:对于已经缺失值的每个数据点,用你的模型公式计算给定当前的模型参数的猜测并给予所观察到的数据(丢失数据的分布注意,你要解决一个分配每个失踪值,而不是预期的数值)。 现在,我们为每个缺失值分布,我们可以对于未观测到的变量计算似然函数的期望 。 如果我们对模型参数的猜测是正确的,这可能预计将是我们观察到的数据的实际可能性; 如果参数是不正确的,这将只是一个下限。
M-步:现在我们已经得到了预期的似然函数与它没有不可观测的变量,发挥最大的功用,你会在完全遵守的情况下,让你的模型参数的新的估计。
重复,直到收敛。
这里是一个直接的食谱,了解期望最大化算法:
1-阅读本EM教程文件由Do和Batzoglou。
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
技术上的术语“EM”是有点得以确认,但我想你指的是高斯混合模型群分析技术,这是一般原则EM的实例 。
其实,EM聚类分析是不是一个分类 。 我知道有些人认为集群是“无监督分类”,但实际上聚类分析是完全不同的东西。
关键的区别,以及很大的误区分类人们总是有聚类分析是: 在集群江苏实际,没有“正确的解决方案”。 它是一个知识发现方法,它实际上意味着找到新的东西! 这使得评估非常棘手。 它是利用已知分类为基准评估的频次,但是这并不总是恰当的:你有分类可能会或可能不会反映的是数据。
让我给你举个例子:你有一个大的数据集的客户,包括性别数据。 当您使用现有的类进行比较的是分裂这个数据集为“男性”和“女性”的方法是最佳的。 这种思想的“预测”的方式是好的,为新用户现在可以预测他们的性别。 这种思想的“知识发现”的方式实际上是不好的,因为你想发现数据中的一些新的结构 。 这将如分割数据为老人和孩子的方法,但是将比分为糟糕的,因为它可以相对于男/女班获得 。 然而,这将是一个很好的聚类结果(如果不给的年龄)。
现在回到EM。 本质上,它假定您的数据是由多个多元正态分布的(注意,这是一个非常强的假设,特别是当你修复簇的数目!)。 然后,它试图通过交替改进模型和对象分配到模型中找到这个局部最优模式。
对于在分类方面最好的结果,选择集群比班数的数量变大 ,甚至应用集群单类只(以找出是否有类中的一些结构!)。
说你要训练一个分类,以分辨“汽车”,“自行车”和“卡车”。 有一个在假设数据准确组成3正态分布的用处不大。 但是,你可能会认为有不止一种类型的车 (和卡车和自行车) 的 。 因此,而不是训练分类为这三个类,您群集轿车,卡车和摩托车为10组每个(或者10辆汽车,3辆卡车和3辆自行车,等等),然后训练分类器来分辨这30类,然后合并类结果回到原来的班级。 你也可能会发现有一个集群是特别难以分类,例如三轮车。 他们有些车,有些自行车。 或者送货车,这更像是超大的汽车比卡车。
其他的答案是好的,我会尽量提供另一种视角和解决问题的直观的一部分。
EM(期望最大化)算法是一类使用迭代算法的一种变型对偶
节选(重点煤矿):
在数学中,对偶,一般来说,通过对合操作来转换的概念,定理或数学结构到其他概念,定理或结构,在一对一的方式,通常(但不总是):如果双的A是B,则双B的是A.这样对合有时有固定点 ,从而使该双A的是A本身
一般的对象物 A的双乙以某种方式,可以保留一些对称或兼容性涉及一种。 例如AB = 常数
的迭代算法的例子,采用对偶(在前面的意义上)是:
以类似的方式, EM算法也可以被看作是两个双最大化步骤 :
.. [EM]被看作是最大化的参数的关节功能和分布在未观察到的变量。该E-步骤最大化相对于分布在未观察到的变量此功能; M步骤相对于所述参数..
在迭代算法使用对偶有一个平衡(或固定)的收敛点的显式的(或隐式的)假设(为EM这是使用Jensen不等式证明)
所以,这样的算法的轮廓:
注意 ,当这样的算法收敛到(全局)最佳,它已找到一个配置,它在两个方向最好 (即在x域/参数和y域/参数)。 然而该算法能够找到一个局部最优而不是全局最优。
我会说这是该算法的外形的直观描述
对于统计参数和应用程序,其他的答案已经给出很好的解释(也为您在这个答案参考)
接受答案引用CHUONG EM纸 ,其中做了体面的工作解释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的偏见非常好......但该值只是一个猜测,所以我们不能肯定。
考虑到这一点,我喜欢把EM原液是这样的:
这可能是一个过于简单(甚至是根本错误的一些水平),但我希望这有助于在一个直观的水平!
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步骤:给定的电流估计的Z计算预期的对数似然函数
间步骤:找到THETA最大化这个问答
GMM实施例:
E步骤:给定当前GMM-参数估计为每个数据点估计标签分配
M-步:最大限度地提高新THETA给出的新标签assigments
K-手段也是EM算法有很多解释上的K-手段动画。
使用Do和Batzoglou在同一篇文章中Zhubarb的回答引用,我实现了EM在Java这个问题。 他的回答表明,该算法被卡在局部最优,这也发生我的执行如果参数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实现EM的解决(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 (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);
}
}