统计学习方法 第九章习题答案

习题9.1

如例9.1的三硬币模型.假设观测数据不变，试选择不同的初值，例如${\pi }^0=0.46,p^0=0.55,q^0=0.67$求模型参数$\theta =(\pi ,p,q)$的极大似然估计。
例9.1（三硬币模型）假设有3枚硬币，分别记作A，B，C.这些硬币正面出现的概率分别是π，p和q.进行如下掷硬币试验：先掷硬币A，根据其结果选出硬币B或硬币C，正面选硬币B，反面选硬币C；然后掷选出的硬币，掷硬币的结果，出现正面记作1，出现反面记作0；独立地重复n次试验（这里，n=10），观测结果如下：
1，1，0，1，0，0，1，0，1，1假设只能观测到掷硬币的结果，不能观测掷硬币的过程.问如何估计三硬币正面出现的概率，即三硬币模型的参数.
解答

由公式9.5得到
${\mu }_j^1=\{\begin{array}{c}0.537y_j=0\\ 0.412y_j=1\end{array}$
利用迭代公式9.6~9.8迭代
${\pi }^1=0.462,p^1=0.535,q^1=0.656$
同理
${\mu }_j^2=\{\begin{array}{c}0.537y_j=0\\ 0.412y_j=1\end{array}$
${\pi }^2=0.462,p^2=0.535,q^2=0.656$
得到模型参数$\theta$的极大似然估计
$\hat{\pi }=0.462,\hat{p}=0.535,\hat{q}=0.656$
代码验证一下

import numpy as np
theta = np.array([0.46, 0.55, 0.67])
Y = np.array([1,1,0,1,0,0,1,0,1,1])
mu = np.zeros(Y.shape[0])
print("初值为pi=%f,p=%f,q=%f"%(theta[0], theta[1], theta[2]))
for j in range(2):
    print("第%d次迭代"%(j+1))
    # E step
    for i in range(Y.shape[0]):
        temp = theta[0]*theta[1]**Y[i]*(1-theta[1])**(1-Y[i])
        temp1 = (1 - theta[0])*theta[2]**Y[i]*(1-theta[2])**(1-Y[i])
        mu[i] = temp / (temp + temp1)
    print("mu=",mu)
    # M step
    theta[0] = 1.0 / Y.shape[0] * sum(mu)
    theta[1] = sum(mu * Y) / sum(mu)
    theta[2] = sum((1 - mu)*Y) / sum(1-mu)
    print("pi=%f,p=%f,q=%f" % (theta[0], theta[1], theta[2]))

执行结果：

习题9.2

证明引理9.2.
若${\tilde{P}}_{\theta }(Z)=P(Z|Y,\theta )$，则$F(\tilde{P},\theta )=logP(Y|\theta )$
解答
$$\begin{gathered} F(\tilde{P},\theta ) \\ =E_p[logP(Y,Z|\theta )]+H(\tilde{P}) \\ =E_p[logP(Y,Z|\theta )]-E_{\tilde{P}}log\tilde{P}(Z) \\ ={\sum }_z{\tilde{P}}_{\theta }(Z)logP(Y,Z|\theta )-{\sum }_Z\tilde{P}(Z)log\tilde{P}(Z) \\ ={\sum }_{z}P(Z|Y,\theta )logP(Y,Z|\theta )-{\sum }_{Z}P(Z|Y,\theta )logP(Z|Y,\theta ) \\ ={\sum }_{z}P(Z|Y,\theta )log\frac{P(Y,Z|\theta )}{P(Z|Y,\theta )} \\ ={\sum }_{z}P(Z|Y,\theta )logP(Y|\theta ) \\ =logP(Y|\theta ) \end{gathered}$$

习题9.3

已知观测数据
-67，-48，6，8，14，16，23，24，28，29，41，49，56，60，75
试估计两个分量的高斯混合模型的5个参数.

其中

解答
先简单理一下题目的意思，这个模型由两个高斯模型混合而成，5个参数指的是两个高斯模型的参数${\mu }_0,{\sigma }_0,{\mu }_1,{\sigma }_1$，另外还有其对应的系数${\alpha }_0,{\alpha }_1({\alpha }_0+{\alpha }_1=1)$，相当于只需要求一个系数。
初值设置为$\sigma =1.0,\mu =0.5,\alpha =0.5$

import numpy as np
from scipy.stats import norm

y = np.array([-67, -48, 6, 8, 14, 16, 23, 24, 28, 29, 41, 49, 56, 60, 75])
K = 2  # 两个高斯
N = 15  # y有15个数据

# 参数初始化
mu = np.array([0.5, 0.5])
sigma = np.array([1.0, 1.0]) * 10
alpha = np.array([0.5, 0.5])

for i in range(10):
    gm = np.zeros((N, K))

    # E 步
    for j in range(N):
        for k in range(K):
            gm[j, k] = alpha[k] * norm(mu[k], sigma[k]).pdf(y[j]) #使用scipy实现高斯分布
        gm[j, :] /= sum(gm[j, :])  # gm[j,:] = gm[j,:] /sum(gm[j,:])

    # M 步
    mu2 = y.dot(gm) / sum(gm)
    alpha2 = sum(gm) / N
    sigma2 = np.zeros((2,))
    sigma2[0] = sum(gm[:, 0] * (y - mu[0]) ** 2) / sum(gm[:, 0])
    sigma2[1] = sum(gm[:, 1] * (y - mu[1]) ** 2) / sum(gm[:, 1])
    #判断是否收敛
    if sum((mu - mu2) ** 2 + (sigma - sigma2) ** 2 + (alpha - alpha2) ** 2) < 0.01:
        break

    mu = mu2
    sigma = sigma2
    alpha = alpha2
    print("第%d次迭代\nalpha_0=%f,mu_0=%f,sigma_0=%f\nalpha_1=%f,mu_1=%f,sigma_1=%f\n"
          %(i+1, alpha[0], mu[0], sigma[0], alpha[1], mu[1], sigma[1]))

运行结果

习题9.4

EM算法可以用到朴素贝叶斯法的非监督学习.试写出其算法.
参考Blog