正态分布的均值估计

模拟k=2个正态分布的均值估计。其中，ini_data(Sigma,Mu1,Mu2,k,N)函数用于生成训练样本,此训练样本是从两个正态分布中随机生成的,其中正态分布1的均值Mu1=55,方差Sigma=5;正态分布2的均值Mu2=22,方差为5;
由于本问题的视线无法直接从样本中获取两个正态分布参数,因此需要使用EM算法估计出具体的Mu1/Mu2的取值

from numpy import *
import numpy as np
import copy


# EM算法步骤1
def em_e(Sigma, k, N, P, Mu, X):
    for i in range(0, N):
        Denom = 0
        for j in range(0, k):
            Denom += np.exp(-1.0 / (2.0 * Sigma ** 2) * (X[i] - Mu[j]) ** 2)
        for j in range(0, k):
            Numer = np.exp(-1.0 / (2.0 * Sigma ** 2) * (X[i] - Mu[j]) ** 2)
            P[i, j] = Numer / Denom


# EM算法:步骤2
def em_m(k, N, P, X):
    for j in range(0, k):
        Numer = 0
        Denom = 0
        for i in range(0, N):
            Numer += P[i, j] * X[i]
            Denom += P[i, j]
        Mu[j] = Numer / Denom


Sigma = 5
Mu1 = 55
Mu2 = 22
k = 2
N = 1200
iter_num = 500
Epsilon = 0.0001
X = mat(zeros((N, 1)))
Mu = np.random.random(2)
P = np.zeros((N, k))
for i in range(0, N):
    if np.random.random(1) > 0.5:
        X[i] = np.random.normal() * Sigma + Mu1
    else:
        X[i] = np.random.normal() * Sigma + Mu2

print(f'初始<u1,u2>: {Mu}')
for i in range(iter_num):
    ole_Mu = copy.deepcopy(Mu)
    em_e(Sigma, k, N, P, Mu, X)
    em_m(k, N, P, X)
    print(i, Mu)
    if sum(abs(Mu - ole_Mu)) < Epsilon:
        break