EM算法例子简单理解（例题+基本思想+python实现

最新推荐文章于 2022-11-27 20:36:47 发布

小果一粒沙

最新推荐文章于 2022-11-27 20:36:47 发布

阅读量6.6k

点赞数 6

分类专栏： python 统计学文章标签：算法

本文链接：https://blog.csdn.net/qq_35167821/article/details/106167820

版权

python 同时被 2 个专栏收录

26 篇文章 3 订阅

订阅专栏

统计学

18 篇文章 1 订阅

订阅专栏

是研究生复试的时候问到了一个问题，我不会，导致复试成绩不好。复试完了，打算好好理解一下，于是有了下文（搞了差不多一个月也是很不好意思）：
为了不致使正文太过于臃肿，将讲述原理的啰嗦话放到另一篇文章中了。

题目

已知：
给定一个城市中成年男性和女性的身高分布， $h_0\thicksim N(175,6^2)$ , $h_1\thicksim N(165,5^2)$ ， $h_0, h_1$ 分别代表男性身高和女性身高，且男女的比例为5：4。
问题：
（1）从城市中任意取一个人，求他的身高分布。
（2）给定1000个人的身高资料，但是一个城市的身高比例是未知的，请利用样本的身高数据来估计男女的比例。

分析

一、第一问

这个被抽取到的人是不知道是男性还是女性的，想知道这个人的身高分布，需要同时考虑到男女身高混合加成。
这样，我们就可以使用条件概率的思想，加权计算就可以得到一个混合正态分布，它的曲线是一个双峰分布。
在这里插入图片描述
即便他们之间相互独立，分布函数这样加权计算之后得到的不是一个单峰的正态分布，而是一个双峰的混合正态分布。
其中他们的分布函数为：

$f(h)=\cfrac{5}{9} \cdot \dfrac{e^{-\cfrac{(h-175)^2}{2\cdot6^2}}}{\sqrt{2\pi}\cdot6} + \cfrac{4}{9} \cdot \dfrac{e^{-\cfrac{(h-165)^2}{2\cdot5^2}}}{\sqrt{2\pi}\cdot5}$

在这里插入图片描述
这里的双峰不太明显，主要是因为两个总体的均值相差太近了，调节一下均值参数即可，比如将均值变为180和160，就显现出来了。

二、第二问

男女比例不知道，仅知道又1000个人的身高资料，然而计算极大似然函数需要依赖于男女比例，这所以就无法直接使用极大似然函数来求得。这里的男还是女是一个隐含变量（只知道身高信息，不知道性别信息），无法求得极大似然函数，EM估计就是为估计这种存在隐含变量的参数而生的。

步骤：

给定参数的初始值，比如说比例男：女=1：1.以下的是 $\alpha_0(1)=0.5$ 和 $\alpha_1(1)=0.5$ ，其中 $\alpha_0(t)+\alpha_1(t)=1$ .
计算第 $t$ 步

EM步骤

给定数据
E步，对隐变量求期望，求得Q函数。
M步，在给定数据的条件下，求得新的迭代值。
重复2-3步，Q函数或者迭代值不再发生显著变化，停止迭代。

高斯混合模型参数估计的EM算法

Input: 观测数据 $y_1,y_2,\cdots,y_n$ 。
Output: 估计的参数 $\mu_k,\sigma_k,\alpha_k(k=1,2,\cdots,K).$

