使用python如何實現(xiàn)一個em算法

'''
數(shù)據(jù)集：偽造數(shù)據(jù)集（兩個高斯分布混合）
數(shù)據(jù)集長度：1000
------------------------------
運行結(jié)果：
----------------------------
the Parameters set is:
alpha0:0.3, mu0:0.7, sigmod0:-2.0, alpha1:0.5, mu1:0.5, sigmod1:1.0
----------------------------
the Parameters predict is:
alpha0:0.4, mu0:0.6, sigmod0:-1.7, alpha1:0.7, mu1:0.7, sigmod1:0.9
----------------------------
'''

import numpy as np
import random
import math
import time

def loadData(mu0, sigma0, mu1, sigma1, alpha0, alpha1):
  '''
  初始化數(shù)據(jù)集
  這里通過服從高斯分布的隨機(jī)函數(shù)來偽造數(shù)據(jù)集
  :param mu0: 高斯0的均值
  :param sigma0: 高斯0的方差
  :param mu1: 高斯1的均值
  :param sigma1: 高斯1的方差
  :param alpha0: 高斯0的系數(shù)
  :param alpha1: 高斯1的系數(shù)
  :return: 混合了兩個高斯分布的數(shù)據(jù)
  '''
  # 定義數(shù)據(jù)集長度為1000
  length = 1000

  # 初始化第一個高斯分布，生成數(shù)據(jù)，數(shù)據(jù)長度為length * alpha系數(shù)，以此來
  # 滿足alpha的作用
  data0 = np.random.normal(mu0, sigma0, int(length * alpha0))
  # 第二個高斯分布的數(shù)據(jù)
  data1 = np.random.normal(mu1, sigma1, int(length * alpha1))

  # 初始化總數(shù)據(jù)集
  # 兩個高斯分布的數(shù)據(jù)混合后會放在該數(shù)據(jù)集中返回
  dataSet = []
  # 將第一個數(shù)據(jù)集的內(nèi)容添加進(jìn)去
  dataSet.extend(data0)
  # 添加第二個數(shù)據(jù)集的數(shù)據(jù)
  dataSet.extend(data1)
  # 對總的數(shù)據(jù)集進(jìn)行打亂（其實不打亂也沒事，只不過打亂一下直觀上讓人感覺已經(jīng)混合了
  # 讀者可以將下面這句話屏蔽以后看看效果是否有差別）
  random.shuffle(dataSet)

  #返回偽造好的數(shù)據(jù)集
  return dataSet

def calcGauss(dataSetArr, mu, sigmod):
  '''
  根據(jù)高斯密度函數(shù)計算值
  依據(jù)：“9.3.1 高斯混合模型” 式9.25
  注：在公式中y是一個實數(shù)，但是在EM算法中(見算法9.2的E步)，需要對每個j
  都求一次yjk，在本實例中有1000個可觀測數(shù)據(jù)，因此需要計算1000次?？紤]到
  在E步時進(jìn)行1000次高斯計算，程序上比較不簡潔，因此這里的y是向量，在numpy
  的exp中如果exp內(nèi)部值為向量，則對向量中每個值進(jìn)行exp，輸出仍是向量的形式。
  所以使用向量的形式1次計算即可將所有計算結(jié)果得出，程序上較為簡潔
  :param dataSetArr: 可觀測數(shù)據(jù)集
  :param mu: 均值
  :param sigmod: 方差
  :return: 整個可觀測數(shù)據(jù)集的高斯分布密度（向量形式）
  '''
  # 計算過程就是依據(jù)式9.25寫的，沒有別的花樣
  result = (1 / (math.sqrt(2*math.pi)*sigmod**2)) * np.exp(-1 * (dataSetArr-mu) * (dataSetArr-mu) / (2*sigmod**2))
  # 返回結(jié)果
  return result


def E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1):
  '''
  EM算法中的E步
  依據(jù)當(dāng)前模型參數(shù)，計算分模型k對觀數(shù)據(jù)y的響應(yīng)度
  :param dataSetArr: 可觀測數(shù)據(jù)y
  :param alpha0: 高斯模型0的系數(shù)
  :param mu0: 高斯模型0的均值
  :param sigmod0: 高斯模型0的方差
  :param alpha1: 高斯模型1的系數(shù)
  :param mu1: 高斯模型1的均值
  :param sigmod1: 高斯模型1的方差
  :return: 兩個模型各自的響應(yīng)度
  '''
  # 計算y0的響應(yīng)度
  # 先計算模型0的響應(yīng)度的分子
  gamma0 = alpha0 * calcGauss(dataSetArr, mu0, sigmod0)
  # 模型1響應(yīng)度的分子
  gamma1 = alpha1 * calcGauss(dataSetArr, mu1, sigmod1)

  # 兩者相加為E步中的分布
  sum = gamma0 + gamma1
  # 各自相除，得到兩個模型的響應(yīng)度
  gamma0 = gamma0 / sum
  gamma1 = gamma1 / sum

  # 返回兩個模型響應(yīng)度
  return gamma0, gamma1

