EM算法与聚类应用

反思总结：
  （1）EM算法分为通俗来说分为E步和M步两部分，E步是通过当前假定的模型参数计算求得当前数据属于每个模型的概率，M步是通过E步求得的数据属于每个模型的概率更新之前假定的模型参数，这样一直循环迭代下去，直到结果趋于收敛为止。
  （2）采用迭代次数截止和数据之差小于特定值截止两次得到参数基本没有差别，但是当改变假定初值时，最终得到的参数变化较大，EM算法和初值选择有较大关系。
  （3）KMeans聚类是通过EM算法进行聚类的一种特殊情况，EM算法和GMM聚类除了能够区分空间距离上的不同，更重要的是能够对不同分布参数的数据进行区分，因此比KMeans算法具有更加强的聚类功能.

import numpy as np
import random
import math
import time

1.利用不同分布参数生成数据

def loadData(mu0, sigma0, mu1, sigma1, alpha0, alpha1):#利用不同分布参数生成数据
    length = 1000
    data0 = np.random.normal(mu0, sigma0, int(length * alpha0))
    data1 = np.random.normal(mu1, sigma1, int(length * alpha1))
    dataSet = []
    dataSet.extend(data0)
    dataSet.extend(data1)
    random.shuffle(dataSet)
    return dataSet

2.返回数据的高斯分布概率结果

def calcGauss(dataSetArr, mu, sigmod):#返回数据的高斯分布概率结果
    result = (1 / (math.sqrt(2*math.pi)*sigmod**2)) * np.exp(-1 * (dataSetArr-mu) * (dataSetArr-mu) / (2*sigmod**2))
    return result

3.E步计算数据属于每个模型的概率

def E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1):#E步计算数据属于每个模型的概率
    gamma0 = alpha0 * calcGauss(dataSetArr, mu0, sigmod0)
    gamma1 = alpha1 * calcGauss(dataSetArr, mu1, sigmod1)
    sum = gamma0 + gamma1
    gamma0 = gamma0 / sum
    gamma1 = gamma1 / sum
    return gamma0, gamma1

4.M步利用E步的结果更新每个模型的参数

def M_step(muo, mu1, gamma0, gamma1, dataSetArr):#M步利用E步的结果更新每个模型的参数
    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

5.重复E-M步，训练数据直到达到迭代次数为止

def EM_Train0(dataSetList, iter=500):#重复E-M步，训练数据直到达到迭代次数为止
    dataSetArr = np.array(dataSetList)
    alpha0 = 0.5
    mu0 = 0
    sigmod0 = 1
    alpha1 = 0.5
    mu1 = 1
    sigmod1 = 1
    step = 0
    while (step < iter):
        step += 1
        gamma0, gamma1 = E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1)
        mu0, mu1, sigmod0, sigmod1, alpha0, alpha1 = M_step(mu0, mu1, gamma0, gamma1, dataSetArr)
    return alpha0, mu0, sigmod0, alpha1, mu1, sigmod1

6.重复E-M步，训练数据直到两次更新参数之差的总和小于0.01为止

def EM_Train1(dataSetList):#重复E-M步，训练数据直到两次更新参数之差的总和小于0.01为止
    dataSetArr = np.array(dataSetList)
    alpha0 = 0.5
    mu0 = 0
    sigmod0 = 1
    alpha1 = 0.5
    mu1 = 1
    sigmod1 = 1
    res=[np.array([alpha0,mu0,sigmod0,alpha1,mu1,sigmod1])]
    step = 0
    while True:
        step += 1
        gamma0, gamma1 = E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1)
        mu0, mu1, sigmod0, sigmod1, alpha0, alpha1 = M_step(mu0, mu1, gamma0, gamma1, dataSetArr)
        res.append(np.array([alpha0,mu0,sigmod0,alpha1,mu1,sigmod1]))
        if np.abs(res[-1]-res[-2]).sum()<0.01:
            break
    return alpha0, mu0, sigmod0, alpha1, mu1, sigmod1

7.更改初始参数赋值后，训练数据返回最终参数

def EM_Train2(dataSetList):#更改初始参数赋值后，训练数据返回最终参数
    dataSetArr = np.array(dataSetList)
    alpha0 = 0.2
    mu0 = 0.5
    sigmod0 = 2
    alpha1 = 0.8
    mu1 = 0.6
    sigmod1 = 2
    res=[np.array([alpha0,mu0,sigmod0,alpha1,mu1,sigmod1])]
    step = 0
    while True:
        step += 1
        gamma0, gamma1 = E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1)
        mu0, mu1, sigmod0, sigmod1, alpha0, alpha1 = M_step(mu0, mu1, gamma0, gamma1, dataSetArr)
        res.append(np.array([alpha0,mu0,sigmod0,alpha1,mu1,sigmod1]))
        if np.abs(res[-1]-res[-2]).sum()<0.01:
            break
    return alpha0, mu0, sigmod0, alpha1, mu1, sigmod1

8.执行函数，生成数据，计算参数并打印

if __name__ == '__main__':
    start = time.time() 
    alpha0 = 0.3 
    mu0 = -2  
    sigmod0 = 0.5 

    alpha1 = 0.7  
    mu1 = 0.5  
    sigmod1 = 1  

    # 初始化数据集
    dataSetList = loadData(mu0, sigmod0, mu1, sigmod1, alpha0, alpha1)

    #打印设置的参数
    print('---------------------------')
    print('the Parameters set is:')
    print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f'