给定均值 $\mu_k$ 、标准差 $\sigma_k$ 、各类比重 $\alpha_k(k=1,2,\cdots, K)$ 的初值。
E步：根据当前模型参数，计算各分模型 $k$ 对观测数据 $y_j$ 的响应度。即 $\hat{\gamma}_{jk}=\cfrac{\alpha_k \varphi(y_j|\theta_k)}{\sum_{k=1}^{K} \alpha_k\varphi(y_j|\theta_k)}$
M步：计算新一轮的迭代的模型参数：
$\begin{aligned} \hat{\mu}_k & =\cfrac{\sum_{j=1}^{N}\hat{\gamma}_{jk} y_j}{\sum_{j=1}^{N}\hat{\gamma}_{jk}} \\ & \hat{\sigma_k}^2=\cfrac{\sum_{j=1}^{N}\hat{\gamma}_{jk} (y_j-\mu_k)^2}{\sum_{k=1}^{N}\hat{\gamma}_{jk}} \\ & \hat{\alpha}_k = \cfrac{\sum_{j=1}^{N}\hat{\gamma}_{jk}}{N} \end{aligned}$
重复2和3步，直至收敛。

其中，
$j(j=1,2,\cdots,N)$ 代表第几个样本，
$k(k=1,2,\cdots,K)$ 代表第几类。

代码

1-混合高斯分布的代码

R的代码

library(ggplot2) #画图的包
x <-seq(150,190,by=0.01) # 自变量，这里的height
a=5/9 # 男性比例
f_x<-a*exp(-(x-175)^2/(2*6^2))/(sqrt(2*pi)*6)+(1-a)*exp(-(x-165)^2/(2*5^2))/(sqrt(2*pi)*5) # 概率密度
df_x_fx<-data.frame(x,f_x) #数据框，方便画图
p<-ggplot(df_x_fx,aes(x,f_x))
p+geom_line() # 画连续的线

或者是python

# 参考别人自己修改的，但是忘记出处了，我看到一定贴上地址！！
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import math
 
#正态分布的概率密度函数。可以理解成 x 是 mu（均值）和 sigma（标准差）的函数
def normfun(x,mu,sigma):
    pdf = np.exp(-((x - mu)**2)/(2*sigma**2)) / (sigma * np.sqrt(2*np.pi))
    return pdf

mu_1=175
mu_2=165
sigma_1=6
sigma_2=5
# Python实现正态分布
# 绘制正态分布概率密度函数
x = np.linspace(150, 190, 1000)
y_sig = np.exp(-(x - mu_1) ** 2 /(2* sigma_1 **2))/(math.sqrt(2*math.pi)*sigma_1)\
+np.exp(-(x - mu_2) ** 2 /(2* sigma_2 **2))/(math.sqrt(2*math.pi)*sigma_2)
plt.plot(x, y_sig, "r-", linewidth=2)
plt.xlabel('height')
plt.ylabel('probability distribution function')
plt.grid(True)
plt.show()

2-男性女性身高比例估计

这是初步的小白写法，是比较简单的代码，还有一个三类模型的、比较正规的代码在后面一章。

EM_data.py

生成数据，分别放在两个文件中。后来看其实这样过于麻烦，其实混合放在一个文件中也可以。

'''
1. 给定数据，需要生成样本的正态分布参数
2. 使用numpy中的random方法生成数据
3. 存储数据，写进CSV文件中。
'''
import numpy as np
import csv

# -------------------1. 给定数据，需要生成样本的正态分布参数-------------------------
# 初始化试验数据。
# normal distribution with mu and sigma
mu = [165, 175]
sigma = [6, 5]
# sample_num_0, sample_num_1 = sample_num
sample_num = [2400, 1600]
N = sum(sample_num) # 样本容量

# -------------------2. 使用numpy中的random方法生成数据-------------------------
def create_normal(mu, sigma, size_in):
    """
    params: one-dimensional data.
    """
    temp = np.random.normal(loc=mu, scale=sigma, size=size_in)
    return temp # 返回的是一个列表
# 生成数据
def create_data(mu, sigma, size_in):
    """
    paras:
    all the parameters have the same dimension, 
    not one-dimensional,
    all of them are list.
    mu: the mean of the normal distribution.
    sigma: the standard deviation of the normal distribution.
    size_in: the number of each category.
    """
    # 生成服从正态分布的数据


    # 依次按照不同的分布，生成数据
    tag_data = []
    for i in range(len(sample_num)):
        tag_data.append(list(create_normal(mu[i],sigma[i],sample_num[i]))) # list(): np.ndarray -> list        
    return tag_data # 一个列表，里面的元素为列表。

