EM算法学习

EM算法学习记录

算法思路

实现代码

import math
import numpy

# 函数原型:
# Jensen 不等式
# f_2(x) >= 0 (凹函数);E(F(X)) >= F(E(x))

class gaussian_distribution:
def __init__(self,mu,sigma):
# 高斯模型
gaussian_distribution.mu = mu
# sigma代表标准差（math.sqrt(方差)）
gaussian_distribution.sigma = sigma
# 可能性
gaussian_distribution.probability = 0.0

def get_value(self,x):
# 高斯模型
if type(x) == type(list()):
value = numpy.exp(-((numpy.mat(x) - numpy.mat(self.mu))*(numpy.mat(x) - numpy.mat(self.mu)).transpose()).sum()/(2*self.sigma**2))
return value/(math.sqrt(2*numpy.pi)*self.sigma)
# 这里可以增加对于矩阵的判断
else :
# 二维
return numpy.exp(-(x - self.mu)**2/2*self.sigma**2)/(math.sqrt(2*numpy.pi)*self.sigma)

def change_mu(self,mu):
self.mu = mu

def change_sigma(self,sigma):
self.sigma = sigma

# receive a list of the value
# return a list include average value
def get_average_value(self,value_list):
sum = numpy.zeros((1,len(value_list[0])))
for list_reader in value_list:
sum += numpy.mat(list_reader)
sum /= len(value_list)
return sum.tolist()[0]

def get_sigma(self,value_list):
sum = 0.0
for list_reader in value_list:
sum += numpy.mat(list_reader)*numpy.mat(list_reader).transpose()
sum /= len(value_list)
return math.sqrt(sum)

class EM_gauss:
def __init__(self,k,db):
EM_gauss.gauss_dis_list = list()
EM_gauss.data_list = db
# 初始化各个模型
# k代表初始有k个高斯模型
for reader in range(k):
gauss_dis = gaussian_distribution(0,1)
gauss_dis.mu = gauss_dis.get_average_value(db)
gauss_dis.sigma = gauss_dis.get_sigma(db)
gauss_dis.probability = 1/k
EM_gauss.gauss_dis_list.append(gauss_dis)

def get_data_probability(self,data_index,model_index):
value = self.gauss_dis_list[model_index].probability*self.gauss_dis_list[model_index].get_value(self.data_list[data_index])
sum = 0.0
for gauss_dis in self.gauss_dis_list:
sum += gauss_dis.probability*gauss_dis.get_value(self.data_list[data_index])
return value/sum

# 模型存在的可能性
def get_model_probability(self,model_index):
sum = 0.0
for i in range(len(self.data_list)):
sum += self.get_data_probability(i-1,model_index)
return sum

def round(self):
# E
model_probability_list = list()
# 每个模型获取它存在的可能性
for i in range(len(self.gauss_dis_list)):
model_probability_list.append(self.get_model_probability(i-1))
# M
# 维护μ
for i in range(len(self.gauss_dis_list)):
mu = numpy.zeros((1,len(self.gauss_dis_list[0].mu)))
for j in range(len(self.data_list)):
mu += self.get_data_probability(j-1,i-1)*numpy.mat(self.data_list[j-1])
self.gauss_dis_list[i-1].mu = (mu/model_probability_list[i-1]).tolist()[0]
# 维护σ
for i in range(len(self.gauss_dis_list)):
sigma = 0.0
for j in range(len(self.data_list)):
sigma += self.get_data_probability(j-1,i-1)*(numpy.mat(self.data_list[j-1])-self.gauss_dis_list[i-1].mu)*(numpy.mat(self.data_list[j-1])-self.gauss_dis_list[i-1].mu).transpose()
self.gauss_dis_list[i-1].sigma = (sigma/model_probability_list[i-1]).tolist()[0][0]
# 维护每个数据符合高斯模型的可能性
for i in range(len(self.gauss_dis_list)):
self.gauss_dis_list[i-1].probability = model_probability_list[i-1]/len(self.data_list)

def rounds(self,error = 0.0001):
# 多轮运算，直到获得误差小于1-error的模型
record_last = list()
null_list = list()
for list_reader in self.gauss_dis_list:
if list_reader.sigma == 0:
null_list.append(elf.gauss_dis_list.index(list_reader))
continue
record_last += list_reader.mu
record_last.append(list_reader.sigma)
while null_list != []:
index = null_list.pop()
del self.gauss_dis_list[index]
self.round()
while not self.check_gauss_equal(record_last):
# for i in self.gauss_dis_list:
#     print(i.mu)
#     print(i.sigma)
#     print(i.probability)
record_last = []
for list_reader in self.gauss_dis_list:
if list_reader.sigma == 0:
null_list.append(self.gauss_dis_list.index(list_reader))
continue
record_last += list_reader.mu
record_last.append(list_reader.sigma)
while null_list != []:
index = null_list.pop()
del self.gauss_dis_list[index]
self.round()
# for i in self.gauss_dis_list:
#     print(i.mu)
#     print(i.sigma)
#     print(i.probability)

# 需保证两个list长度相同
# 检查两个高斯模型是否等价
def check_gauss_equal(self,list_a,error = 0.0001):
list_b = list()
for list_reader in self.gauss_dis_list:
list_b += list_reader.mu
list_b.append(list_reader.sigma)
for i in range(len(list_a)):
if abs(list_a[i-1] - list_b[i-1]) > error:
return False
return True

# 数据的标准格式
# num_x num_y num_z ...
# ...
# (均为数字)
class point_data_reader:
# 获取点类数据
file_name = str()
def __init__(self,file_name):
point_data_reader.file_name = file_name

def get_data_list(self,num_lost):
file_ = open(self.file_name,'r+')
db = list()
for line in file_.readlines():
reader_list = list()
for reader_ in line.split():
try:
eval(reader_)
except:
reader_list.append(num_lost)
else:
reader_list.append(eval(reader_))
db.append(reader_list)
file_.close()
return db

if __name__ == '__main__':
# num_lost为填充丢失数据
db = point_data_reader('text.dat').get_data_list(num_lost = 1.0)
# 2为模型个数
em = EM_gauss(2,db)
em.rounds()