期望最大化（Expectation Maximization）算法簡(jiǎn)介和Python代碼實(shí)現(xiàn)（附代碼）

來(lái)源：DeepHub IMBA
本文約3400字，建議閱讀5分鐘
本文中通過(guò)幾個(gè)簡(jiǎn)單的示例解釋期望最大化算法是如何工作的。

這個(gè)算法最流行的例子（互聯(lián)網(wǎng)上討論最多的）可能來(lái)自這篇論文

（http://www.nature.com/nbt/journal/v26/n8/full/nbt1406.html）。這是一個(gè)非常簡(jiǎn)單的例子，所以我們也從這里開(kāi)始。

m?=?10?#?number?of?flips?in?experimentn?=?5?#?number?of?experiments
p_1 = 0.8p_2 = 0.3
xs = [] # number of heads in each experimentzs = [0,0,1,0,1] # which coin to flip
np.random.seed(5)for i in zs:  if i==0:      xs.append(stats.binom(n=m, p=p_1).rvs()) # flip coin 1  elif i==1:      xs.append(stats.binom(n=m, p=p_2).rvs()) # flip coin 2       xs = np.array(xs)print(xs)
exp1 = [0,1,3] # indexes of experiments with coin 1exp2 = [2,4] # indexes of experiments with coin 2
print('Analytical solutions:')print('p1: ', xs[exp1].sum() / (m * len(exp1)))print('p2: ', xs[exp2].sum() / (m * len(exp2)))

ef neg_log_likelihood(probs, m, xs, zs):  '''  compute negative log-likelihood  '''  ll = 0  for x,z in zip(xs,zs):      ll += stats.binom(p=probs[z], n=m).logpmf(x)         return -ll res = optimize.minimize(neg_log_likelihood, [0.5,0.5], bounds=[(0,1),(0,1)], args=(m,xs,zs), method='tnc')
print('Numerical solution:')print('p1: ', res.x[0])print('p2: ', res.x[1])

1. 期望步驟（E-step）：計(jì)算完整對(duì)數(shù)似然函數(shù)相對(duì)于 Z 給定數(shù)據(jù) X 的當(dāng)前條件分布和當(dāng)前參數(shù)估計(jì) theta 的條件期望；
2. 最大化步驟（M-step）：找到最大化該期望的參數(shù) theta 的值。

m = 10 # number of flips in each samplen = 5 # number of samplesxs = np.array([5,9,8,4,7])
theta = [0.6, 0.5] # initial guess p_1, p_2
for i in range(10):  p_1,p_2 = theta  T_1s = []  T_2s = []     # E-step  for x in xs:      T_1 = stats.binom(n=m,p=theta[0]).pmf(x) / (stats.binom(n=m,p=theta[0]).pmf(x) +                                                  stats.binom(n=m,p=theta[1]).pmf(x))      T_2 = stats.binom(n=m,p=theta[1]).pmf(x) / (stats.binom(n=m,p=theta[0]).pmf(x) +                                                  stats.binom(n=m,p=theta[1]).pmf(x))             T_1s.append(T_1)      T_2s.append(T_2)            # M-step  T_1s = np.array(T_1s)  T_2s = np.array(T_2s)  p_1 = np.dot(T_1s, xs) / (m * np.sum(T_1s))  p_2 = np.dot(T_2s, xs) / (m * np.sum(T_2s))  print(f'step:{i}, p1={p_1:.2f}, p2={p_2:.2f}')  theta = [p_1, p_2]

# model parametersp_1 = 0.1p_2 = 0.8tau_1 = 0.3tau_2 = 1-tau_1
m = 10 # number of flips in each samplen = 10 # number of samples
# generate datanp.random.seed(123)dists = [stats.binom(n=m, p=p_1), stats.binom(n=m, p=p_2)]xs = [dists[x].rvs() for x in np.random.choice([0,1], size=n, p=[tau_1,tau_2])]
# random initial guessnp.random.seed(123)theta = [np.random.rand() for _ in range(3)]
last_ll = 0max_iter = 100
for i in range(max_iter):  tau,p_1,p_2 = theta  T_1s = []  T_2s = []     # E-step  lls = []  for x in xs:      denom = (tau * stats.binom(n=m,p=p_1).pmf(x) + (1-tau) * stats.binom(n=m,p=p_2).pmf(x))      T_1 = tau * stats.binom(n=m,p=p_1).pmf(x) / denom      T_2 = (1-tau) * stats.binom(n=m,p=p_2).pmf(x) / denom      T_1s.append(T_1)      T_2s.append(T_2)      lls.append(T_1 * np.log(tau * stats.binom(n=m,p=p_1).pmf(x)) +                  T_2 * np.log(tau * stats.binom(n=m,p=p_2).pmf(x)))     # M-step  T_1s = np.array(T_1s)  T_2s = np.array(T_2s)  tau = np.sum(T_1s) / n  p_1 = np.dot(T_1s, xs) / (m * np.sum(T_1s))  p_2 = np.dot(T_2s, xs) / (m * np.sum(T_2s))  print(f'step:{i}, tau={tau}, p1={p_1:.2f}, p2={p_2:.2f}, log_likelihood={sum(lls):.2f}')  theta = [tau, p_1, p_2]     # stop when likelihood doesn't improve anymore  if abs(sum(lls) - last_ll) < 0.001:      break  else:      last_ll=sum(lls)

https://github.com/financialnoob/misc/blob/305bf8bc7cbdddaf47d40078100ba27935ff4452/6.introduction_to_em_algorithm.ipynb

編輯：于騰凱
校對(duì)：林亦霖