二维高斯模型

二维高斯混合模型GMM图形化简单明了
2014-11-25 20:48:35
在matlab上实现的二维高斯混合模型是基于界面形式的直接画出高斯分布更方便的了解GMM 你值得拥有机器学习
收起

EM算法在二维高斯混合模型参数估计中的应用

2019-04-26 09:53:46

高斯混合模型高斯混合模型是指具有如下形式的概率分布模型：其中，是系数，，；是高斯分布密度，， ...称为第k个分模型。...对于二维高斯混合模型，d=2，y和都是二维的数据，用矩阵表示就是一行...


	高斯混合模型
	
高斯混合模型是指具有如下形式的概率分布模型：



其中，是系数，，；是高斯分布密度，，



称为第k个分模型。

参考：《统计学习方法》9.3 EM算法在高斯混合模型学习中的应用


	多维高斯混合模型
	
多维高斯混合模型具有如下形式的概率分布模型：



其中d为数据的维度，为均值，为协方差矩阵。

对于二维高斯混合模型，d=2，y和都是二维的数据，用矩阵表示就是一行两列，则是两行两列的协方差矩阵。


	二维高斯混合模型参数估计的EM算法
	
输入：观测数据，二维高斯混合模型；

输出：二维高斯混合模型参数。

算法步骤如下：

（1）取参数的初始值开始迭代

（2）E步：依据当前模型参数，计算分模型k对观测数据的响应度：



（3）M步：计算新一轮迭代的模型参数：







（4）重复第(2)步和第(3)步，直到收敛。


	实验步骤
	
1.确定二维高斯混合模型的参数

我们的二维高斯混合模型为两个权重不同的二维单高斯分布组合而成，所以可以简化为一个参数，即第一个二维单高斯分布的权重为π，第二个二维单高斯分布的权重为1-π。选择π值为0.8，第一个二维单高斯分布的均值为，协方差矩阵为，第二个二维单高斯分布的均值为，协方差矩阵为。

import numpy as np
import matplotlib.pyplot as plt


N = 1000
pi = 0.8
mean_1 = np.array([0.0, 0.0])
cov_1 = np.mat([[1.0, 0.0], [0.0, 1.0]])
mean_2 = np.array([3.0, 3.0])
cov_2 = np.mat([[1.0, 0.5], [0.5, 1.0]])


class GMM:
    def __init__(self, N, pi, mean_1, cov_1, mean_2, cov_2):
        # 采样点数
        self.N = N
        # 高斯分布的权重
        self.pi = pi
        # 第一个二维高斯分布的均值
        self.mean_1 = mean_1
        # 第一个二维高斯分布的协方差矩阵
        self.cov_1 = cov_1
        # 第二个二维高斯分布的均值
        self.mean_2 = mean_2
        # 第二个二维高斯分布的协方差矩阵
        self.cov_2 = cov_2

 该二维高斯混合模型的大致图像如下：



 2.采样

选取采样点为1000个。采样方法为在1000次循环中，每次产生一个0-1的随机数，如果小于0.8，则以第一个二维单高斯分布的概率密度函数生成一个样本点；如果大于0.8，则以第二个二维单高斯分布的概率密度函数生成一个样本点，最后以生成的1000个样本点作为观测数据集。

# 生成观测数据集
def dataset(self):
    # 数据集大小为(N,2)
    D = np.zeros(shape=(self.N, 2))
    for i in range(self.N):
        # 产生0-1之间的随机数
        j = np.random.random()
        if j < self.pi:
            # pi的概率以第一个二维高斯分布产生一个样本点
            x = np.random.multivariate_normal(mean=self.mean_1, cov=self.cov_1, size=1)
        else:
            # 1-pi的概率以第二个二维高斯分布产生一个样本点
            x = np.random.multivariate_normal(mean=self.mean_2, cov=self.cov_2, size=1)
        D[i] = x
    return D

 生成的观测数据集的分布图如下：

gmm = GMM(N, pi, mean_1, cov_1, mean_2, cov_2)
plt.figure()
plt.gca().set_aspect('equal')  # 令x轴和y轴的同一区间的刻度相同
D = gmm.dataset()
plt.scatter(D[:, 0], D[:, 1])
plt.show()



3.利用EM算法进行参数估计

（1）参数初始化

选取参数的初始值为π为0.5，第一个均值为，第一个协方差矩阵为，第二个均值为，第二个协方差矩阵为。

def EM(self, D, N):
    p = 0.5
    m_1 = np.array([[0.0, 0.0]])
    c_1 = np.mat([[1.0, 0.0], [0.0, 1.0]])
    m_2 = np.array([[1.0, 1.0]])
    c_2 = np.mat([[1.0, 0.0], [0.0, 1.0]])

（2）E步

依据当前模型参数，计算分模型k对观测数据的响应度，由于该模型只有两个二维单高斯分布组成，故K=2：



其中的表达式为：



    # 高斯分布密度phi
    def phi(x, mean, cov):
        return np.exp(-(x - mean) * np.linalg.pinv(cov) * (x - mean).T / 2) /
                (2 * np.pi * np.sqrt(np.linalg.det(cov)))

    # 分模型对观测数据的响应度gamma
    def gamma(D, j, k, p, m_1, c_1, m_2, c_2):
        # 第一个分模型
        if k == 0:
            return p * phi(D[j], m_1, c_1) / (p * phi(D[j], m_1, c_1) + (1 - p) * phi(D[j], m_2, c_2))
        # 第二个分模型
        elif k == 1:
            return (1 - p) * phi(D[j], m_2, c_2) / (p * phi(D[j], m_1, c_1) + (1 - p) * phi(D[j], m_2, c_2))

（3）M步

计算新一轮迭代的模型参数：



    def mu(D, g, N):
        # 第一个分模型的均值的分子部分
        mu_1a = np.array([[0.0, 0.0]])
        # 第一个分模型的均值的分母部分
        mu_1b = np.array([[0.0, 0.0]])
        # 第二个分模型的均值的分子部分
        mu_2a = np.array([[0.0, 0.0]])
        # 第二个分模型的均值的分母部分
        mu_2b = np.array([[0.0, 0.0]])
        for j in range(N):
            mu_1a += g[j][0] * D[j]
            mu_1b += g[j][0]
            mu_2a += g[j][1] * D[j]
            mu_2b += g[j][1]
        # 返回第一个分模型的均值和第二个分模型的均值，都是一行两列的矩阵
        return mu_1a / mu_1b, mu_2a / mu_2b



    # 模型的协方差矩阵sigma
    def sigma(D, m1, m2, g, N):
        # 第一个分模型的协方差矩阵的分子部分
        sigma_1a = np.mat([[0.0, 0.0], [0.0, 0.0]])
        # 第一个分模型的协方差矩阵的分母部分
        sigma_1b = np.mat([[0.0, 0.0], [0.0, 0.0]])
        # 第二个分模型的协方差矩阵的分子部分
        sigma_2a = np.mat([[0.0, 0.0], [0.0, 0.0]])
        # 第二个分模型的协方差矩阵的分母部分
        sigma_2b = np.mat([[0.0, 0.0], [0.0, 0.0]])
        for j in range(N):
            sigma_1a += g[j][0] * ((D[j] - m1).T * (D[j] - m1))
            sigma_1b += g[j][0]
            sigma_2a += g[j][1] * ((D[j] - m2).T * (D[j] - m2))
            sigma_2b += g[j][1]
        # 返回第一个分模型的协方差矩阵和第二个分模型的协方差矩阵，都是两行两列的矩阵
        return sigma_1a / sigma_1b, sigma_2a / sigma_2b



由于该模型只有两个二维单高斯分布组成，故只需要求第一个分模型的权重π，第二个分模型的权重1-π即可得到。

即只需要求下式：



    # 模型的权重alpha
    def alpha(g, N):
        a = 0
        for j in range(N):
            # 只需要求第一个分模型的权重
            a += g[j][0]
        return a / N

