机器学习实战EM算法python3实现

Import  math

Import  copy

Import  numpy  as  np

import  matplotlib.pyplot  as  plt

isdebug = False

# 指定k个高斯分布参数，这里指定k=2。注意2个高斯分布具有相同均方差Sigma，分别为Mu1,Mu2。

def ini_data(Sigma, Mu1, Mu2, k, N):

    global X  # X产生的数据 ,k维向量

    global Mu  # 初始均值

    global Expectations

    X = np.zeros((1, N))

    Mu = np.random.random(2)  # 随机产生一个初始均值。

    ExpectatiOns= np.zeros((N, k))  # k个高斯分布，100个二维向量组成的矩阵。

    for i in range(0, N):

        if np.random.random(1) > 0.5:

# 随机从均值为Mu1,Mu2的分布中取样。

            X[0, i] = np.random.normal() * Sigma + Mu1

        else:

            X[0, i] = np.random.normal() * Sigma + Mu2

        if isdebug:

            print("***********")

    print(u"初始观测数据X：")

    print(X)

# EM算法：步骤1，计算E[zij]

def e_step(Sigma, k, N):

# 求期望。sigma协方差，k高斯混合模型数，N数据个数。

    global Expectations  # N个k维向量

    global Mu

    global X

    for i in range(0, N):

        Denom = 0

    for j in range(0, k):

        Denom += math.exp((-1 / (2 * (float(Sigma ** 2)))) * (float(X[0, i] - Mu[j])) ** 2)

# Denom  分母项  Mu(j)第j个高斯分布的均值。

    for j in range(0, k):

        Numer = math.exp((-1 / (2 * (float(Sigma ** 2)))) * (float(X[0, i] - Mu[j])) ** 2)  # 分子项

    Expectations[i, j] = Numer / Denom  # 期望，计算出每一个高斯分布所占的期望，即该高斯分布以多大比例形成这个样本

# EM算法：步骤2，求最大化E[zij]的参数Mu

def m_step(k, N):

# 最大化

    global Expectations  # 期望值

    global X  # 数据

    for j in range(0, k):

# 遍历k个高斯混合模型数据

        Numer = 0  # 分子项

        Denom = 0  # 分母项

    for i in range(0, N):

        Numer += Expectations[i, j] * X[0, i]  # 每一个高斯分布的期望*该样本的值。

        Denom += Expectations[i, j]  # 第j个高斯分布的总期望值作为分母

        Mu[j] = Numer / Denom  # 第j个高斯分布新的均值，

# 算法迭代iter_num次，或达到精度Epsilon停止迭代

def run(Sigma, Mu1, Mu2, k, N, iter_num, Epsilon):

    ini_data(Sigma, Mu1, Mu2, k, N)

    print(X)

    print(u"初始:", Mu)  # 初始均值

    for i in range(iter_num):

        Old_Mu = copy.deepcopy(Mu)  # 算法之前的MU

        e_step(Sigma, k, N)

        m_step(k, N)

        print(i, Mu)  # 经过EM算法之后的MU，

        if sum(abs(Mu - Old_Mu))

            break

if __name__ == '__main__':

    run(6, 40, 20, 2, 100, 100, 0.001)

    plt.hist(X[0, :], 50)

    plt.show()

结果：

初始观测数据X：

[[ 34.78479349  40.8762506   12.9259573   17.16456673  16.3469559

   14.5410094   25.78720329  17.6779438   35.82670264  41.90071202

   50.52081174  26.38782308  21.04963529  14.56343039  38.21448941

   49.69993673  48.18442754  44.51268208  26.25675349  40.86309978

   44.02344275  41.49166553  21.55141972   5.75166566  34.16528679

   19.70743537  30.42653305  39.23467249  40.34117738  23.60869141

   28.5130148   17.287984    17.33970054  16.01856562  38.22256625

   34.81288304  15.24771646  42.2952032   49.09309169  53.95353144

   32.16791165  29.0789197   48.66267963  41.61831626  21.31745232

   35.886364    15.29341564  19.08278906  18.612558    50.07119199

   31.17027178  34.35832737  19.79392945  37.1648494   13.38859369

   15.28348077  21.86257621  28.53368242  18.00679612  29.3180723

   19.50236702  42.7221693   32.27174091  13.13910741  14.47212264

   29.98649409  46.30665712  32.14165658   2.84063103  25.67421301

   32.7497982   42.17391758  33.22993139  52.13921332  35.17490892

   45.82663241  41.61183305  20.31732587  15.17501376  22.13357363

   45.05877176  45.66095107  33.53623831  17.80499695  27.53976863

   24.60028367  27.87441091  18.43420556  30.05336577  16.23903851

   13.13443783  43.42828628  21.37174492  27.25854313  37.4277017 ]]

[[ 34.78479349  40.8762506   12.9259573   17.16456673  16.3469559

   14.5410094   25.78720329  17.6779438   35.82670264  41.90071202

   50.52081174  26.38782308  21.04963529  14.56343039  38.21448941

   49.69993673  48.18442754  44.51268208  26.25675349  40.86309978

   25.95215462  25.69322242  24.25158242  36.06655722  24.90611346

   44.02344275  41.49166553  21.55141972   5.75166566  34.16528679

   19.70743537  30.42653305  39.23467249  40.34117738  23.60869141

   28.5130148   17.287984    17.33970054  16.01856562  38.22256625

   34.81288304  15.24771646  42.2952032   49.09309169  53.95353144

   32.16791165  29.0789197   48.66267963  41.61831626  21.31745232

   35.886364    15.29341564  19.08278906  18.612558    50.07119199

   31.17027178  34.35832737  19.79392945  37.1648494   13.38859369

   15.28348077  21.86257621  28.53368242  18.00679612  29.3180723

   19.50236702  42.7221693   32.27174091  13.13910741  14.47212264

   29.98649409  46.30665712  32.14165658   2.84063103  25.67421301

   32.7497982   42.17391758  33.22993139  52.13921332  35.17490892

   45.82663241  41.61183305  20.31732587  15.17501376  22.13357363

   45.05877176  45.66095107  33.53623831  17.80499695  27.53976863

   24.60028367  27.87441091  18.43420556  30.05336577  16.23903851

   13.13443783  43.42828628  21.37174492  27.25854313  37.4277017 ]]

初始: [ 0.1492249   0.53474567]

0 [  0.1492249  37.4277017]

1 [  0.1492249  37.4277017]