在结束了有参估计，无参估计后，现在记录混合模型（Mixture models）。这里附一张有参和无参的对比图（本来应该附在Part 2的，不想回去改了。。）：
字面意思，混合模型就是有参模型和无参模型的混合。
举个例子，高斯模型的混合（Mixture of Gaussians，MoG）。
现有三个高斯模型如下：
我们可以将其视为：
其概率密度可以近似表示为：
p
(
x
)
=
∑
j
=
1
M
p
(
x
∣
j
)
p
(
j
)
p(x)=\sum^M_{j=1}p(x|j)p(j)
p(x)=j=1∑Mp(x∣j)p(j)各变量有如下关系：
大概解释一下：
p
(
x
∣
j
)
p(x|j)
p(x∣j)表示第
j
j
j个高斯模型的概率密度，
p
(
j
)
=
π
j
p(j)=\pi_j
p(j)=πj表示第
j
j
j个高斯模型出现的概率，是一种先验概率。
这里记住两点：
即是说，混合模型概率密度积分为1，然后对于混合模型的每个子模型，我们都要求出对应的三个值（均值，方差和模型先验），这是对于MoG而言的。
对于简单的（单个）高斯模型，我们可以用MLE近似求解，但在高斯混合中无法使用该方法，因为我们在求参数时并不知道这个参数属于高斯混合的哪个子分布。如下：
这个
p
(
j
∣
x
)
p(j|x)
p(j∣x)我们是无法观测到的，毕竟如果这个都知道，那高斯混合就只是多个高斯模型的简单叠加而已。
也就是说单个高斯模型的参数会取决于所有高斯模型的参数。所以MLE不可用。
分为两种：
简单说的话，hard assignments是根据Label划分的，也就是说每个数据点只属于一个label，比如说K-mean；soft assignment是根据概率区分的，每个数据点可能属于多个label，但概率不同。如下图：
这里只介绍soft assignment，毕竟两个几乎是一样的。这就要用到这章的重点了，
E
M
EM
EM算法。
E
M
EM
EM算法是一种迭代算法，本质上它也是一种最大似然估计的方法，其特点是，当我们的数据不完整时，比如说只有观测数据，缺乏隐含数据时，可以用EM算法进行迭代推导。
该算法包括
E
E
E步骤和
M
M
M步骤，这里不具体介绍一般性的
E
M
EM
EM算法，主要说明其在高斯混合中如何运用。其步骤大体可以概括为：
现在说一下更加详细的步骤。
假设现在有两组数据，即是观测数据和隐藏数据：
组合后形参完整数据：
联合密度为
p
(
Z
)
=
p
(
X
,
Y
)
=
p
(
Y
∣
X
)
p
(
X
)
p(Z)=p(X,Y)=p(Y|X)p(X)
p(Z)=p(X,Y)=p(Y∣X)p(X)即是
p
(
Z
∣
θ
)
=
p
(
X
,
Y
∣
θ
)
=
p
(
Y
∣
X
,
θ
)
p
(
X
∣
θ
)
p(Z|\theta)=p(X,Y|\theta)=p(Y|X,\theta)p(X|\theta)
p(Z∣θ)=p(X,Y∣θ)=p(Y∣X,θ)p(X∣θ)在高斯混合中：
p
(
X
∣
θ
)
p(X|\theta)
p(X∣θ)是混合模型的似然
p
(
Y
∣
X
,
θ
)
p(Y|X,\theta)
p(Y∣X,θ)是混合模型中子分布的估计
对于不完整的数据（观测数据），其似然为：
L
(
θ
∣
X
)
=
p
(
X
∣
θ
)
=
∏
n
=
1
N
p
(
X
n
∣
θ
)
L(\theta|X)=p(X|\theta)=\prod_{n=1}^{N}p(X_n|\theta)
L(θ∣X)=p(X∣θ)=n=1∏Np(Xn∣θ)对于完整数据（Z），其似然为：
在这里我们虽然不知道
Y
Y
Y，但如果我们知道当前的参数猜测
θ
i
−
1
\theta^{i-1}
θi−1，我们就能用它来预测
Y
Y
Y。
在这里我们计算完整数据的对数似然的期望，如下：
其中
X
X
X和
θ
i
−
1
\theta^{i-1}
θi−1是已知的。更进一步展开如下：
这个等式是根据均值和积分的关系写出来的。即是如下关系：
E
[
x
]
=
∫
X
x
f
(
x
)
d
x
E[x]=\int_Xxf(x)dx
E[x]=∫Xxf(x)dx其中
f
(
x
)
f(x)
f(x)是概率密度（这部分不太确定，有错的话麻烦大家指出）。
我们需要最大化这个
Q
Q
Q函数。
接下来是
E
M
EM
EM算法：
E-step(expectation): 计算
p
(
y
∣
X
,
θ
i
−
1
)
p(y|X,\theta^{i-1})
p(y∣X,θi−1)以便计算
Q
(
θ
,
θ
i
−
1
)
Q(\theta,\theta^{i-1})
Q(θ,θi−1);
M-step(maximization): 最大化
Q
Q
Q函数求出
θ
\theta
θ
θ
^
=
a
r
g
m
a
x
θ
Q
(
θ
,
θ
i
−
1
)
\hat{\theta}=arg max_\theta Q(\theta,\theta^{i-1})
θ^=argmaxθQ(θ,θi−1)
这是一种迭代运算，我们要确保每次迭代中，第
i
i
i次的结果至少和第
i
−
1
i-1
i−1的一样好，即是：
Q
(
θ
i
,
θ
i
−
1
)
≥
Q
(
θ
i
−
1
,
θ
i
−
1
)
Q(\theta^i,\theta^{i-1})\geq Q(\theta^{i-1},\theta^{i-1})
Q(θi,θi−1)≥Q(θi−1,θi−1)
若该期望值对于
θ
\theta
θ来说是最大的，则可以认为(这部分未理解)：
L
(
θ
i
∣
X
)
≥
L
(
θ
i
−
1
∣
X
)
L(\theta^i|X)\geq L(\theta^{i-1}|X)
L(θi∣X)≥L(θi−1∣X)
也就是说，在每次迭代中观测数据的对数似然都会不断增大（或者至少保持不变），最终达到局部最大值。
所以，在实际运用中，初始化对于
E
M
EM
EM算法很重要，一个不好的初始化可能会使结果停在一个不好的局部最优值中。
步骤：
至此，高斯混合的
E
M
EM
EM算法就结束了，然后还有最后一个问题，这部分不太理解，但把结论贴上来吧，以后再探究：
也就是给出数据点，然后用EM算法进行高斯拟合：
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
def load_data(path):
return np.loadtxt(path)
def init_params(k, d = 2):
pi = np.ones(k) * 1/k
mu = [np.random.rand(2, 1) for i in range(k)]
cov = [np.eye(2) for i in range(k)]
return pi, mu, cov
def EM(iter_times, data, k):
pi, mu, cov = init_params(k)
N = data.shape[0]
alpha = np.zeros((N, k))
likelihood_list = []
def update_alpha(alpha, N, k):
for i in range(N):
for j in range(k):
p_ij = pi[j] * multivariate_normal.pdf(data[i], mu[j].flatten(), cov[j])
alpha[i][j] = p_ij
likelihood[i] += p_ij
alpha[i] = alpha[i] / np.sum(alpha[i])
return alpha
sigma = np.empty((N, 2, 2))
for i in range(N):
d = data[i, :].reshape(2, -1)
sigma[i] = np.dot(d, d.T)
for i in range(iter_times):
likelihood = np.zeros(N)
# E-step: update alpha
alpha = update_alpha(alpha, N, k)
likelihood_list.append(np.sum(np.log(likelihood)))
# M-step: update mu, cov, pi
N_j = np.sum(alpha, axis = 0)
for j in range(k):
# update mu
alpha_x = 0
alpha_x_mu = 0
for n in range(N):
alpha_x += (alpha[n][j] * data[n])
alpha_x = alpha_x.reshape(2,1)
mu[j] = alpha_x / N_j[j]
# update pi
pi[j] = N_j[j] / N
# update cov
for n in range(N):
alpha_x_mu += alpha[n][j]*(data[n] - mu[j].T)*(data[n] - mu[j].T).T
cov[j] = alpha_x_mu / N_j[j]
plot_contour(iter_times, data, mu, cov, k)
if iter_times == 30:
plt.figure()
plt.plot(np.arange(1, 31), np.array(likelihood_list))
plt.title('log-likelihood for every iteration')
plt.xlabel('iteration')
plt.ylabel('log-likelihood')
plt.grid()
plt.show()
def plot_contour(iter_num, data, mu, cov, k):
plt.figure()
plt.title("iter_num = %d" % iter_num)
plt.scatter(data[:,0], data[:,1])
x_min, x_max = plt.gca().get_xlim()
y_min, y_max = plt.gca().get_ylim()
num = 50
x = np.linspace(x_min, x_max, num)
y = np.linspace(y_min, y_max, num)
X, Y = np.meshgrid(x, y)
for sub_k in range(k):
Z = np.zeros_like(X)
for i in range(len(x)):
for j in range(len(y)):
z = multivariate_normal.pdf(np.array([[x[i]], [y[j]]]).flatten(), mu[sub_k].flatten(), cov[sub_k])
Z[j, i] = z
plt.contour(X, Y, Z)
plt.show()
path = "./dataSets/gmm.txt"
def main():
k = 4
iter_lst = [1, 3, 5, 10, 30]
pi, mu, cov = init_params(k)
data = load_data(path)
for i in iter_lst:
EM(i, data, k)
if __name__ == '__main__':
main()
效果如下图所示：
迭代1次：
迭代3次：
迭代5次：
迭代10次：
迭代30次：
每次迭代的对数似然：
