使用 EM 算法求混合高斯正态分布

使用 EM 算法求混合高斯正态分布

1.问题的提出

假设学校男女生的身高服从不同的正态分布，男生的身高分布服从 $N({\mu }_1,{\sigma }_1^2)$ ， 女生的身高分布服从 $N({\mu }_2,{\sigma }_2^2)$ 。当使用样本数据估计总体参数 $\theta =({\mu }_1,{\mu }_2,{\sigma }_1^2,{\sigma }_2^2{)}^T$ 时，但是不知道每个样本是男生身高数据还是女生的身高数据，那么便不能使用通常的极大似然估计方法去估计总体参数了。

2.EM（Expectation Maximization） 算法

第一步，猜测男生和女生身体分布的参数 ${\theta }_0=({{\mu }_1}_0,{{\mu }_2}_0,{{\sigma }_1^2}_0,{{\sigma }_2^2}_0{)}^T$，然后计算出每个人更可能是属于第一个分布还是第二分布。比如一个样本身高值为1.90，那么这个样本可能更大概率来自男生总体。
第二步，根据计算出来的概率，使用最大似然估计求出 ${\theta }_1=({{\mu }_1}_1,{{\mu }_2}_1,{{\sigma }_1^2}_1,{{\sigma }_2^2}_1{)}^T$。然后重新第一步的过程，使用 ${\theta }_1=({{\mu }_1}_1,{{\mu }_2}_1,{{\sigma }_1^2}_1,{{\sigma }_2^2}_1{)}^T$ 去计算出每个人可能属于的分布，概率计算出来以后，继续第二步使用极大似然估计计算出 ${\theta }_2=({{\mu }_1}_2,{{\mu }_2}_2,{{\sigma }_1^2}_2,{{\sigma }_2^2}_2{)}^T$；往复循环直到计算最有的参数。

3.EM 算法的数学推导

样本集 { x 1 , x 2 , . . . , x N } \{x_1,x_2,...,x_N\} {x1​,x2​,...,xN​} 包含有 N 个样本， 假设每个样本可能的类别集为 Z = { z 1 , z 2 , . . . , z m } Z = \{z_1,z_2,...,z_m\} Z={z1​,z2​,...,zm​}，样本 $x_i$ 在每个类别的概率设为 $Q_i(z)$ , 其中 $z\in Z$ 且 ${\sum }_{z\in Z}{Q}_i=1$ 。那么似然函数的求法如下：
$$\begin{array}{rl}L(\theta )& =\sum _{i=1}^N\ln p(x_i;\theta )\\ & =\sum _{i=1}^N\ln \sum _{z\in Z}p(x_i,z；\theta )\\ & =\sum _{i=1}^N\ln \sum _{z\in Z}{Q}_i(z)\frac{p(x_i,z；\theta )}{Q_i(z)}\\ & \ge \sum _{i=1}^N\sum _{z\in Z}{Q}_i(z)\ln \frac{p(x_i,z；\theta )}{Q_i(z)}\end{array}$$
​上述不等式使用的是 Jensen 不等式 ：如果 $f$ 是凹函数， $x$ 是随机变量， 那么$E[f(X)]\le f(E[X])$ ，当且仅当 $X$ 是常量 $(E(X)=X)$ 时等式成立。要保证上述不等式成立，那么应有：
$$\begin{array}{rl}& \frac{p(x_i,z；\theta )}{Q_i(z)}=c\\ & ⟹p(x_i,z;\theta )=c{Q}_i(z)\\ & ⟹\sum _{z\in Z}p(x_i,z;\theta )=c\sum _{z\in Z}{Q}_i(z)\\ & ⟹\sum _{z\in Z}p(x_i,z;\theta )=c\end{array}$$
​从而
$$\begin{array}{rl}{Q}_i(z)& =\frac{p(x_i,z;\theta )}{{\sum }_{z\in Z}p(x_i,z;\theta )}\\ & =p(z|x_i;\theta )\end{array}$$
​因此我们的最大化目标转化为：
$$\max \sum _{i=1}^N\sum _{z\in Z}{Q}_i(z)\ln \frac{p(x_i,z；\theta )}{Q_i(z)}$$
​在计算混合正态时，假设总体服从 $N({\mu }_z,{\sigma }_z^2)$, 其中 $z=1,2$。那么求解的目标函数为:
$$\begin{array}{rl}& \max \sum _{i=1}^N\sum _{z=1}^2{Q}_i(z)\ln \frac{(2\pi {\sigma }_z^2{)}^{-1/2}\exp (-(x_i-{\mu }_z{)}^2/2{\sigma }_z^2)}{Q_i(z)}\\ & ⟹\sum _{i=1}^N\sum _{z=1}^2{Q}_i(z)[-\frac{1}{2}\ln (2\pi )-\frac{1}{2}\ln {\sigma }_z^2-\frac{(x_i-{\mu }_z{)}^2}{2{\sigma }_z^2}]\end{array}$$
​关于 ${\mu }_z$ 求导：
$$\begin{array}{rl}& \sum _{i=1}^N{Q}_i(z)(x_i-{\mu }_z)=0\\ & ⟹{\mu }_z=\frac{{\sum }_{i=1}^N{Q}_i(z)x_i}{{\sum }_{i=1}^2{Q}_i(z)}\end{array}$$
关于 ${\sigma }_z^2$ 求导：
$$\begin{array}{rl}& \sum _{i=1}^2{Q}_i(z)[-\frac{1}{2{\sigma }_z^2}+\frac{(x_i-{\mu }_z{)}^2}{2{\sigma }_z^4}]=0\\ & ⟹{\sigma }_z^2=\frac{{\sum }_{i=1}^N{Q}_i(z)(x_i-{\mu }_z{)}^2}{{\sum }_{i=1}^N{Q}_i(z)}\end{array}$$