（4）重复第2步和第3步，直到收敛

标准做法是重复以上计算，直到对数似然函数值不再有明显变化为止。这里为了计算方便，当权重π的变化量小于时，停止迭代。

    # 权重pi的变化量初始化为1
    delta_p = 1
    # 迭代次数i
    i = 0
    # 当权重pi的变化量小于10e-7时停止迭代
    while delta_p > 10e-7:
        Gamma = np.zeros(shape=(N, 2))
        for j in range(N):
            for k in range(2):
                # 计算分模型k对观测数据的响应度gamma
                Gamma[j][k] = gamma(D, j, k, p, m_1, c_1, m_2, c_2)
            
        # 更新模型的均值
        m_1, m_2 = mu(D, Gamma, N)
        # 更新模型的协方差矩阵
        c_1, c_2 = sigma(D, m_1, m_2, Gamma, N)
        # 计算权重pi的变化量
        delta_p = abs(p - alpha(Gamma, N))
        # 更新模型的权重
        p = alpha(Gamma, N)

        i += 1
        # 每五次迭代打印权重值
        if i % 5 == 0:
            print(i, "steps' pi:", p)

返回EM算法得到的五个参数值、

    # 返回EM算法学习到的五个参数
    return p, m_1, c_1, m_2, c_2

4.将学习到的参数与设定值比较

pi_learn, mean_1_learn, cov_1_learn, mean_2_learn, cov_2_learn = gmm.EM(D, N)
print("权重值：", pi_learn)
print("第一个分模型的均值：\n", mean_1_learn)
print("第一个分模型的协方差矩阵：\n", cov_1_learn)
print("第二个分模型的均值：\n", mean_2_learn)
print("第二个分模型的协方差矩阵：\n", cov_2_learn)

50 steps' pi: 0.7930182035003274
权重值： 0.7930182035003274
第一个分模型的均值： 
[[-0.08556625 -0.03833228]]
第一个分模型的协方差矩阵： 
[[ 0.97774587 -0.04718663]
 [-0.04718663  0.88868208]]
第二个分模型的均值： 
[[2.89290133 3.01754719]]
第二个分模型的协方差矩阵： 
[[1.06842288 0.45923016]
 [0.45923016 0.93972666]]

经过50次迭代后停止学习，50 steps' pi: 0.7930182035003274

得到权重值π为0.7930182035003274

第一个二维单高斯分布的均值为[[-0.08556625 -0.03833228]]

协方差矩阵为

[[ 0.97774587 -0.04718663]

 [-0.04718663  0.88868208]]

第二个二维单高斯分布的均值为[[2.89290133 3.01754719]]

协方差矩阵为

[[1.06842288 0.45923016]

 [0.45923016 0.93972666]]

可以看到，EM算法学习到的二维高斯混合模型的参数接近于设定值。


	考虑采样点数N对参数估计的影响
	
比较N取100,1000,10000时参数估计的准确度，由于一次参数估计的随机性太大，故不同的N值各进行10次参数估计并计算权重值π的均值和方差。

# 计算steps次参数估计的权重值的均值和方差
def pi_N(N, pi, mean_1, cov_1, mean_2, cov_2, steps):
    gmm = GMM(N, pi, mean_1, cov_1, mean_2, cov_2)
    pi_steps = np.zeros(shape=steps)
    # 均值
    pi_mu = 0
    # 方差
    pi_sigma = 0

    # 计算steps次估计值和均值
    for i in range(steps):
        D = gmm.dataset()
        pi_learn, _, _, _, _ = gmm.EM(D, N)
        pi_mu += pi_learn
        pi_steps[i] = pi_learn

    pi_mu /= steps

    # 计算steps次方差
    for i in range(steps):
        pi_sigma += (pi_steps[i]-pi_mu)**2

    pi_sigma /= steps

    return pi_mu, pi_sigma


（1）取N=100



# 参数估计次数steps
steps = 10

Y_mu = []
Y_sigma = []

pi_mu_100, pi_sigma_100 = pi_N(N1, pi, mean_1, cov_1, mean_2, cov_2, steps)
Y_mu.append(pi_mu_100)
Y_sigma.append(pi_sigma_100)
print("N=100时学习", steps, "次得到的权重值均值：", pi_mu_100, ",方差：", pi_sigma_100)

N=100时学习 10 次得到的权重值均值： 0.7482191862720133 ,方差： 0.0024508318225135105

（2）取N=1000



pi_mu_1000, pi_sigma_1000 = pi_N(N2, pi, mean_1, cov_1, mean_2, cov_2, steps)
Y_mu.append(pi_mu_1000)
Y_sigma.append(pi_sigma_1000)
print("N=1000时学习", steps, "次得到的权重值均值：", pi_mu_1000, ",方差：", pi_sigma_1000)

N=1000时学习 10 次得到的权重值均值： 0.8033618327297114 ,方差： 0.0001821331641820217

（3）取N=10000



pi_mu_10000, pi_sigma_10000 = pi_N(N3, pi, mean_1, cov_1, mean_2, cov_2, steps)
Y_mu.append(pi_mu_10000)
Y_sigma.append(pi_sigma_10000)
print("N=10000时学习", steps, "次得到的权重值均值：", pi_mu_10000, ",方差：", pi_sigma_10000)


N=10000时学习 10 次得到的权重值均值： 0.7980447763980503 ,方差： 1.5680138676242733e-05

 

将不同的N值所学习到的权重值的均值与设定值(0.8)相比较，并画在柱状图中，可以看到N越大，其学习到的权重值的均值越接近于设定值0.8。

Y_mu.append(0.8)
X = [1, 2, 3, 4]
plt.bar(X, Y_mu, align='center', tick_label=['N=100', 'N=1000', 'N=10000', 'standard value'], width=0.5)
plt.ylim((0.77, 0.82))
plt.xlabel('different N')
plt.ylabel('mean of pi')
plt.show()



由于N=10000时方差接近为0，不方便画图，从数据上直观来看，N越大，其学习到的权重值的方差越小且接近于0.

总体上来说，样本点数N取的越大，使用EM算法学习高斯混合模型的参数效果越好。

收起

展开全文

EM算法--二维高斯混合模型(GMM)

万次阅读 2017-04-10 19:19:16

EM算法是一种迭代算法，1977年由Dempster等人总结出，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望(expectation);... 首先准备一些预备知识，如：二维

   参考文章
   http://blog.163.com/baolong_zhu/blog/static/196311091201421185531966/
   《统计学习方法》 李航
   EM算法是一种迭代算法，1977年由Dempster等人总结出，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望(expectation);M步，求极大(maximization).所以这一算法称为期望极大算法(expectation maximization algorithm),简称EM算法。
   首先准备一些预备知识，如：二维高斯混合函数的形式，协方差的定义
高斯混合函数
        对于一维混合高斯函数，如(1)式所示。(所谓的混合高斯，就是将一些高斯函数利用一些权值系数组合起来)  （1）
        二维高斯混合函数，如(2)式所示。
   （2）
       可以看出二维高斯混合模型，其中的d为数据的维度，在这里等于2；∑  代表协方差矩阵，两行两列；这里的x以及μ都是二维的数据，用矩阵表示就是一行两列。这里提到了协方差矩阵，那么什么是协方差矩阵，与方差有什么区别？
协方差矩阵
       理解协方差矩阵的关键在于牢记它计算的是不同维度之间的协方差，而不是不同样本之间。拿到一个样本矩阵，我们最先明确的是一行是一个样本还是一个维度（重要）。
        方差表示的是每个样本值与全体样本值的平均数之差的平方值的平均数，反映了一组数据的离散程度。
 
         既然有了方差还要协方差干啥？显然，方差只能衡量一维数据的离散程度，对于二维数据，就用到了协方差，如下式所示。其中X，Y代表一组数据的两个维度。COV(X,Y)的值大于零，表示两个维度的数据成正相关，小于零则表示成负相关，等于零表示两个维度的数据相互独立。显然当X=Y时，COV(X,X)计算的就是方差。由协方差的定义可以得到COV(X,Y)=COV(Y,X)。

   二维数据的协方差矩阵表示为 ，显然协方差矩阵为对角矩阵。
    三维数据的协方差矩阵表示为，同理知道N维数据的协方差矩阵。