male, female = create_data(mu,sigma,sample_num)
# print(type(male)) # 列表

# 数据存储地址：
path_male = r"D:\lagua\study\coding\PythonStudy\EM_data\male_1.csv"
path_female = r"D:\lagua\study\coding\PythonStudy\EM_data\female_1.csv"


# -------------------3. 存储数据，写进CSV文件中。-------------------------
# 将生成的数据放到CSV文件中。
def write_data(li, path_li):
    # 将以上生成的数据写入到CSV文件中。
    for li_i in li:
        with open(path_li, "w+", encoding="utf-8", newline='') as f:
            writer = csv.writer(f)
            to_be_write = []
            # 判断应该写入哪个文件。
            if path_li == path_male:
                tag_ = [0 for i in range(len(li))]
            elif path_li == path_female:
                tag_ = [1 for i in range(len(li))]
            else:
                print('The path that you have entered does not exist.')            
            # 将生成的数据，按照列表打包，写入到CSV文件中。
            for tag_i, li_h in zip(tag_, li):
                to_be_write.append([tag_i,li_h])
            writer.writerows(to_be_write)    
    print('The data has been wirtten to the csv file.')

# 执行函数，写入数据。
write_data(male,path_male)
write_data(female,path_female)

EM.py

读出数据，设置函数进行迭代，输出结果。

"""
目标-步骤分解：
1. 给出初始条件 完成。
2. 读出数据，完成。
3. 设计函数，E步，完成。
4. 设计函数，M步，完成。
5. 设计函数，终止条件，完成。-2020/6/10/10:00,功能实现。
TODO 使用类，在类中定义方法来写代码
- 主要随机生成
- 测定对初值的敏感性

- 已生成于CSV文件中，2020/6/5/15：43
"""
import numpy as np
import csv
import timeit

# ----------------------1、给定初始条件 ----------------------
# 正态总体的均值、方差
mu = [165, 175]
sigma = [6, 5]

# 两类的比例，分别是男性：女性
# 后面的默认参数，定义为不变对象--元组
alpha_init = (0.5, 0.5) # 实际数据是6：4

# 类别的个数
J = len(alpha_init)

# ----------------------2、读出数据 ----------------------
path_male = r"D:\lagua\study\coding\PythonStudy\EM_data\male_1.csv"
path_female = r"D:\lagua\study\coding\PythonStudy\EM_data\female_1.csv"
# 获取数据。
def get_data(path_in):
    height_data = []
    with open(path_in,'r',encoding='utf-8') as f:
        reader = csv.reader(f)
        for line in reader:
            # 从CSV文件中获取的，是字符串数据，需要自己转化为浮点型数据。
            height_data.append(float(line[1])) # 每一行是一个列表，从列表中取出数据。
    print('Height Data collected.')
    return height_data

# TODO 可以设置一个函数，直接获取数据，而不用分两步。
male_height = get_data(path_male)
female_height = get_data(path_female)
Mixed_height = male_height + female_height # 列表相加/合并

# 样本数据个数
N = len(Mixed_height)


# 返回一个数 正态分布函数的密度
def normal_pdf(x, mu_num, sigma_num):
    temp = np.exp(-((x-mu_num)**2)/(2*sigma_num**2)) / (sigma_num * np.sqrt(2*np.pi))
    return temp


def EucliDis(x, y):
    '''
    作用：计算二范数，即两个列表/向量之间的距离
    '''
    d = np.sqrt(sum([(argu1-argu2)**2 for argu1,argu2 in zip(x,y)]))
    return d
def alpha_update(gamma):
    '''
    作用：更新响应度
    原理：根据上一步的响应度，来计算新的迭代值，依照书上公式打的
    '''
    alpha_new = gamma.mean(axis=0)
    return alpha_new