4. Python 代码实现

import baostock as bs # 获取股票相关信息的库 
import datetime
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# 获取沪深300收盘价
if __name__ == '__main__':
    login_result = bs.login(user_id='anonymous', password='123456')
    stockcode = 'sh.000300'
    startdate = '1990-01-01'
    today = datetime.datetime.now()
    delta = datetime.timedelta(days=1)
    predate = today - delta
    strpredate = datetime.datetime.strftime(predate, '%Y-%m-%d')
    rs = bs.query_history_k_data("%s" % stockcode, "date,close", start_date=startdate, end_date=strpredate, frequency="d", adjustflag="3")
    # 获取沪深300从1990年迄今的收盘价格
    result_list=[]
    while (rs.error_code == '0') & rs.next():
        # 获取一条记录，将记录合并在一起
        result_list.append(rs.get_row_data())
        hushen_close=pd.DataFrame(result_list, columns = rs.fields)

# 计算对数收益率
np_array = np.array(hushen_close['close'], dtype=np.float64) # 读取数据框为 numpy array
log_np_array = np.log(np_array) # 去对数
log_return1 = np.diff(log_np_array) # 一阶差分，获取当期的对数收益率
log_return = 100*log_return1 # 百分数形式的对数收益率

# 绘制图形的一些设定
plt.rcParams['font.sans-serif']=['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False # 用来正常显示负号


# 用混合高斯正态分布拟合收益率的分布
# EM 算法
# 1.设定参数的初始值
mu1 = np.array([0.1])
mu2 = np.array([-0.05]) 
sigma1 = np.array([2])
sigma2 = np.array([1.2])

max_loop = 500 # 设定的最大迭代数
def gaussian(array, mu, sigma): # 定义正态分布的概率密度函数　
    return np.exp((-(array-mu)**2)/2*(sigma**2))/np.sqrt(2*np.pi*(sigma**2))
# 2.根据设定的参数值，去求解优化后的参数值
for i in range(max_loop):
    mu1_old = mu1.copy() 
    mu2_old = mu2.copy()
    p_1 = gaussian(log_return, mu1, sigma1)
    p_2 = gaussian(log_return, mu2, sigma2)
    # 求隐变量的概率分布，其中p1 + p2 = 1
    p1 = p_1/(p_1+p_2) # 样本属于第一种概率密度分布的概率
    p2 = p_2/(p_1+p_2) # 样本属于第二种概率密度分布的概率
    # 计算最大化目标函数的参数
    mu1 = np.dot(p1, log_return)/np.sum(p1)
    mu2 = np.dot(p2, log_return)/np.sum(p2)
    sigma1 = np.sqrt((np.sum(p1*(log_return-mu1)**2))/np.sum(p1))
    sigma2 = np.sqrt((np.sum((p2)*(log_return-mu2)**2))/np.sum(p2))
    # 设置收敛准则
    epsilon = 0.00001
    if ((np.abs(mu1-mu1_old)<epsilon).all()) and ((np.abs(mu2-mu2_old)<epsilon).all()):
        break
# 3.最优化的混合高斯分布概率
p_1_final = np.sum(p1)/p1.size
p_2_final = np.sum(p2)/p2.size

# 绘制图形
def mix_gaussian(array, mu1, sigma1, mu2, sigma2, p1, p2): # 定义混合正态分布的概率密度函数
    return p1*gaussian(array, mu1, sigma1)+p2*gaussian(array, mu2, sigma2)
plt.hist(log_return, 100, density=True,facecolor="blue", edgecolor="black", alpha=0.7)
x1 = np.linspace(-10,10,10000)
plt.plot(x1, mix_gaussian(x1, mu1, sigma1, mu2, sigma2, p_1_final, p_2_final),"r")  # 标准正态分布
plt.title('用混合正态分布拟合')
plt.show()

5. 绘制图形