EM算法
   假设一组数据x1,x2,...,xn是来自上述二维高斯混合模型，这里我们利用EM算法确定各个高斯部件的参数，充分拟合给定的数据。算法目的具体的理论知识与算法步骤参考《统计学习方法》李航一书中对于EM算法的讲解(李航的这一本书，绝对值得学习)。
   我主要写一下自己的理解，也为自己以后复习用。
   EM算法其实就是概率模型参数的极大似然估计法，只不过存在一些隐变量，从而不能像一般的概率模型参数估计一样直接求似然函数的最大值，因为这里的似然函数存在隐变量，因而不存在解析解，就不能直接对 似然函数求偏导来求极大值，只能迭代求解。
   基本思想：
   E步，假设参数已知，去估计隐藏变量的期望；
   M步，利用E步求得的隐藏变量的期望，根据最大似然估计求得参数，然后新得到的参数值重新被用于E步。。。直至收敛到局部最优解。
  需要注意的一点是：理论中，E步最重要的是求Q函数（李航 9.29式），即完全数据的对数似然函数的期望，M步要做的是用Q函数对各参数求偏导，从而确定满足Q函数最大的各参数的值。而在程序中，因为M步Q函数对参数求偏导后，各参数的值由Q函数中的Gamma(j,k)决定，所以E步只需计算Gamma(j,k)即可，Gamma(j,k)代表第j个样本属于第k个高斯模型的概率，这样M步各参数的值也就确定了，所以切记实例理论与程序实现的区别，好像是程序只是利用了理论的精简结果。
学习完理论之后，最好写个程序练习一下，加深自己的理解。下面的程序(matlab编写)，是这样一个过程：
1.首先产生两组数据点。每一组数据点都是符合给定均值Mu和协方差矩阵∑的。程序和图如下所示，
%产生2个二维正态数据
MU1 = [1 2];
SIGMA1 = [1 0; 0 0.5];
MU2 = [-1 -1];
SIGMA2 = [1 0; 0 1];
X = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];
scatter(X(:,1),X(:,2),10,'*');


  2.对于上述产生的数据，利用EM算法预测这组数据对应的均值μ和协方差矩阵∑
   对于各参数初值的设置：两个均值初始值各设置为两组数据点中的随机一个样本点，权值系数初始值设置为1/k,k是高斯混合模型中模型的个数，协方差矩阵初值设置为所有数据的协方差。
 程序：GMM.m
%程序中Psi对应二维高斯函数，其中的(x-μ)与其转置的顺序与上述提到的高斯函数不同，这样是为了保证矩阵可乘性，不影响结果.
%Gamma 为隐变量的值，Gamma(i,j)代表第i个样本属于第j个模型的概率。
%Mu为期望，Sigma为协方差矩阵%Pi为各模型的权值系数%产生2个二维正态数据
%产生2个二维正态数据
MU1 = [1 2];
SIGMA1 = [1 0; 0 0.5];
MU2 = [-1 -1];
SIGMA2 = [1 0; 0 1];
%生成1000行2列(默认)均值为mu标准差为sigma的正态分布随机数
X = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];
%显示
scatter(X(:,1),X(:,2),10,'.');
%====================
K=2;
[N,D]=size(X);
Gamma=zeros(N,K);
Psi=zeros(N,K);
Mu=zeros(K,D);
LM=zeros(K,D);
Sigma =zeros(D, D, K); 
Pi=zeros(1,K);
%选择随机的两个样本点作为期望迭代初值
Mu(1,:)=X(randint(1,1,[1 N/2]),:);
Mu(2,:)=X(randint(1,1,[N/2 N]),:);
%所有数据的协方差作为协方差初值
for k=1:K
  Pi(k)=1/K;
  Sigma(:, :, k)=cov(X);
