記錄：EM 算法估計混合高斯模型參數

當概率模型依賴於無法觀測的隱性變量時，使用普通的極大似然估計法無法估計出概率模型中參數。此時需要利用優化的極大似然估計：EM算法。

在這里我只是想要使用這個EM算法估計混合高斯模型中的參數。由於直觀原因，采用一維高斯分布。

一維高斯分布的概率密度函數表示為：

多個高斯分布疊加在一起形成混合高斯分布：

其中：k 表示一共有 k 個子分布，。為什么累加之和為 1？因為哪怕是混合模型也表示一個概率密度，從負無窮到正無窮積分概率為 1，所以只有累加之和為 1才能保證，很簡單的推導。

設總體 ξ，總體服從混合高斯分布。 是一個取自總體的樣本。罷了，公式編輯實在慢到令人發指，簡單記錄而已，手寫。

以下是關於一維混合高斯分布的參數估計推導過程：

參考：周志華《機器學習》

簡單代碼實現一下，代碼很丑：

import numpy as np
import matplotlib.pyplot as plt


# 使用 numpy 生成兩組符合高斯分布（正態分布）的數據，然后將他們累加成混合模型，使用 EM 算法求解其中參數
# 假設兩個分布累加的系數 α1=0.6，α2=0.4
# 假設 N1 分布的均值 μ1=1.7，方差 δ1²=0.57²=0.3249
# 假設 N2 分布的均值 μ2=3.5，方差 δ2²=0.33²=0.1089
np.random.seed(77)
num1 = 6000
num2 = 4000
X1 = np.random.normal(1.7, 0.57, num1).astype(np.float32)
X2 = np.random.normal(3.5, 0.33, num2).astype(np.float32)
X = np.hstack((X1, X2))  # 其中包含兩個高斯分布的數據
np.random.shuffle(X)  # 混洗數據

re_tuple = plt.hist(X, 300, density=1, facecolor='r')
plt.show()

# 設置 EM 算法的初始值，任意設置
modulus = np.array([0.2, 0.8])
mean = np.array([1.1, 2.1])
var = np.array([1.2, 1.5])

# 首先計算每個樣本點由每一個獨立分布產生的概率，然后通過推導公式去更新參數
gamma_j_i = np.zeros((2, num1 + num2), dtype=np.float32)

# 設置迭代次數
epochs = 100
for epoch in range(epochs):
    print('開始第 %d 次迭代 ...' % (epoch + 1))
    # E 步
    part_1 = 1 / np.sqrt(2 * np.pi * var[0])
    part_2 = 1 / np.sqrt(2 * np.pi * var[1])
    for i in range(2):
        part_i = 1 / np.sqrt(2 * np.pi * var[i])
        for j in range(num1 + num2):
            p_m = (modulus[0] * (part_1 * np.exp(-1 * ((X[j] - mean[0]) ** 2) / (2 * var[0]))) +
                   modulus[1] * (part_2 * np.exp(-1 * ((X[j] - mean[1]) ** 2) / (2 * var[1]))))
            p_i = modulus[i] * (part_i * np.exp(-1 * ((X[j] - mean[i]) ** 2) / (2 * var[i])))
            gamma_j_i[i, j] = p_i / p_m

    # 中間計算步驟
    sum_gamma_j_i = np.sum(gamma_j_i, axis=1)
    sum_for_mean = np.matmul(gamma_j_i, X)
    sum_for_var = np.sum(gamma_j_i * np.square(np.broadcast_to(X, (2, num1 + num2)) - mean.reshape((2, 1))), axis=1)

    # M 步
    for i in range(2):
        mean[i] = sum_for_mean[i] / sum_gamma_j_i[i]
        modulus[i] = sum_gamma_j_i[i] / (num1 + num2)
        var[i] = sum_for_var[i] / sum_gamma_j_i[i]

    print('迭代 %d 次后得到的 N1 分布的比率、均值和方差分別為：%s %s %s' % (epoch + 1, modulus[0], mean[0], var[0]))
    print('迭代 %d 次后得到的 N2 分布的比率、均值和方差分別為：%s %s %s' % (epoch + 1, modulus[1], mean[1], var[1]))
    print()

# 迭代 100 次后得到的結果是：
# N1: 0.59798 1.69166 0.33037
# N2: 0.40202 3.49959 0.11023
# 總之，結果還不錯