《统计学习方法（第二版）》-李航-第九章EM算法及其推广 思维导图笔记及课后练习题目答案（使用python3编写EM算法及在高斯混合模型中的应用）

第九章思维导图总结：

import numpy as np

9.1

y = np.array([[1, 1, 0, 1, 0, 0, 1, 0, 1, 1]]).T
m = y.shape[0]
theta = np.array([[0.46, 0.55, 0.67]]).T  # initialization

for i in range(100):
    theta_old = theta.copy()
    q_theta = (theta[0] * theta[1] ** y * (1 - theta[1]) ** (1 - y)) / (theta[0] * theta[1] ** y * (1 - theta[1]) ** (1 - y) + 
        ((1 - theta[0]) * theta[2] ** y * (1 - theta[2]) ** (1 - y)))
    theta[0] = q_theta.sum() / m
    theta[1] = (q_theta.T @ y) / q_theta.sum()
    theta[2] = ((1 - q_theta.T) @ y) / (1 - q_theta).sum()
    print(f'iteration {i} with theta:',theta.T)
    if np.linalg.norm(theta - theta_old) < 0.000001:
        print(f'converge at iteragtion: {i}')
        break
print(f'final theta: {theta.T}')

iteration 0 with theta: [[0.46186284 0.534595   0.65613464]]
iteration 1 with theta: [[0.46186284 0.534595   0.65613464]]
converge at iteragtion: 1
final theta: [[0.46186284 0.534595   0.65613464]]

9.2

$$\begin{array}{l}F(\tilde{P},\theta )\\ =E_p[\log P(Y,Z\mid \theta )]+H(\tilde{P})\\ =E_p[\log P(Y,Z\mid \theta )]-E_{\tilde{P}}\log \tilde{P}(Z)\\ ={\sum }_z{\tilde{P}}_{\theta }(Z)\log P(Y,Z\mid \theta )-{\sum }_Z\tilde{P}(Z)\log \tilde{P}(Z)\\ ={\sum }_{z}P(Z\mid Y,\theta )\log P(Y,Z\mid \theta )-{\sum }_{Z}P(Z\mid Y,\theta )\log P(Z\mid Y,\theta )\\ ={\sum }_{z}P(Z\mid Y,\theta )\log \frac{P(Y,Z\mid \theta )}{P(Z\mid Y,\theta )}\\ ={\sum }_{z}P(Z\mid Y,\theta )\log P(Y\mid \theta )\\ =\log P(Y\mid \theta )\end{array}$$

9.3

from scipy.stats import norm as Gaussian


y = np.array([[-67, -48, 6, 8, 14, 16, 23, 24, 28, 29, 41, 49, 56, 60, 75]]).T
m = y.shape[0]
k = 2  # number of sub_distribution
alpha = np.array([1.0 / k for _ in range(k)]).T # initialize parameters
miu = np.array([[y.mean() for _ in range(k)]]).T
sigma = np.array([[np.sqrt(y.var()) for _ in range(k)]]).T


for i in range(100):
    alpha_old = alpha.copy()
    miu_old = miu.copy()
    sigma_old = sigma.copy()
    gamma_hat = np.zeros((m, k))
    for j in range(k):
        gamma_hat[:, j, None] = alpha[j] * Gaussian.pdf(y, miu[j], sigma[j]**2)
    gamma_hat = gamma_hat / gamma_hat.sum(axis=1).reshape((-1, 1))
    miu = (gamma_hat.T @ y) / gamma_hat.sum(axis=0).T
    sigma = (gamma_hat * (y - miu_old.T)**2).sum(axis=0).T / gamma_hat.sum(axis=0).T
    sigma = np.sqrt(sigma)
    alpha = gamma_hat.sum(axis=0).T / m
    if np.sum((miu - miu_old)**2 + (sigma - sigma_old)**2 + (alpha - alpha_old)**2) < 1e-15:
        print(f'converge at iteragtion: {i}')
        break