# -------------------------EM 主要程序-------------------
def EM(mu, sigma, alpha=alpha_init, EPSILON=0.001):
    """
    mu：各总体的均值
    sigma：各总体的标准差
    alpha：隐含变量的初值
    EPSILON: 两次迭代，距离小于这个数，跳出循环
    """
    # 每个身高，每一个种类的概率密度函数值，直接根据身高数据全部算出来了，没有设定函数。
    PDF = np.zeros((N,J))
    for i,height in enumerate(Mixed_height):
        for j in range(J):
            PDF[i,j] = normal_pdf(height,mu[j],sigma[j])
    # ----------------------3、E步，求期望。 ----------------------
    # 不同类，不同参数，具体的概率密度函数
    ###########################没法将此函数声明在外面，因为其中有一个PDF需要计算。
    ##################TODO 将此函数声明在外面
    def gammajk(i, j, alpha):
        '''
        作用：计算分模型j对观测数据height[i]的响应度
        i: 第i个身高数据
        j: 假设这个身高数据所属的模型类别j
        '''
        p_multi = 0
        for jj in range(J):
            p_multi += alpha[jj]*PDF[i,jj]
        temp = alpha[j]*PDF[i,j]
        return temp/p_multi

    # 迭代，第t步和第t+1步 之间是怎么继承上一步的？
    alpha_cur = list(alpha_init) # 将之前的元组 --> 列表
    alpha_new = [0.4, 0.6] # 随便设置的跟上一步不同的数

    # ----------------------5、中止条件 ----------------------
    while 1:
        '''
        作用：设置迭代循环，找到最优的参数
        停止条件：两次迭代之间的距离小于EPSILON
        '''
        # 计算所有的响应度，方便以后计算，以及下次的更新
        gamma = np.zeros((N,J))
        for i in range(N):
            for j in range(J):
                gamma[i,j] = gammajk(i,j,alpha_cur)

        # ----------------------4、M步 ----------------------
        alpha_new = alpha_update(gamma) # 更新参数值

        # 判断条件
        if EucliDis(alpha_cur,alpha_new) < EPSILON:
            break
        alpha_cur = alpha_new # 将新的参数传到alpha_cur中，方便进行下次迭代
    
    return alpha_new

# 算法耗费的时间
# 算法的参数估计结果
t_start = timeit.default_timer()
print('last result is',EM(mu,sigma))
t_end = timeit.default_timer()

cost = t_end - t_start
print('Time cost of EM is %f'%cost)
# last result is [0.61507698 0.38492302]
# Time cost of EM is 0.194936

参考

LaTex 把上下标符号放在正上和正下方公式介绍
主要参考上面latex写公式下面的符号，结果还是当作下标写了，不过还是记录在这里，因为它也给了我帮助。
LaTex中把下标置于文本正下方的方法
将一个符号放在公式下面，如果不是数学符号，需要使用\mathop进行转换，即可使用\limits.
Python实现机器学习算法：EM算法
原来估计是可以有误差的，并且样本量越大，误差越小。

小果一粒沙

关注

6
点赞
踩
55

收藏

觉得还不错? 一键收藏
2
评论
EM算法例子简单理解（例题+基本思想+python实现

是研究生复试的时候问到了一个问题，我不会，导致复试成绩不好。复试完了，打算好好理解一下，于是有了下文：简介极大似然估计根据样本的信息，求得未知参数的估计。主要思想：样本信息出现的概率最大，也就是说我们想要找到一个参数，使得我们所有的样本出现的概率最大。步骤：写出关于似然函数。取对数（方便后来求极值运算）。对未知参数求偏导，令其为0（求极值），求得参数值为极大似然估计值。补充：当然，有的并不能够使用数学上求极值的方法来求得极大似然估计值，可以使用单调性求得。不管是哪种解法，都是要求找.
复制链接

扫一扫