2022 高教社杯(国赛数学建模)思路解析

2022高教社杯ABCD赛题思路解析：

https://blog.51cto.com/u_15734606/5574428
<br/>

算法背景

假设：

在这里插入图片描述

但现在我们让情况复杂一点，就是这 50 个男生和 50 个女生混在一起了。我们拥有 100 个人的身高数据，却不知道这 100 个人每一个是男生还是女生。

这时候情况就有点尴尬，因为通常来说，我们只有知道了精确的男女身高的正态分布参数我们才能知道每一个人更有可能是男生还是女生。但从另一方面去考量，我们只有知道了每个人是男生还是女生才能尽可能准确地估计男女各自身高的正态分布的参数。

这个时候有人就想到我们必须从某一点开始，并用迭代的办法去解决这个问题：我们先设定男生身高和女生身高分布的几个参数（初始值），然后根据这些参数去判断每一个样本（人）是男生还是女生，之后根据标注后的样本再反过来重新估计参数。之后再多次重复这个过程，直至稳定。这个算法也就是 EM 算法。

python算法实现

#! /usr/bin/env python
#! -*- coding=utf-8 -*-

#模拟两个正态分布的均值估计

from numpy import *
import numpy as np
import random
import copy

SIGMA = 6
EPS = 0.0001
#生成方差相同,均值不同的样本
def generate_data():    
    Miu1 = 20
    Miu2 = 40
    N = 1000
    X = mat(zeros((N,1)))
    for i in range(N):
        temp = random.uniform(0,1)
        if(temp > 0.5):
            X[i] = temp*SIGMA + Miu1
        else:
            X[i] = temp*SIGMA + Miu2
    return X

#EM算法
def my_EM(X):
    k = 2
    N = len(X)
    Miu = np.random.rand(k,1)
    Posterior = mat(zeros((N,2)))
    dominator = 0
    numerator = 0
    #先求后验概率
    for iter in range(1000):
        for i in range(N):
            dominator = 0
            for j in range(k):
                dominator = dominator + np.exp(-1.0/(2.0*SIGMA**2) * (X[i] - Miu[j])**2)
                #print dominator,-1/(2*SIGMA**2) * (X[i] - Miu[j])**2,2*SIGMA**2,(X[i] - Miu[j])**2
                #return
            for j in range(k):
                numerator = np.exp(-1.0/(2.0*SIGMA**2) * (X[i] - Miu[j])**2)
                Posterior[i,j] = numerator/dominator            
        oldMiu = copy.deepcopy(Miu)
        #最大化    
        for j in range(k):
            numerator = 0
            dominator = 0
            for i in range(N):
                numerator = numerator + Posterior[i,j] * X[i]
                dominator = dominator + Posterior[i,j]
            Miu[j] = numerator/dominator
        print (abs(Miu - oldMiu)).sum() 
            #print '\n'
        if (abs(Miu - oldMiu)).sum() < EPS:
            print Miu,iter
            break

if __name__ == '__main__':
    X = generate_data()
    my_EM(X)

<br/>