print(f'final alpha: {alpha.T} \nfinal miu: {miu.T} \nfinal sigma:{sigma.T}')

converge at iteragtion: 0
final alpha: [0.5 0.5] 
final miu: [[20.93333333 20.93333333]
 [20.93333333 20.93333333]] 
final sigma:[36.46453376 36.46453376]

看一下10个子模型（高斯分布）的结果：

from scipy.stats import norm as Gaussian


y = np.array([[-67, -48, 6, 8, 14, 16, 23, 24, 28, 29, 41, 49, 56, 60, 75]]).T
m = y.shape[0]
k = 10  # number of sub_distribution
alpha = np.array([1.0 / k for _ in range(k)]).T # initialize parameters
miu = np.array([[y.mean() for _ in range(k)]]).T
sigma = np.array([[np.sqrt(y.var()) for _ in range(k)]]).T


for i in range(100):
    alpha_old = alpha.copy()
    miu_old = miu.copy()
    sigma_old = sigma.copy()
    gamma_hat = np.zeros((m, k))
    for j in range(k):
        gamma_hat[:, j, None] = alpha[j] * Gaussian.pdf(y, miu[j], sigma[j]**2)
    gamma_hat = gamma_hat / gamma_hat.sum(axis=1).reshape((-1, 1))
    miu = (gamma_hat.T @ y) / gamma_hat.sum(axis=0).T
    sigma = (gamma_hat * (y - miu_old.T)**2).sum(axis=0).T / gamma_hat.sum(axis=0).T
    sigma = np.sqrt(sigma)
    alpha = gamma_hat.sum(axis=0).T / m
    if np.sum((miu - miu_old)**2 + (sigma - sigma_old)**2 + (alpha - alpha_old)**2) < 1e-15:
        print(f'converge at iteragtion: {i}')
        break
print(f'final alpha: {alpha.T} \nfinal miu: {miu.T} \nfinal sigma:{sigma.T}')

converge at iteragtion: 0
final alpha: [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1] 
final miu: [[20.93333333 20.93333333 20.93333333 20.93333333 20.93333333 20.93333333
  20.93333333 20.93333333 20.93333333 20.93333333]
 [20.93333333 20.93333333 20.93333333 20.93333333 20.93333333 20.93333333
  20.93333333 20.93333333 20.93333333 20.93333333]
 [20.93333333 20.93333333 20.93333333 20.93333333 20.93333333 20.93333333
  20.93333333 20.93333333 20.93333333 20.93333333]
 [20.93333333 20.93333333 20.93333333 20.93333333 20.93333333 20.93333333
  20.93333333 20.93333333 20.93333333 20.93333333]
 [20.93333333 20.93333333 20.93333333 20.93333333 20.93333333 20.93333333
  20.93333333 20.93333333 20.93333333 20.93333333]
 [20.93333333 20.93333333 20.93333333 20.93333333 20.93333333 20.93333333
  20.93333333 20.93333333 20.93333333 20.93333333]
 [20.93333333 20.93333333 20.93333333 20.93333333 20.93333333 20.93333333
  20.93333333 20.93333333 20.93333333 20.93333333]
 [20.93333333 20.93333333 20.93333333 20.93333333 20.93333333 20.93333333
  20.93333333 20.93333333 20.93333333 20.93333333]
 [20.93333333 20.93333333 20.93333333 20.93333333 20.93333333 20.93333333
  20.93333333 20.93333333 20.93333333 20.93333333]
 [20.93333333 20.93333333 20.93333333 20.93333333 20.93333333 20.93333333
  20.93333333 20.93333333 20.93333333 20.93333333]] 
final sigma:[36.46453376 36.46453376 36.46453376 36.46453376 36.46453376 36.46453376
 36.46453376 36.46453376 36.46453376 36.46453376]

9.4

还没看非监督呢