高斯混合模型

郝鸿涛 / 2024-11-21

混合分布 #

假设一所艺术类高校测量了一千名学生的身高。男生占 40%,其平均身高为 175cm,标准差为 7cm。女生占 60%,平均身高为 162cm,标准差为 6cm。

把数据合并,我们看看最后的直方分布:

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.stats import norm
import json

np.random.seed(42)
Show Code

n_female = 600
mu_female = 162
sd_female = 6
n_male = 1000 - n_female
mu_male = 175
sd_male = 7

female_heights = np.random.normal(mu_female, sd_female, n_female)
male_heights = np.random.normal(mu_male, sd_male, n_male)
heights = np.concatenate([female_heights, male_heights])
np.random.shuffle(heights)

plt.figure(figsize=(12, 6))
plt.hist(heights, bins = 30, density=True, alpha=0.7, color='gray')
plt.title('Histogram of Heights')
plt.xlabel("Height (cm)")
plt.ylabel("Density")
plt.show()

png

我们看到,上面的分布大概有两个峰顶,一个在女生平均身高,另一个在男生平均身高。

混合在一起的身高分布可以用「高斯混合模型」(Gaussian Mixture Model) 来描述:

$$0.4 * \mathcal{N}(175, 7^2) + 0.6 * \mathcal{N}(162, 6^2)$$

Show Code

x = np.linspace(140, 200, 100)
female_dist = norm.pdf(x, mu_female, sd_female)
male_dist = norm.pdf(x, mu_male, sd_male)
plt.plot(x, female_dist, 'r--', label="Female (60%)")
plt.plot(x, male_dist, 'b--', label="Male (40%)")
plt.plot(x, 0.6*female_dist + 0.4*male_dist, 'k-', label="Mixture")
plt.xlabel("Height (cm)")
plt.ylabel("Density")
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

png

E-M 算法 #

假设现在我只把一千名学生的身高数据给你了,但是不知道男、女生各自的比例,也不知道各组的分布参数(平均值、标准差)。请问,如何把男女比例以及各组的分布参数求出来?

我们的思路大概是这样:

def em_algo(heights, initial_means, initial_stds, n_iter = 100):
    mean1, mean2 = initial_means 
    std1, std2 = initial_stds
    pi1 = pi2 = 0.5

    history = []

    history.append({
            'iteration': 0,
            'mean1': mean1,
            'mean2': mean2,
            'std1': std1,
            'std2': std2,
            'pi1': pi1,
            'pi2': pi2
        })

    for iter in range(n_iter):
        prob1 = norm.pdf(heights, mean1, std1)*pi1
        prob2 = norm.pdf(heights, mean2, std2)*pi2

        # responsibility
        resp1 = prob1/(prob1+prob2)
        resp2 = prob2/(prob1+prob2)

        pi1 = np.mean(resp1)
        pi2 = np.mean(resp2)

        mean1 = np.sum(resp1*heights)/np.sum(resp1)
        mean2 = np.sum(resp2*heights)/np.sum(resp2)

        std1 = np.sqrt(np.sum(resp1 * (heights - mean1)**2) / np.sum(resp1))
        std2 = np.sqrt(np.sum(resp2 * (heights - mean2)**2) / np.sum(resp2))

        history.append({
            'iteration': iter + 1,
            'mean1': mean1,
            'mean2': mean2,
            'std1': std1,
            'std2': std2,
            'pi1': pi1,
            'pi2': pi2
        })

    return history
res = em_algo(
    heights, 
    initial_means = [173, 165], 
    initial_stds = [5, 4.5], 
    n_iter = 100
)
res[-1]
{'iteration': 100,
 'mean1': 174.98737549406312,
 'mean2': 161.59379516925364,
 'std1': 7.155329609964584,
 'std2': 5.570256355396785,
 'pi1': 0.42928333161106164,
 'pi2': 0.5707166683889384}

我们看到,这个已经很接近了。让我们把 iteration 增加一下:

res = em_algo(
    heights, 
    initial_means = [173, 165], 
    initial_stds = [5, 4.5], 
    n_iter = 1000
)
res[-1]
{'iteration': 1000,
 'mean1': 174.89513198090734,
 'mean2': 161.5536415719893,
 'std1': 7.193088560338676,
 'std2': 5.553978958610935,
 'pi1': 0.4339690847005224,
 'pi2': 0.5660309152994776}

好像并没有好到哪里去。

那我们想,如果初始值非常不靠谱呢?结果会怎么样?

res = em_algo(
    heights, 
    initial_means = [199, 177], 
    initial_stds = [10, 9], 
    n_iter = 100
)
res[-1]
{'iteration': 100,
 'mean1': 176.84518000880422,
 'mean2': 162.49234107631176,
 'std1': 6.421576228504393,
 'std2': 5.971350004395543,
 'pi1': 0.33798852615801583,
 'pi2': 0.6620114738419841}

我们看到,结果依然不错。

E-M 算法背后的数学原理 #

E-M 算法很好理解,但是我们怎么知道随着迭代次数的增加,结果越来越好呢?这就需要数学证明。

这部分我只懂了一部分而已。下面我把我懂的部分解释一下。

对于一个单独的身高 $x_i$,其发生的概率是:

$$P(x_i) = \pi_1 \cdot \mathcal{N}(x_i | \mu_1, \sigma_1^2) + \pi_2 \cdot \mathcal{N}(x_i | \mu_1, \sigma_2^2)$$

我们引入一个「隐变量」$z$,表示 $x_i$ 来自哪个分布。上面的公式等价于

$$P(x_i | \theta) = \sum_{z} P(x_i | z, \theta) P(z | \theta)$$

这里需要理解的是 $z$$\pi_1, \pi_2$ 不同,因为 $z$ 是对于每一个身高来说的,而$\pi_1, \pi_2$是对于所有身高来说的。

因为每个身高数据是独立的,所以 $P(X|\theta)$ 是一个联合概率:

$$P(x_i | \theta) = \prod_{i=1}^{n} \sum_{z} P(x_i | z, \theta) P(z | \theta)$$

为了便于计算,也为了防止数字下溢,我们取对数,这样就把连乘转化为加法:

$$\log P(x_i | \theta) = \sum_{i=1}^n \log \left(\sum_{z} P(x_i | z, \theta) P(z | \theta)\right)$$

我们把上面的公式称为「对数似然」。我们的目标是找到一组 $\theta$ 让对数似然最大化,这也代表着概率最大,也就是说当参数为 $\theta$ 是,这些身高数据被我们观测到的概率最大。

但是这个优化起来,以我们现有的技术,是非常困难的。因为 $\log\left(\sum_z \right)$ 没办法直接优化。另外,$z$ 是未知的,也让这个求解非常困难。

剩下的我就不太懂了。之后再学吧。

#统计

最后一次修改于 2024-11-21