end
LMu=Mu;
LSigma=Sigma;
LPi=Pi;
while true
%Estimation Step  
for k = 1:K
  Y = X - repmat(Mu(k,:),N,1);
  Psi(:,k) = (2*pi)^(-D/2)*det(Sigma(:,:,k))^(-1/2)*diag(exp(-1/2*Y/(Sigma(:,:,k))*(Y')));      %Psi第一列代表第一个高斯分布对所有数据的取值
end
Gamma_SUM=zeros(D,D);
for j = 1:N
  for k=1:K
  Gamma(j,k) = Pi(1,k)*Psi(j,k)/sum(Psi(j,:)*Pi');                                               %Psi的第一行分别代表两个高斯分布对第一个数据的取值
  end
end
%Maximization Step
for k = 1:K
%update Mu
  Mu_SUM= zeros(1,D);
  for j=1:N
     Mu_SUM=Mu_SUM+Gamma(j,k)*X(j,:);   
  end
  Mu(k,:)= Mu_SUM/sum(Gamma(:,k));
%update Sigma
 Sigma_SUM= zeros(D,D);
 for j = 1:N
   Sigma_SUM = Sigma_SUM+ Gamma(j,k)*(X(j,:)-Mu(k,:))'*(X(j,:)-Mu(k,:));
 end
 Sigma(:,:,k)= Sigma_SUM/sum(Gamma(:,k));
 %update Pi
 Pi_SUM=0;
 for j=1:N
   Pi_SUM=Pi_SUM+Gamma(j,k);
 end
 Pi(1,k)=Pi_SUM/N;
end

R_Mu=sum(sum(abs(LMu- Mu)));
R_Sigma=sum(sum(abs(LSigma- Sigma)));
R_Pi=sum(sum(abs(LPi- Pi)));
R=R_Mu+R_Sigma+R_Pi;
if ( R<1e-10)
    disp('期望');
    disp(Mu);
    disp('协方差矩阵');
    disp(Sigma);
    disp('权值系数');
    disp(Pi);
    obj=gmdistribution(Mu,Sigma);
    figure,h = ezmesh(@(x,y)pdf(obj,[x,y]),[-8 6], [-8 6]);
    break;
end
 LMu=Mu;
 LSigma=Sigma;
 LPi=Pi;
end
%=====================


%%%%%%%%运行结果%%%%%%%
%
% 期望
%    1.0334    2.0479
%    -0.9703   -0.9879
% 
% 协方差矩阵
% 
% (:,:,1) =
% 
%     0.9604   -0.0177
%    -0.0177    0.4463
% 
% 
% (:,:,2) =
% 
%     0.9976    0.0667
%     0.0667    1.0249
% 
% 权值系数
%     0.4935    0.5065
%%%%%%%%%%%%%%%%%%%%%%%%%

收起

展开全文

EM算法--二维高斯混合模型(GMMs)

2019-03-30 20:58:21

参考文章 ...EM算法是一种迭代算法，1977年由Dempster等人总结出，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望(expectation);M步...

参考文章

   http://blog.163.com/baolong_zhu/blog/static/196311091201421185531966/

   《统计学习方法》 李航

   EM算法是一种迭代算法，1977年由Dempster等人总结出，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望(expectation);M步，求极大(maximization).所以这一算法称为期望极大算法(expectation maximization algorithm),简称EM算法。

   首先准备一些预备知识，如：二维高斯混合函数的形式，协方差的定义

高斯混合函数

        对于一维混合高斯函数，如(1)式所示。(所谓的混合高斯，就是将一些高斯函数利用一些权值系数组合起来)  （1）



二维高斯混合函数，如(2)式所示。



   （2）

       可以看出二维高斯混合模型，其中的d为数据的维度，在这里等于2；∑  代表协方差矩阵，两行两列；这里的x以及μ都是二维的数据，用矩阵表示就是一行两列。这里提到了协方差矩阵，那么什么是协方差矩阵，与方差有什么区别？

协方差矩阵

       理解协方差矩阵的关键在于牢记它计算的是不同维度之间的协方差，而不是不同样本之间。拿到一个样本矩阵，我们最先明确的是一行是一个样本还是一个维度（重要）。

        方差表示的是每个样本值与全体样本值的平均数之差的平方值的平均数，反映了一组数据的离散程度。

 

         既然有了方差还要协方差干啥？显然，方差只能衡量一维数据的离散程度，对于二维数据，就用到了协方差，如下式所示。其中X，Y代表一组数据的两个维度。COV(X,Y)的值大于零，表示两个维度的数据成正相关，小于零则表示成负相关，等于零表示两个维度的数据相互独立。显然当X=Y时，COV(X,X)计算的就是方差。由协方差的定义可以得到COV(X,Y)=COV(Y,X)。

   二维数据的协方差矩阵表示为 ，显然协方差矩阵为对角矩阵。

    三维数据的协方差矩阵表示为，同理知道N维数据的协方差矩阵。

EM算法

   假设一组数据x1,x2,...,xn是来自上述二维高斯混合模型，这里我们利用EM算法确定各个高斯部件的参数，充分拟合给定的数据。算法目的具体的理论知识与算法步骤参考《统计学习方法》李航一书中对于EM算法的讲解(李航的这一本书，绝对值得学习)。

   我主要写一下自己的理解，也为自己以后复习用。

   EM算法其实就是概率模型参数的极大似然估计法，只不过存在一些隐变量，从而不能像一般的概率模型参数估计一样直接求似然函数的最大值，因为这里的似然函数存在隐变量，因而不存在解析解，就不能直接对 似然函数求偏导来求极大值，只能迭代求解。

   基本思想：

   E步，假设参数已知，去估计隐藏变量的期望；

   M步，利用E步求得的隐藏变量的期望，根据最大似然估计求得参数，然后新得到的参数值重新被用于E步。。。直至收敛到局部最优解。

  需要注意的一点是：理论中，E步最重要的是求Q函数（李航 9.29式），即完全数据的对数似然函数的期望，M步要做的是用Q函数对各参数求偏导，从而确定满足Q函数最大的各参数的值。而在程序中，因为M步Q函数对参数求偏导后，各参数的值由Q函数中的Gamma(j,k)决定，所以E步只需计算Gamma(j,k)即可，Gamma(j,k)代表第j个样本属于第k个高斯模型的概率，这样M步各参数的值也就确定了，所以切记实例理论与程序实现的区别，好像是程序只是利用了理论的精简结果。

学习完理论之后，最好写个程序练习一下，加深自己的理解。下面的程序(matlab编写)，是这样一个过程：

1.首先产生两组数据点。每一组数据点都是符合给定均值Mu和协方差矩阵∑的。程序和图如下所示，



%产生2个二维正态数据

MU1 = [1 2];

SIGMA1 = [1 0; 0 0.5];

MU2 = [-1 -1];

SIGMA2 = [1 0; 0 1];

X = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];

scatter(X(:,1),X(:,2),10,'*');



  2.对于上述产生的数据，利用EM算法预测这组数据对应的均值μ和协方差矩阵∑

   对于各参数初值的设置：两个均值初始值各设置为两组数据点中的随机一个样本点，权值系数初始值设置为1/k,k是高斯混合模型中模型的个数，协方差矩阵初值设置为所有数据的协方差。

 程序：GMM.m



%程序中Psi对应二维高斯函数，其中的(x-μ)与其转置的顺序与上述提到的高斯函数不同，这样是为了保证矩阵可乘性，不影响结果.

%Gamma 为隐变量的值，Gamma(i,j)代表第i个样本属于第j个模型的概率。

%Mu为期望，Sigma为协方差矩阵%Pi为各模型的权值系数%产生2个二维正态数据

%产生2个二维正态数据

MU1 = [1 2];

SIGMA1 = [1 0; 0 0.5];

MU2 = [-1 -1];

SIGMA2 = [1 0; 0 1];

%生成1000行2列(默认)均值为mu标准差为sigma的正态分布随机数

X = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];

%显示

scatter(X(:,1),X(:,2),10,'.');

%====================

K=2;

[N,D]=size(X);

Gamma=zeros(N,K);

Psi=zeros(N,K);

Mu=zeros(K,D);

LM=zeros(K,D);

Sigma =zeros(D, D, K); 

Pi=zeros(1,K);

%选择随机的两个样本点作为期望迭代初值

Mu(1,:)=X(randint(1,1,[1 N/2]),:);

Mu(2,:)=X(randint(1,1,[N/2 N]),:);

%所有数据的协方差作为协方差初值

for k=1:K

  Pi(k)=1/K;

  Sigma(:, :, k)=cov(X);

end

LMu=Mu;

LSigma=Sigma;

LPi=Pi;

while true

%Estimation Step  

for k = 1:K

  Y = X - repmat(Mu(k,:),N,1);

  Psi(:,k) = (2*pi)^(-D/2)*det(Sigma(:,:,k))^(-1/2)*diag(exp(-1/2*Y/(Sigma(:,:,k))*(Y')));      %Psi第一列代表第一个高斯分布对所有数据的取值

end

Gamma_SUM=zeros(D,D);

for j = 1:N

  for k=1:K

  Gamma(j,k) = Pi(1,k)*Psi(j,k)/sum(Psi(j,:)*Pi');                                               %Psi的第一行分别代表两个高斯分布对第一个数据的取值

  end

end

%Maximization Step

for k = 1:K

%update Mu

  Mu_SUM= zeros(1,D);

  for j=1:N

     Mu_SUM=Mu_SUM+Gamma(j,k)*X(j,:);   

  end

  Mu(k,:)= Mu_SUM/sum(Gamma(:,k));

%update Sigma

 Sigma_SUM= zeros(D,D);

 for j = 1:N

   Sigma_SUM = Sigma_SUM+ Gamma(j,k)*(X(j,:)-Mu(k,:))'*(X(j,:)-Mu(k,:));

 end

 Sigma(:,:,k)= Sigma_SUM/sum(Gamma(:,k));

 %update Pi

 Pi_SUM=0;

 for j=1:N

   Pi_SUM=Pi_SUM+Gamma(j,k);

 end

 Pi(1,k)=Pi_SUM/N;

end

R_Mu=sum(sum(abs(LMu- Mu)));

R_Sigma=sum(sum(abs(LSigma- Sigma)));

R_Pi=sum(sum(abs(LPi- Pi)));

R=R_Mu+R_Sigma+R_Pi;

if ( R<1e-10)

    disp('期望');

    disp(Mu);

    disp('协方差矩阵');

    disp(Sigma);

    disp('权值系数');

    disp(Pi);

    obj=gmdistribution(Mu,Sigma);

    figure,h = ezmesh(@(x,y)pdf(obj,[x,y]),[-8 6], [-8 6]);

    break;

end

 LMu=Mu;

 LSigma=Sigma;

 LPi=Pi;

end

%=====================



%%%%%%%%运行结果%%%%%%%

%

% 期望

%    1.0334    2.0479

%    -0.9703   -0.9879

% 

% 协方差矩阵

% 

% (:,:,1) =

% 

%     0.9604   -0.0177

%    -0.0177    0.4463

% 

% 

% (:,:,2) =

% 

%     0.9976    0.0667

%     0.0667    1.0249

% 

% 权值系数

%     0.4935    0.5065

%%%%%%%%%%%%%%%%%%%%%%%%%

 

参考文章

   http://blog.163.com/baolong_zhu/blog/static/196311091201421185531966/

   《统计学习方法》 李航

   EM算法是一种迭代算法，1977年由Dempster等人总结出，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望(expectation);M步，求极大(maximization).所以这一算法称为期望极大算法(expectation maximization algorithm),简称EM算法。

   首先准备一些预备知识，如：二维高斯混合函数的形式，协方差的定义

高斯混合函数

        对于一维混合高斯函数，如(1)式所示。(所谓的混合高斯，就是将一些高斯函数利用一些权值系数组合起来)  （1）

        二维高斯混合函数，如(2)式所示。

   （2）

       可以看出二维高斯混合模型，其中的d为数据的维度，在这里等于2；∑  代表协方差矩阵，两行两列；这里的x以及μ都是二维的数据，用矩阵表示就是一行两列。这里提到了协方差矩阵，那么什么是协方差矩阵，与方差有什么区别？

协方差矩阵

       理解协方差矩阵的关键在于牢记它计算的是不同维度之间的协方差，而不是不同样本之间。拿到一个样本矩阵，我们最先明确的是一行是一个样本还是一个维度（重要）。

        方差表示的是每个样本值与全体样本值的平均数之差的平方值的平均数，反映了一组数据的离散程度。

 

         既然有了方差还要协方差干啥？显然，方差只能衡量一维数据的离散程度，对于二维数据，就用到了协方差，如下式所示。其中X，Y代表一组数据的两个维度。COV(X,Y)的值大于零，表示两个维度的数据成正相关，小于零则表示成负相关，等于零表示两个维度的数据相互独立。显然当X=Y时，COV(X,X)计算的就是方差。由协方差的定义可以得到COV(X,Y)=COV(Y,X)。

   二维数据的协方差矩阵表示为 ，显然协方差矩阵为对角矩阵。

    三维数据的协方差矩阵表示为，同理知道N维数据的协方差矩阵。

EM算法

   假设一组数据x1,x2,...,xn是来自上述二维高斯混合模型，这里我们利用EM算法确定各个高斯部件的参数，充分拟合给定的数据。算法目的具体的理论知识与算法步骤参考《统计学习方法》李航一书中对于EM算法的讲解(李航的这一本书，绝对值得学习)。

   我主要写一下自己的理解，也为自己以后复习用。

   EM算法其实就是概率模型参数的极大似然估计法，只不过存在一些隐变量，从而不能像一般的概率模型参数估计一样直接求似然函数的最大值，因为这里的似然函数存在隐变量，因而不存在解析解，就不能直接对 似然函数求偏导来求极大值，只能迭代求解。

   基本思想：

   E步，假设参数已知，去估计隐藏变量的期望；

   M步，利用E步求得的隐藏变量的期望，根据最大似然估计求得参数，然后新得到的参数值重新被用于E步。。。直至收敛到局部最优解。

  需要注意的一点是：理论中，E步最重要的是求Q函数（李航 9.29式），即完全数据的对数似然函数的期望，M步要做的是用Q函数对各参数求偏导，从而确定满足Q函数最大的各参数的值。而在程序中，因为M步Q函数对参数求偏导后，各参数的值由Q函数中的Gamma(j,k)决定，所以E步只需计算Gamma(j,k)即可，Gamma(j,k)代表第j个样本属于第k个高斯模型的概率，这样M步各参数的值也就确定了，所以切记实例理论与程序实现的区别，好像是程序只是利用了理论的精简结果。

学习完理论之后，最好写个程序练习一下，加深自己的理解。下面的程序(matlab编写)，是这样一个过程：

1.首先产生两组数据点。每一组数据点都是符合给定均值Mu和协方差矩阵∑的。程序和图如下所示，



%产生2个二维正态数据

MU1 = [1 2];

SIGMA1 = [1 0; 0 0.5];

MU2 = [-1 -1];

SIGMA2 = [1 0; 0 1];

X = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];

scatter(X(:,1),X(:,2),10,'*');



  2.对于上述产生的数据，利用EM算法预测这组数据对应的均值μ和协方差矩阵∑

   对于各参数初值的设置：两个均值初始值各设置为两组数据点中的随机一个样本点，权值系数初始值设置为1/k,k是高斯混合模型中模型的个数，协方差矩阵初值设置为所有数据的协方差。

 程序：GMM.m



%程序中Psi对应二维高斯函数，其中的(x-μ)与其转置的顺序与上述提到的高斯函数不同，这样是为了保证矩阵可乘性，不影响结果.

%Gamma 为隐变量的值，Gamma(i,j)代表第i个样本属于第j个模型的概率。

%Mu为期望，Sigma为协方差矩阵%Pi为各模型的权值系数%产生2个二维正态数据

%产生2个二维正态数据

MU1 = [1 2];

SIGMA1 = [1 0; 0 0.5];

MU2 = [-1 -1];

SIGMA2 = [1 0; 0 1];

%生成1000行2列(默认)均值为mu标准差为sigma的正态分布随机数

X = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];

%显示

scatter(X(:,1),X(:,2),10,'.');

%====================

K=2;

[N,D]=size(X);

Gamma=zeros(N,K);

Psi=zeros(N,K);

Mu=zeros(K,D);

LM=zeros(K,D);

Sigma =zeros(D, D, K); 

Pi=zeros(1,K);

%选择随机的两个样本点作为期望迭代初值

Mu(1,:)=X(randint(1,1,[1 N/2]),:);

Mu(2,:)=X(randint(1,1,[N/2 N]),:);

%所有数据的协方差作为协方差初值

for k=1:K

  Pi(k)=1/K;

  Sigma(:, :, k)=cov(X);

end

LMu=Mu;

LSigma=Sigma;

LPi=Pi;

while true

%Estimation Step  

for k = 1:K

  Y = X - repmat(Mu(k,:),N,1);

  Psi(:,k) = (2*pi)^(-D/2)*det(Sigma(:,:,k))^(-1/2)*diag(exp(-1/2*Y/(Sigma(:,:,k))*(Y')));      %Psi第一列代表第一个高斯分布对所有数据的取值

end

Gamma_SUM=zeros(D,D);

for j = 1:N

  for k=1:K

  Gamma(j,k) = Pi(1,k)*Psi(j,k)/sum(Psi(j,:)*Pi');                                               %Psi的第一行分别代表两个高斯分布对第一个数据的取值

  end

end

%Maximization Step

for k = 1:K

%update Mu

  Mu_SUM= zeros(1,D);

  for j=1:N

     Mu_SUM=Mu_SUM+Gamma(j,k)*X(j,:);   

  end

  Mu(k,:)= Mu_SUM/sum(Gamma(:,k));

%update Sigma

 Sigma_SUM= zeros(D,D);

 for j = 1:N

   Sigma_SUM = Sigma_SUM+ Gamma(j,k)*(X(j,:)-Mu(k,:))'*(X(j,:)-Mu(k,:));

 end

 Sigma(:,:,k)= Sigma_SUM/sum(Gamma(:,k));

 %update Pi

 Pi_SUM=0;

 for j=1:N

   Pi_SUM=Pi_SUM+Gamma(j,k);

 end

 Pi(1,k)=Pi_SUM/N;

end

R_Mu=sum(sum(abs(LMu- Mu)));

R_Sigma=sum(sum(abs(LSigma- Sigma)));

R_Pi=sum(sum(abs(LPi- Pi)));

R=R_Mu+R_Sigma+R_Pi;

if ( R<1e-10)

    disp('期望');

    disp(Mu);

    disp('协方差矩阵');

    disp(Sigma);

    disp('权值系数');

    disp(Pi);

    obj=gmdistribution(Mu,Sigma);

    figure,h = ezmesh(@(x,y)pdf(obj,[x,y]),[-8 6], [-8 6]);

    break;

end

 LMu=Mu;

 LSigma=Sigma;

 LPi=Pi;

end

%=====================



%%%%%%%%运行结果%%%%%%%

%

% 期望

%    1.0334    2.0479

%    -0.9703   -0.9879

% 

% 协方差矩阵

% 

% (:,:,1) =

% 

%     0.9604   -0.0177

%    -0.0177    0.4463

% 

% 

% (:,:,2) =

% 

%     0.9976    0.0667

%     0.0667    1.0249

% 

% 权值系数

%     0.4935    0.5065

%%%%%%%%%%%%%%%%%%%%%%%%%

 

     



     

--------------------- 

作者：lpkinging 

原文：https://blog.csdn.net/u011582757/article/details/70003351

收起

展开全文

机器学习 EM算法

EM算法估计GMM （二维高斯混合模型）参数　ｐython 2.7（西瓜书版）

2019-05-27 09:08:42

#!/usr/bin/python # coding:utf-8 # 19-5-23 上午10:07 # @File : EM.py import numpy as np import math import matplotlib.pyplot as plt ...p = [0.2, 0.3, 0.5] #选择高斯分布的概率 mean1 = [1.0, 2.0] mean...

#!/usr/bin/python
# coding:utf-8 
# 19-5-23 上午10:07
# @File    : EM.py
import numpy as np
import math
import matplotlib.pyplot as plt
#real parm
p = [0.2, 0.3, 0.5] #选择高斯分布的概率
mean1 = [1.0, 2.0]
mean2 = [2.0, 1.0]
mean3 = [1.0, 1.0]
cov1 = [[1.0, 0], [0, 10.0]]
cov2 = [[10.0, 0], [0, 1.0]]
cov3 = [[3.0, 0], [0, 4.0]]
# estimate parm
Emean1 = np.array([1.0, 1.0])

Emean2 = np.array([1.0, 1.0])
Emean3 = np.array([1.0, 1.0])
Ecov1 = np.array([[1.0, 0.0], [0.0, 10.0]])
Ecov2 = np.array([[2.0, 0.0], [0.0, 3.0]])
Ecov3 = np.array([[3.0, 0.0], [0.0, 4.0]])
Ep = p = [0.2, 0.3, 0.5]
# saver parameter
# E 步保存的第l 个模型的　R(n,l)
parameter_dict={}
# 分别存放第l 个模型的　R(n,l)
parameter_dict["gama_list1"] = []
parameter_dict["gama_list2"] = []
parameter_dict["gama_list3"] = []

def DenityOfNorm(Data, Emean,Ecov):
    # return 　EM估计参数　下　data (x, y)处　的正态函数密度
    # Data 输入的数据的点
    #Emean ,EM算法估计的参数
    #Ecov ,EM在第t步估计的协方差矩阵
    #
    # print Ecov.shape
    #
    Ecov_sqrt = math.sqrt(np.linalg.det(Ecov)) #协方差矩阵行列式值的开平方
    Ecov_inv = np.linalg.inv(Ecov) #协方差矩阵的逆
    #we
    Data =np.reshape(Data,[2,1])
    Emean =np.reshape(Emean,[2,1])
    Minus_Data = Data - Emean
    Minus_Data_T = Minus_Data.T
    res = (1.0 / (2.0 * math.pi * Ecov_sqrt)) * math.exp(
        (-0.5) * (np.dot(np.dot(Minus_Data_T, Ecov_inv), Minus_Data)))
    return res
    # sigma_sqrt = math.sqrt(np.linalg.det(sigma))  # 协方差矩阵绝对值的1/2次
    # sigma_inv = np.l

#　生成二维高斯模型的数据
def ini_Normal(mean, cov):
    # 返回一个由size指定形状的数组，数组中的值服从
    # μ = loc, σ = scale
    data =np.random.multivariate_normal(mean, cov,1)
    return data

#生成高斯混合模型的数据
def ini_data(k=1,parm1=0.2,parm2=0.3):
    data= []
    for i in range(k):
        d = np.random.rand(1)
        if d < parm1:
            data.append(ini_Normal(mean1, cov1))
        elif d < parm1+parm2:
            data.append(ini_Normal(mean2, cov2))
        else:
            data.append(ini_Normal(mean3, cov3))
    return data

#第Ｅ步
def E(Data):
    # !!!!!!!!!!!!!!!!
    # 这里一定要清空，　不然会很尴尬，　之前第三个忘了清空，以为陷入局部最优值，检查了很久
    # !!!!!!!!!!!!!!!
    global parameter_dict
    parameter_dict["gama_list1"] = []
    parameter_dict["gama_list2"] = []
    parameter_dict["gama_list3"] = []
    for point in Data:
        # gama_i = (pw * PDF(point, mu_2, sigma_2)) / (
        #         (1.0 - pw) * PDF(point, mu_1, sigma_1) + pw * PDF(point, mu_2, sigma_2))
        # parameter_dict["gama_list"].append(gama_i)
        gam_sum = Ep[0] * DenityOfNorm(point,Emean1,Ecov1)+Ep[1] * DenityOfNorm(point,Emean2,Ecov2)+Ep[2] * DenityOfNorm(point,Emean3,Ecov3)
        gam1_i = Ep[0] * DenityOfNorm(point,Emean1,Ecov1) / gam_sum
        parameter_dict["gama_list1"].append(gam1_i)
        gam2_i = Ep[1] * DenityOfNorm(point,Emean2, Ecov2) / gam_sum
        parameter_dict["gama_list2"].append(gam2_i)
        parameter_dict["gama_list3"].append(1.0-gam2_i-gam1_i)

#第Ｍ步
def M(Data):
    r1 = 0
    r2 = 0
    r3 = 0
    for i in range(len(parameter_dict['gama_list2'])):
        r1 = r1 + parameter_dict["gama_list1"][i]
        r2 = r2 + parameter_dict["gama_list2"][i]
        r3 = r3 + parameter_dict["gama_list3"][i]
        # r3 += 1- parameter_dict["gama_list1"][i]- parameter_dict["gama_list2"][i]
    # updata estimated p
    # for i in range(3):
    Ep[0] = r1 / len(Data)
    Ep[1] = r2 / len(Data)
    Ep[2] = r3 / len(Data)

    # updata estimated mean
    Emean1_new= np.array([0.0, 0.0])
    Emean2_new= np.array([0.0, 0.0])
    Emean3_new= np.array([0.0, 0.0])
    #
    Emean1_new.shape = (1,2)
    Emean2_new.shape = (1,2)
    Emean3_new.shape = (1,2)
    # print parameter_dict["gama_list1"].shape
    for i in range(len(Data)):
        print i
        # print parameter_dict["gama_list1"][i]
        Emean1_new =Emean1_new + Data[i] * parameter_dict["gama_list1"][i] /r1
        Emean2_new =Emean2_new + Data[i] * parameter_dict["gama_list2"][i] /r2
        Emean3_new =Emean3_new + Data[i] * parameter_dict["gama_list3"][i] /r3
        # Emean3_new += Data[i] * (1-parameter_dict["gama_list1"][i]-parameter_dict["gama_list2"])
    global Emean1
    global Emean2
    global Emean3
    # print "---------------------------------------"
    # print Emean1_new.shape
    Emean1 = Emean1_new
    Emean2 = Emean2_new
    Emean3 = Emean3_new


    # updata estimated cov

    Ecov1_new = np.array([[0, 0], [0, 0]])
    Ecov2_new = np.array([[0, 0], [0, 0]])
    Ecov3_new = np.array([[0, 0], [0, 0]])
    for i in range(len(Data)):
        Data_MiuMean1 = Data[i] - Emean1
        Data_MiuMean2 = Data[i] - Emean2
        Data_MiuMean3 = Data[i] - Emean3
        Data_MiuMean1_T = np.array(Data_MiuMean1).T
        Data_MiuMean2_T = np.array(Data_MiuMean2).T
        Data_MiuMean3_T = np.array(Data_MiuMean3).T
        # 测试
        # MYtry = np.dot(Data_MiuMean2,Data_MiuMean2_T)
        # MYtry1 = np.dot(Data_MiuMean2_T,Data_MiuMean2)
        # print MYtry, MYtry.shape
        # print MYtry1, MYtry1.shape
        #

        Ecov1_new = Ecov1_new + np.dot(Data_MiuMean1_T, Data_MiuMean1) * parameter_dict["gama_list1"][i]
        Ecov2_new = Ecov2_new + np.dot(Data_MiuMean2_T, Data_MiuMean2) * parameter_dict["gama_list2"][i]
        Ecov3_new = Ecov3_new + np.dot(Data_MiuMean3_T, Data_MiuMean3) * parameter_dict["gama_list3"][i]
        # Ecov2_new = Ecov2_new + np.d
    global Ecov1
    global Ecov2
    global Ecov3
    # print "------------------------"
    # print Ecov1_new.shape
    Ecov1 = Ecov1_new / r1
    Ecov2 = Ecov2_new / r2
    Ecov3 = Ecov3_new / r3

#不可以控制迭代次数ＥM算法
def EM_iterate( Data, esp=0.0001):


    # Data:数据
    #param esp:终止约束，第n次优化和第ｎ＋１次优化差距小于0.0001

    # set_parameter(mu_1, sigma_1, mu_2, sigma_2, pi_weight)
     while (True):
         old_m1 = Emean1
         old_m2 = Emean2
         old_m3 = Emean3
         E(Data)
         M(Data)
         #
         # print Emean2.shape
         #
         delta_1 = Emean1 - old_m1
         delta_2 = Emean2 - old_m2
         delta_3 = Emean3 - old_m3
         # print delta_1 ,delta_1.shape
         # delta_2 = E - old_m2
         if math.fabs(delta_1[0][0]) < esp and math.fabs(delta_1[0][1]) < esp and math.fabs(
                 delta_2[0][0]) < esp and math.fabs(delta_2[0][1]) < esp and math.fabs(
             delta_3[0][0]) < esp and math.fabs(delta_3[0][1]) < esp:
             break

#EM 算法可以控制迭代次数
def EM_iterat_N(iterm, Data,esp=0.0001):
    #　记录更新的过程
    mean1_trace = [[], []]
    mean2_trace = [[], []]
    mean3_trace = [[], []]
    # mean2_trace = [[], []]
    # mean3_trace = [[], []]
    for i in range(iterm):
            old_m1 = Emean1
            old_m2 = Emean2
            old_m3 = Emean3
            E(Data)
            M(Data)
            #
            # print Emean2.shape
            #
            mean1_trace[0].append(Emean1[0][0])
            mean1_trace[1].append(Emean1[0][1])
            mean2_trace[0].append(Emean2[0][0])
            mean2_trace[1].append(Emean2[0][1])
            mean3_trace[0].append(Emean3[0][0])
            mean3_trace[1].append(Emean3[0][1])
            delta_1 = Emean1 - old_m1
            delta_2 = Emean2 - old_m2
            delta_3 = Emean3 - old_m3
            # print delta_1 ,delta_1.shape
            # delta_2 = E - old_m2
            if math.fabs(delta_1[0][0]) < esp and math.fabs(delta_1[0][1]) < esp and math.fabs(delta_2[0][0]) < esp and math.fabs(delta_2[0][1]) < esp and math.fabs(delta_3[0][0]) < esp and math.fabs(delta_3[0][1]) < esp:
                break
    # print mean1_trace
    # plt.subplot(121)
    # plt.xlim(xmax=5, xmin=2)
    # plt.ylim(ymax=90, ymin=60)
    # plt.xlabel("eruptions")
    # plt.ylabel("waiting")
    # plt.plot(mean1_trace[:,0], mean1_trace[:,1], 'r-')
    # plt.plot(mean1_trace[:,0],mean1_trace[:,1], 'b^')
    #
    # plt.subplot(122)
    # plt.xlim(xmax=4, xmin=0)
    # plt.ylim(ymax=60, ymin=40)
    # plt.xlabel("eruptions")
    # plt.ylabel("waiting")
    # plt.plot(mean2_trace[:,0], mean2_trace[:,1], 'r-')
    # plt.plot(mean2_trace[:,0], mean2_trace[:,1], 'bo')
    # plt.show()
    return mean1_trace,mean2_trace,mean3_trace

def main():
    # data = gen_clusters()
    # save_data(data, '3clusters.txt')
    # d = load_data('3clusters.txt')
    # show_scatter(d)
    # # show_scatter(d)
    data = ini_data(10000)
    #print data
    Emean1.shape = (1,2)
    Emean2.shape = (1,2)
    Emean3.shape = (1,2)
    # EM_iterate( data)
    trace1,trace2,trace3 = EM_iterat_N(50, data)
    # 画图
    plt.subplot(131)
    plt.xlim(xmax=5, xmin=2)
    plt.ylim(ymax=90, ymin=60)
    plt.xlabel("Emean1_x")
    plt.ylabel("Emean1_y")
    plt.plot(trace1[0], trace1[1], 'r-')
    plt.plot(trace1[0], trace1[1], 'b^')
    plt.axis('tight')

    plt.subplot(132)
    plt.xlim(xmax=4, xmin=0)
    plt.ylim(ymax=60, ymin=40)
    plt.xlabel("Emean2_x")
    plt.ylabel("Emean2_y")
    plt.plot(trace2[0], trace2[1], 'r-')
    plt.plot(trace2[0], trace2[1], 'bo')
    plt.axis('tight')

    plt.subplot(133)
    plt.xlim(xmax=4, xmin=0)
    plt.ylim(ymax=60, ymin=40)
    plt.xlabel("Emean3_x")
    plt.ylabel("Emean3_y")
    plt.plot(trace3[0], trace3[1], 'r-')
    plt.plot(trace3[0], trace3[1], 'bv')
    plt.axis('tight')
    plt.show()

    print "real mean1: ", mean1
    print "real mean1: ", mean2
    print "real mean1: ", mean3
    print "estimate means:", Emean1, Emean2, Emean3

if __name__ == '__main__':
    main()












 更新公式：

收起

展开全文

二维高斯分布（Two-dimensional Gaussian distribution）的参数分析
万次阅读 多人点赞 2018-07-12 23:48:29
为此特意回顾了概率论的二维高斯分布的相关概念，并分析了参数对二维高斯分布曲面的影响。 1 多维高斯分布的概率密度函数多维变量X=(x1,x2,...xn)X=(x1,x2,...xn)X = ({x_1},{x_2},...{x_n})的联合概率密度函数...
收起
一种基于YCbCr色彩空间和二维高斯概率分布的方法，将手部区域与背景区分开来，并且绘制测试图像的二维高斯...
2018-06-05 14:57:48
1.以下分割手部区域的代码般自博客：https://blog.csdn.net/soaringlee_fighting/article/details/72983330%% 基于高斯模型的肤色建模与似然图计算 close; clear all; clc; M = [124.2125 132.9449]' ; %肤色...
收起
【Matlab学习手记】二维高斯曲面拟合法定位光斑中心
千次阅读 2019-06-02 21:12:15
终端探测器探测到的光斑灰度分布函数可近似看做高斯分布，因此可以通过二维高斯函数进行拟合，模型表示为：效果图：为了方便计算，做一步变换，两边取对数，得到：展开并进一步变形为：求解...
收起
最小二乘
混合高斯模型算法
2017-12-27 15:46:29
下面介绍一下几种典型的机器算法首先第一种是高斯混合模型算法: 高斯模型有单高斯模型（SGM）和混合高斯模型（GMM）两种...二维情况如下所示：（2）混合高斯模型：对于(b)图所示的情况
收起
混合高斯模型算法（转）
2015-08-25 15:20:00
下面介绍一下几种典型的机器算法首先第一种是高斯混合模型算法: ...二维情况如下所示：（2）混合高斯模型：对于(b)图所示的情况，很明显，单高斯模型是无法解决的。为了解决这个问...
收起
详解高斯混合模型与EM算法
千次阅读 2019-11-05 15:38:19
GSM）一维高斯分布多维高斯分布混合高斯模型（Gaussian mixture model, GMM）混合高斯模型产生的原因直观理解高斯混合模型一维混合高斯模型二维空间3个高斯模型混合极大似然估计(Maximum Likehood Estimate,...
收起
EM算法
二维均值方差高斯分布图 python_图像分割实战 - K均值算法（K-Means）和高斯混合模型（GMM）
2020-11-23 18:46:11
从几何上讲，单高斯分布模型在二维空间应该近似于椭圆，在三维空间上近似于椭球。遗憾的是在很多分类问题中，属于同一类别的样本点并不满足“椭圆”分布的特性。这就引入了高斯混合模型。高斯混合模型GMM：GMM认为...
收起
机器学习混合高斯模型再述
千次阅读 2013-12-03 19:22:55
下面介绍一下几种典型的机器算法首先第一种是高斯混合模型算法: 高斯模型有单高斯模型（SGM）和混合高斯模型（GMM）两种...二维情况如下所示：（2）混合高斯模型：对于(b)图所示的情况
收起
优化机器学习算法计算机视觉迭代
高斯混合模型
2018-09-06 18:24:00
一、什么是高斯混合模型（GMM）　高斯混合模型（Gaussian Mixed Model）指的是多个高斯分布...如果只用一个的二维高斯分布来描述图中的数据。这显然不太合理，毕竟肉眼一看就应该分成两类　但使用两个二维...
收起
二维ising模型概率c语言_伯努利混合模型（Bernoulli Mixture Model）
2020-12-31 08:17:31
一、前言笔者本学期有一门数据挖掘的课程，教授重点讲解了有限混合模型（Finite Mixture Model），是聚类分析常用的方法之一，其中最经典的就是大家所熟知的高斯混合模型（Gaussian Mixture Model，GMM）。...
收起
高斯混合模型（GMM）和EM算法详解
2021-01-20 11:14:56
关于二维高斯分布的参数设定对为高斯曲面的影响，可以参考这篇文章（二维高斯分布的参数分析）（以上两条是基础，为了下面做铺垫，接下来我将通过例子引出高斯混合模型。） 3.高斯混合模型（GMM）为什么会有高斯
收起
基于matlab的高斯模型肤色检测
千次阅读 2020-05-02 08:19:41
摘要：为了实现人体的肤色检测，本文将原始RGB图像转换为... 原理：计算机采集的图像大都为RGB图像，研究发现，肤色图像在YCbCr空间中，以Cb、Cr为维度的二维空间，肤色图像符合高斯模型，故我们首先把RGB图像转换...
收起
机器学习——概率分类(三)高斯概率密度与混合高斯模型
2020-02-10 11:18:12
机器学习——概率分类(三)高斯概率密度与混合高斯模型 在之前的文章机器学习——概率分类(一)朴素贝叶斯模型一文中，我们主要介绍了机器学习中的概率分类问题。我们提出了简单的朴素贝叶斯模型来进行概率密度的估计...
收起
基于混合高斯模型的物体成分拟合方法
2020-10-16 18:47:44
利用二维高斯混合函数，用高斯成分来拟合物体的边缘图像，使得物体的表征由单一的像素表示转变为利用成分进行表征的方式。为了使得拟合结果更具有健壮性，在算法中还引入了分裂-归约机制来对拟合结果进行修正。实验...
收起
高斯混合模型算法
千次阅读 2012-06-21 15:21:52
下面介绍一下几种典型的机器算法首先第一种是高斯混合模型算法: ...二维情况如下所示：（2）混合高斯模型：对于(b)图所示的情况，很明显，单高斯模型是无法解决的。为了解决这个问
收起
算法 parameters matrix function each
高斯光束对二维导体平板太赫兹雷达散射截面影响
2021-02-25 06:51:00
为研究这种假设给平板RCS计算结果带来的误差，计算了高斯光束入射条件下二维导体平板的RCS。仿真研究了入射光为2.52 THz激光对导体平板雷达散射截面的影响，同时与平面波入射结果进行了比较分析。计算结果显示，当...
收起
高斯混合模型算法:
2015-05-11 15:20:32
下面介绍一下几种典型的机器算法首先第一种是高斯混合模型算法: 高斯模型有单高斯模型（SGM）和混合高斯模型（GMM）两种。...二维情况如下所示：（2）混合高斯模型：对于(b)图所示的
收起
机器学习（四）高斯混合模型
千次阅读 2015-04-24 11:46:37
开始学习高斯混合模型，需要先简单复习一下单高斯模型的参数估计方法，描述一个高斯模型其实就是要计算它的均值、协方差矩阵(一维空间为方差，二维以上称之为协方差矩阵)：假设有数据集X={x1,x2,x3……，xn},那么...
收起
机器学习

二维高斯混合模型GMM图形化简单明了

EM算法在二维高斯混合模型参数估计中的应用

高斯混合模型

多维高斯混合模型

二维高斯混合模型参数估计的EM算法

实验步骤

考虑采样点数N对参数估计的影响

EM算法--二维高斯混合模型(GMM)

EM算法--二维高斯混合模型(GMMs)

EM算法估计GMM （二维高斯混合模型）参数　ｐython 2.7（西瓜书版）

二维高斯分布（Two-dimensional Gaussian distribution）的参数分析

一种基于YCbCr色彩空间和二维高斯概率分布的方法，将手部区域与背景区分开来，并且绘制测试图像的二维高斯...

【Matlab学习手记】二维高斯曲面拟合法定位光斑中心

混合高斯模型算法

混合高斯模型算法（转）

详解高斯混合模型与EM算法

二维均值方差高斯分布图 python_图像分割实战 - K均值算法（K-Means）和高斯混合模型（GMM）

机器学习混合高斯模型再述

高斯混合模型

二维ising模型概率c语言_伯努利混合模型（Bernoulli Mixture Model）

高斯混合模型（GMM）和EM算法详解

基于matlab的高斯模型肤色检测

机器学习——概率分类(三)高斯概率密度与混合高斯模型

基于混合高斯模型的物体成分拟合方法

高斯混合模型算法

高斯光束对二维导体平板太赫兹雷达散射截面影响

高斯混合模型算法:

机器学习（四）高斯混合模型

二维高斯模型

二维高斯混合模型GMM图形化简单明了

EM算法在二维高斯混合模型参数估计中的应用

高斯混合模型

多维高斯混合模型

二维高斯混合模型参数估计的EM算法

实验步骤

考虑采样点数N对参数估计的影响

EM算法--二维高斯混合模型(GMM)

EM算法--二维高斯混合模型(GMMs)

EM算法估计GMM （二维高斯混合模型）参数 ｐython 2.7（西瓜书版）

二维高斯分布（Two-dimensional Gaussian distribution）的参数分析

一种基于YCbCr色彩空间和二维高斯概率分布的方法，将手部区域与背景区分开来，并且绘制测试图像的二维高斯...

【Matlab学习手记】二维高斯曲面拟合法定位光斑中心

混合高斯模型算法

混合高斯模型算法（转）

详解高斯混合模型与EM算法

二维 均值方差高斯分布图 python_图像分割实战 - K均值算法（K-Means）和高斯混合模型（GMM）

机器学习 混合高斯模型再述

高斯混合模型

二维ising模型概率c语言_伯努利混合模型（Bernoulli Mixture Model）

高斯混合模型（GMM）和EM算法详解

基于matlab的高斯模型肤色检测

机器学习——概率分类(三)高斯概率密度与混合高斯模型

基于混合高斯模型的物体成分拟合方法

高斯混合模型算法

高斯光束对二维导体平板太赫兹雷达散射截面影响

高斯混合模型算法:

机器学习（四）高斯混合模型

二维高斯模型

EM算法估计GMM （二维高斯混合模型）参数　ｐython 2.7（西瓜书版）

二维均值方差高斯分布图 python_图像分割实战 - K均值算法（K-Means）和高斯混合模型（GMM）

机器学习混合高斯模型再述