def M_step(muo, mu1, gamma0, gamma1, dataSetArr):
  # 依據(jù)算法9.2計算各個值
  # 這里沒什么花樣，對照書本公式看看這里就好了
  mu0_new = np.dot(gamma0, dataSetArr) / np.sum(gamma0)
  mu1_new = np.dot(gamma1, dataSetArr) / np.sum(gamma1)

  sigmod0_new = math.sqrt(np.dot(gamma0, (dataSetArr - muo)**2) / np.sum(gamma0))
  sigmod1_new = math.sqrt(np.dot(gamma1, (dataSetArr - mu1)**2) / np.sum(gamma1))

  alpha0_new = np.sum(gamma0) / len(gamma0)
  alpha1_new = np.sum(gamma1) / len(gamma1)

  # 將更新的值返回
  return mu0_new, mu1_new, sigmod0_new, sigmod1_new, alpha0_new, alpha1_new


def EM_Train(dataSetList, iter=500):
  '''
  根據(jù)EM算法進(jìn)行參數(shù)估計
  算法依據(jù)“9.3.2 高斯混合模型參數(shù)估計的EM算法” 算法9.2
  :param dataSetList:數(shù)據(jù)集（可觀測數(shù)據(jù)）
  :param iter: 迭代次數(shù)
  :return: 估計的參數(shù)
  '''
  # 將可觀測數(shù)據(jù)y轉(zhuǎn)換為數(shù)組形式，主要是為了方便后續(xù)運算
  dataSetArr = np.array(dataSetList)

  # 步驟1：對參數(shù)取初值，開始迭代
  alpha0 = 0.5
  mu0 = 0
  sigmod0 = 1
  alpha1 = 0.5
  mu1 = 1
  sigmod1 = 1

  # 開始迭代
  step = 0
  while (step < iter):
    # 每次進(jìn)入一次迭代后迭代次數(shù)加1
    step += 1
    # 步驟2：E步：依據(jù)當(dāng)前模型參數(shù)，計算分模型k對觀測數(shù)據(jù)y的響應(yīng)度
    gamma0, gamma1 = E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1)
    # 步驟3：M步
    mu0, mu1, sigmod0, sigmod1, alpha0, alpha1 = M_step(mu0, mu1, gamma0, gamma1, dataSetArr)

  # 迭代結(jié)束后將更新后的各參數(shù)返回
  return alpha0, mu0, sigmod0, alpha1, mu1, sigmod1


if __name__ == '__main__':
  start = time.time()

  # 設(shè)置兩個高斯模型進(jìn)行混合，這里是初始化兩個模型各自的參數(shù)
  # 見“9.3 EM算法在高斯混合模型學(xué)習(xí)中的應(yīng)用”
  # alpha是“9.3.1 高斯混合模型” 定義9.2中的系數(shù)α
  # mu0是均值μ
  # sigmod是方差σ
  # 在設(shè)置上兩個alpha的和必須為1，其他沒有什么具體要求，符合高斯定義就可以
  alpha0 = 0.3 # 系數(shù)α
  mu0 = -2 # 均值μ
  sigmod0 = 0.5 # 方差σ

  alpha1 = 0.7 # 系數(shù)α
  mu1 = 0.5 # 均值μ
  sigmod1 = 1 # 方差σ

  # 初始化數(shù)據(jù)集
  dataSetList = loadData(mu0, sigmod0, mu1, sigmod1, alpha0, alpha1)

  #打印設(shè)置的參數(shù)
  print('---------------------------')
  print('the Parameters set is:')
  print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (
    alpha0, alpha1, mu0, mu1, sigmod0, sigmod1
  ))

  # 開始EM算法，進(jìn)行參數(shù)估計
  alpha0, mu0, sigmod0, alpha1, mu1, sigmod1 = EM_Train(dataSetList)

  # 打印參數(shù)預(yù)測結(jié)果
  print('----------------------------')
  print('the Parameters predict is:')
  print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (
    alpha0, alpha1, mu0, mu1, sigmod0, sigmod1
  ))

  # 打印時間
  print('----------------------------')
  print('time span:', time.time() - start)