转载请标明作者
本文主要使用了EM算法计算隐变量,计算使用的是假数据,使用了三个高斯进行混合,由于时间,代码只是进行了部分优化,所以运算有点慢,还有图的显示以后有时间再改吧,代码如下:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy import stats
np.random.seed(1234)
x_plot_1,x_plot_2=np.mgrid[-10:10:50j,-10:10:50j]
x_plot=np.stack((x_plot_1,x_plot_2),axis=2)
u=[[-5,-3],[3,-1],[-3,5]]#三个二元高斯的均值
#真实数据的均值,方差和alpha值
sigma=[np.identity(2),np.diag([2,3]),np.array([[2,1],[1,2]])]#三个高斯的方差
#data1=np.random.multivariate_normal()
alpha=[0.1,0.5,0.4]
data_num=600
#上面是混合高斯的所有参数
#下面开始制造假的的数据
def fake_data_gen(u,sigma,alpha,number):
norm_list=[]
data_list=[]
for ga_par in zip(u,sigma,alpha):
norm=stats.multivariate_normal(ga_par[0],ga_par[1])
norm_list.append(norm)
for i in range(number):
rand=np.random.uniform(0,1)
if i%200==0:
print("have generrate {} samples".format(i))
if rand<=alpha[0]:
data=norm_list[0].rvs(1)
elif alpha[0]<rand<=(1-alpha[2]):
data=norm_list[1].rvs(1)
else:
data=norm_list[2].rvs(1)
data_list.append(data)
return np.array(data_list)
data=fake_data_gen(u,sigma,alpha,data_num)
#画出假数据的散点图
def plot_mygraph(u,sigma):
plt.cla()
plt.scatter(data[:, 0], data[:, 1])
for i in range(3):
pdf_x = stats.multivariate_normal(mean=u[i], cov=sigma[i]).pdf(x_plot)
plt.contour(x_plot_1,x_plot_2,pdf_x)
#初始化所有的参数,u,sigma,alpha,为以后的迭代做准备
u=np.array([[1,1],[3,3],[6,6]])
sigma=np.array([np.identity(2),np.identity(2),np.identity(2)])
alpha=np.array([0.5,0.1,0.4])
#开始迭代
#是哪个分布的后验概率
def poserier_vec_sum(u,sigam,data,l):
pdf_data_l=stats.multivariate_normal(u[l],sigma[l]).pdf(data)
sum_data=0
for i in range(3):
pdf_data_i=stats.multivariate_normal(u[i],sigma[i]).pdf(data)
sum_data+=pdf_data_i
sum_data=pdf_data_l/sum_data
return sum(sum_data)
def posterier(u,sigam,x_i,l):
pdf_x_i=stats.multivariate_normal(u[l],sigma[l]).pdf(x_i)
p_x_sum=sum([stats.multivariate_normal(u[i],sigma[i]).pdf(x_i) for i in range(3)])
return pdf_x_i/p_x_sum
def com_alpha(u,sigam,l,data,num):
sum1=poserier_vec_sum(u,sigma,data,l)
return sum1/(num)
def com_u(u,sigma,l,data):
sum1 = poserier_vec_sum(u, sigma, data, l)
sum2=0
for x_i in (data):
sum2+=x_i*posterier(u, sigma, x_i, l)
return sum2/sum1
def com_sigma(new_u,u,sigma,l,data):
sum2=0
sum1 = poserier_vec_sum(u, sigma, data, l)
for x_i in data:
x_i_u=np.reshape(x_i,(1,2))-np.reshape(new_u[l],(1,2))
x_i_u_T=np.reshape(x_i,(2,1))-np.reshape(new_u[l],(2,1))
sum2+=np.matmul(x_i_u_T,x_i_u)*posterier(u,sigma,x_i,l)
return sum2/sum1
def gussian_mixture_iter(u,sigma,alpha,data,iter_num):
for i in range(iter_num):
plot_mygraph(u, sigma)
for l in range(3):
alpha[l]=com_alpha(u,sigma,l,data,data_num)
alpha = alpha / sum(alpha)
new_u=[]
for l in range(3):
u_tmp=com_u(u,sigma,l,data)
new_u.append(u_tmp)
sigma_list=[]
for l in range(3):
sigma_tmp=com_sigma(new_u,u,sigma,l,data)
sigma_list.append(sigma_tmp)
u=new_u
sigma=sigma_list
alpha = alpha / sum(alpha)
if i%1==0:
print("第{}次迭代均值为:{}".format(i,u))
print("第{}次迭代方差为:{}".format(i,sigma))
print("第{}次迭代alpha为:{}".format(i, alpha))
plt.show()
return u,sigma,alpha
u,sigma,alpha=gussian_mixture_iter(u,sigma,alpha,data,100)
,第一次初始化后的图像为:
第二次地迭代:
第三次迭代:
第四次迭代:
第6次迭代: