python实现EM聚类算法

python实现EM聚类算法

实验结果

当前的分类概率参数：PI_1：0.33334922292091，PI_2：0.66665077707909
当前的分类值： [0 0 1 0 1 1 1 1 1]
分类1 ['T100' 'T200' 'T400']
分类2 ['T300' 'T500' 'T600' 'T700' 'T800' 'T900']

基于下面的数据，采用EM算法将下面的数据聚成2类

E步公式：

M步公式：

​

1 数据工程：

创建dataset.scv文件，文件内容如下。$TID$为用户购买id，$I_i$表示商品类别，数字1表示用户购买了该商品

2 参数准备

• 该实验将类别分为2个分类，故设定${\pi }_1,{\pi }_2$为观察值来自两个类别中的概率。初始参数为${\pi }_1=0.6,{\pi }_2=0.4$
  • 例如：T100观察值来自类别${\pi }_1$的概率为0.6
• $x(i{)}_j$ 表示用户$i$购买商品$j$ 的情况，取值为0或者1
• ${\theta }_{kj}$ 表示在第$k$ 个分类中观察到$j=1$ 的概率
• $p(k|i)$ 表示第 i 条数据属于第 k 个分量的概率

3 计算

由EM算法步骤可得，EM分为两部分，即E和M部分

• E:对每个对象，计算它们属于各个分布的概率，即求出每一行对象，在分类1和2上面的概率

  • 公式为：$P(C_i|x_j)=\frac{f(x_j|{\mu }_i,{\sigma }_i^2)\bullet P(C_i)}{{\sum }_{a=1}^{k}f(x_j|{\mu }_a,{\sigma }_a^2)\bullet P(C_a)}$
  • 计算公式为：[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-M4s6XbGG-1653982609027)(机器学习作业.assets/20220531091848.png)]
  • 拆解计算
    • $P{C}_1$为当前用户使用分类1的概率乘以 参数分布情况概率的乘积
    • $P{C}_2$为当前用户使用分类2的概率乘以 参数分布情况概率的乘积
    • $P{C}_1=P{C}_1/(P{C}_1+P{C}_2)$ $P{C}_2=1-P{C}_1$ 表示对象在2类分类中的概率
  • 得到所有对象的分类概率$PC$

• M: 根据上步得出的所有对象所属分类的概率，计算能够最大化似然 函数的新的参数估计,直到参数不再改变

  • 计算新的$\pi$
    • 公式：${\pi }_k=\frac{1}{n}{\sum }_{i=1}^{n}P(k|x(i))$
    • 代码$ new_\pi[k]=mean(PC[:,k])$ 求当前每个分类下的均值，作为新的$\pi$
  • 更新 参数取值概率
    • 公式：${\theta }_{kj}^{new}=\frac{{\sum }_{i=1}^{n}p(k|i)x_j(i)}{{\sum }_{i=1}^{n}p(k|i)}$
    • 代码$ \theta[k]=((np.matrix(observations).T*np.matrix(PC)).T.getA())[k]/sum(PC[:,k])$

4 源代码：

# -*- coding: utf-8 -*- 
# @Time : 2022/5/30 21:06
# @Author : PVINCE
# @File : NewEm.py
#@desc:使用EM聚类


import math
import numpy as np
import pandas as pd

def E(observations,pi,nums_k,theta):
    ob_length=len(observations)
    PC = np.zeros((ob_length, nums_k), float)
    i = 0
    for observation in observations:
        j=0
        PC_1=pi[0]
        PC_2=pi[1]
        for colums in observation:
            PC_1*=math.pow(theta[0][j],colums) *math.pow((1-theta[0][j]),1-colums)
            PC_2*=math.pow(theta[1][j],colums) *math.pow((1-theta[1][j]),1-colums)
            j=j+1
        PC[i][0]=PC_1/(PC_1+PC_2)#得到I对象属于类别1的概率
        PC[i][1]=1-PC[i][0]##得到I对象属于类别2的概率
        i=i+1
    return PC#返回当前Pi下所有对象的预测值类别概率

def M(observations,PC,theta):
    # 更新Pi和theta
    new_pi=np.zeros(2,float)
    new_pi[0]=np.mean(PC[:,0])
    new_pi[1]=np.mean(PC[:,1])
    theta[0]=((np.matrix(observations).T*np.matrix(PC)).T.getA())[0]/sum(PC[:,0])
    theta[1]=((np.matrix(observations).T*np.matrix(PC)).T.getA())[1]/sum(PC[:,1])
    return new_pi,theta


# ====以下为预设值========
observations=[[1,1,0,0,1],
              [0,1,0,1,0],
              [0,1,1,0,0],
              [1,1,0,1,0],
              [1,0,1,0,0]]
lable=['T100','T200',"T300","T400","T500"]
# ====预设值ENd========  分类结果：分类1 ['T100' 'T200' 'T400']分类2 ['T300' 'T500']

# ========以下为CSV文件中的数据========
data_ori = pd.read_csv('./dataset.csv', encoding='utf-8')
features_scv = data_ori.columns[1:]
observations = data_ori[features_scv].values
lable=data_ori[data_ori.columns[0]].values
#=========csv结束======   分类结果：分类1 ['T100' 'T200' 'T400']，分类2 ['T300' 'T500' 'T600' 'T700' 'T800' 'T900']

#分类个数
nums_k=2
# 类别概率
pi=[0.6,0.4]
#类别属性选择概率：类别为k，故参数k行
theta=[[0.6,0.6,0.3,0.5,0.3],
        [0.4,0.4,0.7,0.5,0.7]
              ]
new_pi=np.array(pi)
# 比较两个分类参数的精度设置
dot_count=4
# 循环次数
count=0
while True:
    print("{}当前的分类概率参数：PI_1：{}，PI_2：{}".format(count,new_pi[0], new_pi[1]))
    print("{}当前的特征参数概率：".format(count), np.array(theta))
    count=count+1
    last_pi=new_pi
    PC = E(observations, new_pi, nums_k, theta)
    new_pi, theta = M(observations, PC, theta)
    if (new_pi <= 0).any():break
    if(np.round(new_pi,dot_count)==np.round(last_pi,dot_count)).all():break

print("当前的分类概率参数：PI_1：{}，PI_2：{}".format(last_pi[0], last_pi[1]))
print("当前的分类值：",np.argmax(np.array(PC),axis=1))
catat1=[]
catat2=[]
i=0
for idx in np.argmax(np.array(PC),axis=1):
    if idx==0:catat1.append(lable[i])
    if idx==1:catat2.append(lable[i])
    i=i+1
print("分类1",np.array(catat1))
print("分类2",np.array(catat2))

5 实验结果

inputs:k=2,pi=[0.6,0.4]  ,
		theta=[[0.6,0.6,0.3,0.5,0.3],
       		 [0.4,0.4,0.7,0.5,0.7]]
outputs:
当前的分类概率参数：PI_1：0.33334922292091，PI_2：0.66665077707909
当前的分类值： [0 0 1 0 1 1 1 1 1]
分类1 ['T100' 'T200' 'T400']
分类2 ['T300' 'T500' 'T600' 'T700' 'T800' 'T900']

参数变化

类别1的概率变化

类别2的概率变化