精华内容
下载资源
问答
  • MATLAB陡度代码gmm1-用于MATLAB的一维高斯混合模型工具箱 一维高斯混合模型(gmm1)的计算工具箱。 该代码通常很快,但是还有进一步改进的空间(例如,改进的矢量化)。 内容: gmm1cdf.m :gmm1累积分布函数(cdf)...
  • phais = [1.0/k for i in xrange(k)] # 各高斯模型的系数 mus = [i for i in xrange(k)] # 混合高斯的初始均值 sigmas = [1 for i in xrange(k)] # 混合高斯的初始标准差 # 开始学习 for i in xrange(100): Qs...

    https://blog.csdn.net/u012436149/article/details/53557008

    # -*- coding: utf-8 -*-
    # for multi-Gaussian
    __author__ = "KeithYin"
    
    import numpy as np
    def gaussian(x,mu,sigma):
        temp = -np.square(x-mu)/(2*sigma)
        return np.exp(temp)/(np.sqrt(2.0*np.pi*sigma)) # sigma = sigma^2
    def e_step(data, phais, mus, sigmas):
        Qs = []
        for i in xrange(len(data)):
            q = [phai*gaussian(data[i],mu,sigma) for phai,mu,sigma in zip(phais,mus,sigmas)]
            Qs.append(q)
        Qs = np.array(Qs)
        Qs = Qs / np.sum(Qs,axis=1).reshape(-1,1)
        return Qs
    def m_step(data, phais, mus, sigmas, Qs):
        data = np.array(data)
        gama_j = np.sum(Qs,axis=0)
        new_phais = gama_j/len(data)
        print "new_phai:",
        print new_phais
        mu_temp = np.sum(Qs*(data.reshape(-1,1)),axis=0)
        new_mus =mu_temp/gama_j
        X_i_mu_j = np.square(np.array([data]).reshape(-1,1)-np.array([mus]))
        new_sigmas = np.sum(Qs*X_i_mu_j,axis=0)/gama_j
        return new_phais,new_mus,new_sigmas
    def EM(data,k=1):
        # 设置均值
        phais = [1.0/k for i in xrange(k)] # 各高斯模型的系数
        mus = [i for i in xrange(k)]       # 混合高斯的初始均值
        sigmas = [1 for i in xrange(k)]    # 混合高斯的初始标准差
        # 开始学习
        for i in xrange(100):
            Qs = e_step(data,phais,mus,sigmas)
            phais, mus, sigmas= m_step(data,phais,mus,sigmas,Qs)
            print phais,mus,sigmas
    
    if __name__ == "__main__":
        s1 = np.random.normal(19,1,10000)
        s2 = np.random.normal(1,1,10000)
        s3 = np.concatenate((s1,s2),axis=0)
    
        EM(s3,2)
    展开全文
  • 高斯混合模型意味着每个数据点(随机)从 C 类数据之中抽取,概率 p_i 从第 i 类中抽取,并且每个类都分布为具有平均标准差 mu_i 和 sigma_i 的高斯分布。 给定从这种分布中提取的组数据,我们试图估计这些未知...
  • EM最大期望算法是个数值求解似然函数极大值的迭代算法,就好像梯度下降算法是种数值求解损失函数极小值的迭代算法一样。EM算法通常适合于随机变量依赖于另外一些不可观测的随机变量(称之为隐...

    EM最大期望算法是一个数值求解似然函数极大值的迭代算法,就好像梯度下降算法是一种数值求解损失函数极小值的迭代算法一样。

    EM算法通常适合于随机变量依赖于另外一些不可观测的随机变量(称之为隐含变量或者中间变量)的场景

    此时由于似然函数的表示形式较为复杂(含有对隐含变量的累加求和或者积分),难以求导获取似然函数的极大值,也无法方便地应用梯度下降算法进行优化。

    EM算法是一个类似梯度下降算法的迭代算法,它首先给随机变量分布参数赋初始值,然后寻找到了一个便于优化的似然函数的下界 (恰好为似然函数在某个分布下的期望Expectation,期望中消去了隐变量),并通过不断地优化(Maximization) 这个下界求解似然函数的极值

    EM算法在机器学习的许多算法中都有使用到,如

    • KMeans:实际上K-Means是一种Hard EM算法, 隐变量直接取最大概率的位置。

    • 支持向量机的SMO算法

    • LDA主题模型参数估计

    • 混合高斯模型的参数估计

    • HMM隐马尔科夫模型的参数估计

    本篇文章我们将详述EM算法的推导过程,并以一维GMM高斯混合模型为例,示范EM算法的应用方法。

    公众号后台回复关键字:源码,获取本文含全部公式的markdown文件。

    一,EM最大期望算法

    当我们关心的随机变量依赖于另外一些不可观测的随机变量时,通过对我们关心的随机变量采样,我们将难以直接通过最大似然估计的方法推断我们关心的随机变量分布律中的未知参数。

    例如我们有100个学生的身高数据,其中有一些是男生,另外一些是女生。男生和女生的身高服从不同的正态分布,但是我们不知道哪些数据是男生的,哪些是女生的,这是这个模型的隐含变量,是不可以被观测到的。

    那么如何根据这批身高数据估计男女生各自正态分布的均值和方差,以及这批数据中男生的比例呢?

    期望最大化算法 EM (Expectation Maximization)这时候就派上用场了,它能够解决这种含有隐含随机变量的模型的参数估计的问题。

    设观测随机变量为, 隐含随机变量为,待确定参数为

    确定时,的分布函数由给出。

    按照极大似然原理,并使用全概率公式,似然函数可以写成

    对数似然函数可以写成

    对数似然函数中,由于有对的求和,如果尝试对求偏导等于0来计算最优的,将难以得到对应的解析解。这和目标函数非常复杂时,无法直接解析求解只能使用梯度下降这类迭代算法是一样的。

    从原则上说,在一些较为简单的情况下我们也能够使用梯度下降法求解对数似然的最优值,例如当隐藏变量Z是离散随机变量时,且可取值较少,我们很容易将对z的求和表示出来,从而可以计算梯度进而使用梯度下降法。

    但对于一般情况,对z的求和将难以进行,如果Z是连续的随机变量,对z的求和将变成积分,此时使用梯度下降法将更加困难。

    我们可以尝试和梯度下降算法效果相当的迭代算法。最大期望算法EM正是可以实现这个目的。

    大概原理如下,我们首先给赋初始值 ,然后在此基础上,找到一个可以使得对数似然函数变大的,然后再在此基础上找到一个能够使对数似然函数变得更大的,如此便可不断地提高对数似然函数的值。迭代执行n干次后,如果的差值足够小,那么我们认为就找到了比较合适的  作为 的估计值。

    下面阐述最大期望算法的原理推导。

    假设在第n次迭代,我们的对数似然函数取值为

    我们希望找到一个使得

    下面我们开始寻找符合条件的

    构造函数

    由于是严格凹函数,Jensen不等式成立

    存在以下缩放:

    注意到

    因此

    取 

    则有

    即符合我们寻找的条件。

    消去无关量,我们可以得到

    注意到实际上是一个分布,因此右边可以理解成求随机变量

    分布下期望的最大值。

    总结下 EM算法算法的流程:

    (1) 初始化

    注意这里的模型参数要是完备的,即给定这些参数,能够计算联合概率分布函数

    ,对于男女生混合身高的例子,我们的参数应当包括,即男生平均身高和身高标准差,女生平均身高和身高标准差,以及男生的比例。

    (2) 计算E步,又分成两小步,先计算概率分布,再算出期望

    (3) 对E求极大,解出的新的估计,将新的估计值代入第(1)步,直到收敛。

    可以证明EM算法是收敛的,但不能保证它能收敛到全局最优,因此可以尝试多个不同的初始值,计算结果,并挑选能够使似然函数取值最大的结果。

    二,一维GMM高斯混合模型

    高斯分布模型也叫正态分布模型,其概率密度函数如下:

    GMM高斯混合模型的概率密度函数为多个高斯分布的线性组合:

    其中为正数,并且:

    高斯混合模型的参数可以理解为样本属于第类的概率。

    则高斯混合模型的概率密度函数可以表示成如下形式:

    根据EM算法,

    (1)我们首先取初始化参数

    然后执行迭代过程。

    (2)我们先求期望值

    代入贝叶斯公式

    (3)我们求期望极大值对应的 作为  

    考虑到约束

    根据拉格朗日乘子法,作拉格朗日函数

    取极大值时我们有:

    于是我们有:

    解得:

    如此迭代,直到收敛。

    展开全文
  • 它开箱即用,生成一维高斯混合的随机数据集,并将推理过程可视化。 参考: * Carl Edward Rasmussen:无限高斯混合模型 -> http://www.kyb.mpg.de/fileadmin/user_upload/files/publications/pdfs/pdf2299.pdf
  • 它开箱即用,生成二维高斯混合的随机数据集,并使用无限高斯混合模型可视化推理过程。 参考: * Carl Edward Rasmussen:无限高斯混合模型 -> ...
  • 高斯混合模型 高斯混合模型是指具有如下形式的概率分布模型: 其中,是系数,,;是高斯分布密度,, ...称为第k个分模型。...对于二维高斯混合模型,d=2,y和都是二维的数据,用矩阵表示就是行...
    • 高斯混合模型

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

    P(y|\theta )=\sum_{k=1}^{K}\alpha _{k}\phi (y|\theta _{k})

    其中,\alpha _{k}是系数,\alpha _{k}\geqslant 0\sum_{k=1}^{K}\alpha _{k}=1\phi (y|\theta _{k})是高斯分布密度,\theta _{k}=(\mu _{k},\sigma _{k}^{2})

    \phi (y|\theta _{k})=\frac{1}{\sqrt{2\pi }\sigma _{k}}exp(-\frac{(y-\mu _{k})^{2}}{2\sigma _{k}^{2}})

    称为第k个分模型。

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

    • 多维高斯混合模型

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

    P(y)=\sum_{k=1}^{K}\alpha _{k}\frac{1}{(2\pi )^{d/2}\left | \Sigma \right |^{1/2}}exp(-\frac{1}{2}(y-\mu _{k})^{T}\Sigma ^{-1}(y-\mu _{k}))

    其中d为数据的维度,\mu为均值,\Sigma为协方差矩阵。

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

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

    输入:观测数据y_{1},y_{2},...,y_{N},二维高斯混合模型;

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

    算法步骤如下:

    (1)取参数的初始值开始迭代

    (2)E步:依据当前模型参数,计算分模型k对观测数据y_{j}的响应度:

    \widehat{\gamma }_{jk}=\frac{\alpha _{k}\phi (y_{j}|\theta _{k})}{\sum_{k=1}^{K}\alpha _{k}\phi (y_{j}|\theta _{k})},\; \; \; j=1,2,...,N;\; k=1,2,...,K

    (3)M步:计算新一轮迭代的模型参数:

    \widehat{\mu }_{k}=\frac{\sum_{j=1}^{N}\widehat{\gamma }_{jk}y_{j}}{\sum_{j=1}^{N}\widehat{\gamma }_{jk}},\; \; \; k=1,2,...,K

    \widehat{\Sigma }_{k}=\frac{\sum_{j=1}^{N}\widehat{\gamma }_{jk}(y_{j}-\mu _{k})^{T}(y_{j}-\mu _{k})}{\sum_{j=1}^{N}\widehat{\gamma }_{jk}},\; \; \; k=1,2,...,K

    \widehat{\alpha }_{k}=\frac{\sum_{j=1}^{N}\widehat{\gamma }_{jk}}{N},\; \; \; k=1,2,...,K

    (4)重复第(2)步和第(3)步,直到收敛。

    • 实验步骤

    1.确定二维高斯混合模型的参数

    我们的二维高斯混合模型为两个权重不同的二维单高斯分布组合而成,所以\alpha _{k}可以简化为一个参数,即第一个二维单高斯分布的权重为π,第二个二维单高斯分布的权重为1-π。选择π值为0.8,第一个二维单高斯分布的均值为[0.0, 0.0],协方差矩阵为\begin{bmatrix} 1.0 &0.0 \\ 0.0 &1.0 \end{bmatrix},第二个二维单高斯分布的均值为\left [ 3.0,3.0 \right ],协方差矩阵为\begin{bmatrix} 1.0 &0.5 \\ 0.5 &1.0 \end{bmatrix}

    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,第一个均值为[0.0, 0.0],第一个协方差矩阵为\begin{bmatrix} 1.0 &0.0 \\ 0.0 & 1.0 \end{bmatrix},第二个均值为[1.0,1.0],第二个协方差矩阵为\begin{bmatrix} 1.0 &0.0 \\ 0.0 &1.0 \end{bmatrix}

    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对观测数据y_{j}的响应度,由于该模型只有两个二维单高斯分布组成,故K=2:

    \widehat{\gamma }_{jk}=\frac{\alpha _{k}\phi (y_{j}|\theta _{k})}{\sum_{k=1}^{K}\alpha _{k}\phi (y_{j}|\theta _{k})},\; \; \; j=1,2,...,N;\; k=1,2

    其中\phi (y_{j}|\theta _{k})的表达式为:

    \phi (y_{j}|\theta _{k})=\frac{1}{(2\pi )^{d/2}\left | \Sigma \right |^{1/2}}exp(-\frac{1}{2}(y-\mu _{k})^{T}\Sigma ^{-1}(y-\mu _{k}))

        # 高斯分布密度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步

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

    \widehat{\mu }_{k}=\frac{\sum_{j=1}^{N}\widehat{\gamma }_{jk}y_{j}}{\sum_{j=1}^{N}\widehat{\gamma }_{jk}},\; \; \; k=1,2

        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

    \widehat{\Sigma }_{k}=\frac{\sum_{j=1}^{N}\widehat{\gamma }_{jk}(y_{j}-\mu _{k})^{T}(y_{j}-\mu _{k})}{\sum_{j=1}^{N}\widehat{\gamma }_{jk}},\; \; \; k=1,2

        # 模型的协方差矩阵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

    \widehat{\alpha }_{k}=\frac{\sum_{j=1}^{N}\widehat{\gamma }_{jk}}{N},\; \; \; k=1,2

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

    即只需要求下式:

    \widehat{\pi }=\frac{\sum_{j=1}^{N}\widehat{\gamma }_{j0}}{N}

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

    (4)重复第2步和第3步,直到收敛

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

        # 权重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算法学习高斯混合模型的参数效果越好。

    展开全文
  • 可以看出二维高斯混合模型,其中的d为数据的维度,在这里等于2;∑ 代表协方差矩阵,两行两列;这里的x以及μ都是二维的数据,用矩阵表示就是行两列。 协方差矩阵 理解协方差矩阵的关键在于牢记它计算的是不同...

     

    点赞支持下吧,让更多的人也能看到这篇内容(收藏不点赞,都是耍流氓 -_-)

    目录

    高斯混合函数

    协方差矩阵

    什么是高斯混合模型(GMM)

    GMM原理

    高斯混合聚类模型伪代码

    MATLAB实现二维高斯混合模型分类的代码,并生成gif图描述收敛过程

    资源传送门

    「❤️ 感谢大家」


    高斯混合函数

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

    img (2)

    可以看出二维高斯混合模型,其中的d为数据的维度,在这里等于2;∑ 代表协方差矩阵,两行两列;这里的x以及μ都是二维的数据,用矩阵表示就是一行两列。

    协方差矩阵

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

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

    img

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

    img

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

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

    什么是高斯混合模型(GMM)

     高斯混合模型(Gaussian Mixed Model)指的是多个高斯分布函数的线性组合,通常用于解决同一集合下的数据包含多个不同的分布的情况,如解决分类情况

     如下图,明显分成两个聚类。这两个聚类中的点分别通过两个不同的正态分布随机生成而来。如果只用一个的二维高斯分布来描述图中的数据。这显然不太合理,毕竟肉眼一看就应该分成两类

    img

     

     但使用两个二维高斯分布来描述图中的数据,将两个二维高斯分布N(\mu_1,\sum_1),N(\mu_2,\sum_2)做线性组合,用线性组合后的分布来描述整个集合中的数据。这就是高斯混合模型(GMM)

    img

    GMM原理

     设有随机变量XX,则混合高斯模型可以这样表示:p(x) = \sum_{k=1}^{K}\pi_kN(x|\mu_k,\sum_k),其中N(x|\mu_k,\sum_k)称为混合模型中的第k个分量(component)。如前面图中的例子,有两个聚类,可以用两个二维高斯分布来表示,那么分量数K=2,\pi_k是混合系数(mixture coefficient),且满足:\sum_{k=1}^{K}\pi_k,0\leq\pi_k\leq 1,其中πk相当于每个高斯分布的权重

     

    GMM用于聚类时,假设数据服从混合高斯分布(Mixture Gaussian Distribution),例如上图的例子,很明显有两个聚类,可以定义K=2,那么对应的GMM形式如下:

    p(x) = \pi_1N(x|\mu_1,\sum_1) + \pi_2N(x|\mu_2,\sum_2)

     上式中未知参数有六个:(\pi_1,\mu_1,\sum_1, \pi_2,\mu_2,\sum_2),GMM聚类时分为两步,第一步是随机地在这K个分量中选一个,每个分量被选中的概率即为混合系数πk,可以设定π1=π2=0.5,表示每个分量被选中的概率是0.5,即从中抽出一个点,这个点属于第一类的概率和第二类的概率各占一半。但实际应用中事先指定πkπk的值是很笨的做法:当从图中的集合随机选取一个点,怎么知道这个点是来自N(x|\mu_1,\sum_1)还是N(x|\mu_2,\sum_2)呢?

     换言之怎么根据数据自动确定π1,π2的值?这就是GMM参数估计的问题。要解决这个问题,可以使用EM算法。通过EM算法,我们可以迭代计算出GMM中的参数:\pi_k,\mu_k,\sum_k

     

     我们目的是找到这样一组参数,它所确定的概率分布生成这些给定的数据点的概率最大,而这个概率实际上就等于:\prod_{i=1}^{N}p(x_i),我们把这个乘积称作似然函数 (Likelihood Function)。我们通常会对其取对数,把乘积变成加和:\sum_{i=1}^{N}logp(x_i),得到 log-likelihood function 。接下来将这个函数最大化(通常的做法是求导并令导数等于零,然后解方程),我们就认为这是最合适的参数,这样就完成了参数估计的过程。下面让我们来看一看 GMM 的 log-likelihood function :

    \sum_{i=1}^{N}log\{\sum_{i=1}^{K}\pi_kN(x_i|\mu_k,\sum_k)\}

     由于上式太复杂,没法直接用求导求得最大值。为了解决这个问题,我们采取之前从 GMM 中随机选点的办法:分成两步

      1. 估计每个数据由每个 Component 生成的概率(并不是每个 Component 被选中的概率):对于每个数据xixi来说,它由第kk个 Component 生成的概率为:

    img

       由于式子里的\mu_k,\sum_k也是需要我们估计的值,我们采用迭代法,在计算\gamma(i,k)的时候我们假定\mu_k,\sum_k均已知,我们将取上一次迭代所得的值(或者初始值)

     

      2. 估计每个Component 的参数:现在我们假设上一步中得到的γ(i,k)就是正确的"数据xi由Component kk生成的概率",亦可以当做该Component 在生成这个数据上所做的贡献,或者说,我们可以看作xi这个数据其中有γ(i,k)xi这部分是由Component k所生成的,集中考虑所有的数据点,现在实际上可以看作Component 生成了\gamma(1,k)x_1,...,\gamma(N,k)x_N这些点,由于每个 Component 都是一个标准的 Gaussian 分布,可以很容易求出最大似然所对应的参数值:

    img

       其中N_k = \sum_{i=1}^{N}\gamma(i,k)并且πk也顺理成章地可以估计N_k/N

     

      3. 重复迭代前面两步,直到似然函数的值收敛为止

     

    高斯混合聚类模型伪代码

    输入:样本集D ={x1,x2,....xm}

        高斯混合成分个数k

    过程

        初始化高斯混合分布的模型参数{(ai,ui,∑i)} 三个参数含义:分别为混合系数,均值向量,协方差矩阵

        迭代:

        for j = 1,2,......,m do

            根据img计算xj的各混合成分生成的后验概率

            img pm等于上面那个式子

       end for

       for i =1,2,....k do

           计算新均值向量 img

           计算新协方差矩阵 img

           计算新混合系数 img

       end for

       将模型参数{(ai,ui,∑i)|1<=i<=k}更新为{(ai’,ui’,∑i’)|1<=i<=k}

    until 满足停止条件

    Ci = Ø (1<=i<=k)

    for j =1,2,...m do

        根据img(i=1,2,....k) 确定xj的簇标记λj (确定属于哪个i)

        将xj划分相应的簇:

    end for

    输出:簇划分c = {c1,c2.....ck}

    MATLAB实现二维高斯混合模型分类的代码,并生成gif图描述收敛过程

    % 高斯混合模型EM 算法参数估计示例
    % 输入 -----------------------------------------------------------------
    % data: N x D , N 个 dim 维的数据矩阵
    % Alpha: 1 x M  类的先验概率
    % Mu: M x D , M 个 Gauss分布的均值向量
    % Sigma: D x D x M , M 个Gauss分布的协方差矩阵
    % 输出 ----------------------------------------------------------------
    clc;
    clear;
    close all;
    %% 生成随机数据
    N=1000;
    alpha0 = [0.3 0.3 0.4];
    mu = [ -3 4; 0 0; 3 4];
    Sigma(:,:,1) = [ 2 0; 0 1];
    Sigma(:,:,2) = [ 2 0; 0 1];
    Sigma(:,:,3) = [ 2 0; 0 3];
    M = length(alpha0);
    R = mnrnd(N, alpha0);
    % Z 类别
    Z = [];
    for j = 1:M
        Z = [Z; ones([R(j) 1])*j];
    end
    data = zeros([N 2]);
    for j = 1:M
        data(Z == j,:) = mvnrnd(mu(j,:), Sigma(:,:,j), R(j));
    end
    %% EM算法
    prev_mu = [ -5 -1; -1 -1; 4 -1];
    %prev_mu = [ -3 0; 0 7; 6 0];
    %prev_mu=3*randn(3,2);
    prev_Sigma(:,:,1) = [ 1 0; 0 1];
    prev_Sigma(:,:,2) = [ 1 0; 0 1];
    prev_Sigma(:,:,3) = [ 1 0; 0 1];
    prev_alpha = ones([M 1])/M;
    MaxIter = 50;
    pxi = zeros([N M]);
    next_mu = zeros(size(prev_mu));
    next_Sigma = zeros(size(prev_Sigma));
    d = size(data,2);
    for t = 1:MaxIter
        if t > 1
            % Compute the porier probality
            for l = 1:M
                pxi(:,l) = mvnpdf(data,prev_mu(l,:), prev_Sigma(:,:,l));
            end
            pxi_sum = sum(pxi,2);
            pxi = pxi ./ repmat(pxi_sum,[1 3]);
            % update alpha^(i+1)
            next_alpha = sum(pxi)'/N;
            % update mu^(i+1)
            for l = 1:M
                next_mu(l,:) = sum(data .* repmat(pxi(:,l),[1 d])) / sum(pxi(:,l));
            end
            % update Sigma
            for l = 1:M
                zero_mean_data = data - repmat(next_mu(l,:),[N 1]);
                covariances = zeros([d d N]);
                for i = 1:N
                    covariances(:,:,i) = zero_mean_data(i,:)' * zero_mean_data(i,:) * pxi(i,l);
                end
                next_Sigma(:,:,l) = sum(covariances,3)/sum(pxi(:,l));
            end
            prev_mu     = next_mu;
            prev_Sigma  = next_Sigma;
            prev_alpha  = next_alpha;
        else
            next_mu     = prev_mu;
            next_Sigma  = prev_Sigma;
            next_alpha  = prev_alpha;
        end
        plot(data(Z==1,1), data(Z==1,2), 'o');
        hold on;
        plot(data(Z==2,1), data(Z==2,2), 'o');
        plot(data(Z==3,1), data(Z==3,2), 'o');
        axis([min(data(:,1)) max(data(:,1)) min(data(:,2)) max(data(:,2))]);
        % 生成高斯对象
        gmm_pdf = gmdistribution(next_mu,next_Sigma,next_alpha);
        % 高斯混合分布的函数
        ezcontour(@(u,v)pdf(gmm_pdf,[u v]));
        axis([min(data(:,1)) max(data(:,1)) min(data(:,2)) max(data(:,2))]);
        hold off;
        % gif 动图
        [em, map]= rgb2ind(frame2im(getframe), 256);
        if t == 1
            imwrite(em, map, 'em.gif', 'gif', 'Loopcount',inf, 'DelayTime',2)
        else
            imwrite(em, map, 'em.gif', 'gif', 'writeMode','append','DelayTime',0.2)
        end
        % waitforbuttonpress;
    end

    运行结果:

    高斯混合模型:EM 算法参数估计示例

    资源传送门

    • 关注【做一个柔情的程序猿】公众号
    • 在【做一个柔情的程序猿】公众号后台回复 【python资料】【2020秋招】 即可获取相应的惊喜哦!

    「❤️ 感谢大家」

     

    • 点赞支持下吧,让更多的人也能看到这篇内容(收藏不点赞,都是耍流氓 -_-)
    • 欢迎在留言区与我分享你的想法,也欢迎你在留言区记录你的思考过程
    展开全文
  • 使用基于高斯混合模型的机器学习对 2 类和 3 类问题进行一维矩阵分类。 它还包含一个基于矩阵的 AND 门示例和大小为 12 和 3 的输入样本
  • EM算法--二维高斯混合模型(GMM)

    万次阅读 2017-04-10 19:19:16
    EM算法是种迭代算法,1977年由Dempster等人总结出,用于含有隐变量的概率模型参数的极大似然估计,或极大后验概率估计。EM算法的每次迭代由两步组成:E步,求期望(expectation);M步,求极大(maximization).所以这...
  • 该算法基于RGB空间建立码本模型,然后基于码字中的R、G、B分量建立三维高斯模型,从而使整个码本具有三维高斯混合模型的特征.实验结果表明:该算法具有较高的实时性(该算法的平均帧率约23.0帧/s,而iGMM( improved ...
  • EM算法--二维高斯混合模型(GMMs)

    千次阅读 2019-03-30 20:58:21
    参考文章 ...EM算法是种迭代算法,1977年由Dempster等人总结出,用于含有隐变量的概率模型参数的极大似然估计,或极大后验概率估计。EM算法的每次迭代由两步组成:E步,求期望(expectation);M步...
  • 高斯混合模型

    2021-05-09 17:13:53
    先给出高斯混合模型的定义,高斯混合模型是指具有如下形式的概率分布模型: (1) 其中,是系数,且,,而是高斯分布密度,,对于随机变量y是一维数据时, (2) 称为第k个分模型。 理解 高斯混合模型属于生成...
  • 详解高斯混合模型与EM算法

    千次阅读 2019-11-05 15:38:19
    详解高斯混合模型与EM算法 详解高斯混合模型与EM算法高斯混合模型单高斯模型(Gaussian single model, GSM)一维高斯分布多维高斯分布混合高斯模型(Gaussian mixture model, GMM)混合高斯模型产生的原因直观理解...
  • 高斯混合模型的和详细的EM算法推导见《统计学习方法》 这里说明一点: EM算法叫期望极大算法,是先在当前参数下求得完全分布对于隐变量的期望,然后求解对数似然的最大化问题,以获得新轮迭代的参数。 其中核心...
  • 一维高斯模型(One-dimensional Gaussian Model) 若随机变量X服从一个数学期望为,标准方差为的高斯分布,记为: x~N(,)。 则概率密度函数为: 高斯分布的期望值决定了其位置,标准方差决定了其幅度。 ...
  • 小白枚,接触到GMM和EM,现将学习到的整理出来,如有错误,欢迎指正,文中涉及到公式的推导比较...(以上两条是基础,为了下面做铺垫,接下来我将通过例子引出高斯混合模型。) 3.高斯混合模型(GMM) 为什么会有高斯
  • 一维高斯分布的概率密度函数如下: 多维变量X=(x1,x2,…xn)的联合概率密度函数为: 这里引用李航《统计学习方法》书中的定义 简而言之,GMM是多个高斯分布的加权和,并且权重α之和等于1 。 Sklearn ...
  • 高斯混合模型算法

    千次阅读 2012-06-21 15:21:52
    首先第种是高斯混合模型算法: 高斯模型有单高斯模型(SGM)和混合高斯模型(GMM)两种。 (1)单高斯模型: 为简单起见,阈值t的选取一般靠经验值来设定。通常意义下,我们一般取t=0.7-0.75之间。 二...
  • 新颖检测中,可应用高斯混合模型建立已知数据模型,拟合数据分布,但当数据数较高时,自由参数太多,训练需要巨大的数据采样,而ICA搜寻数据的最大统计独立表示,可以将数据从高维空间投影到低空间。提出种...
  • 篇文章学习了极大似然估计和EM算法,这章节就是学习EM算法的具体应用场景,高斯混合模型,是一中聚类的算法。 高斯混合模型(GMM) 我们将个分布复杂的数据分布,用多个高斯概率分布来共同表示,如上图...
  • 高斯混合模型算法:

    2015-05-11 15:20:32
    首先第种是高斯混合模型算法: 高斯模型有单高斯模型(SGM)和混合高斯模型(GMM)两种。 (1)单高斯模型: 为简单起见,阈值t的选取一般靠经验值来设定。通常意义下,我们一般取t=0.7-0.75之间。 ...
  • 高斯混合模型(GMM)

    2019-11-13 11:24:08
    混合模型(Mixture Model) 混合模型是一个可以用来表示在总体分布(distribution)中含有 K 个子分布的概率模型,换句话说,混合模型表示了观测数据在...当样本数据 X 是一维数据(Univariate)时,高斯分布遵从下...
  • 目录高斯模型单高斯模型高斯混合模型 高斯模型 单高斯模型 当样本数据 x∈Rx\in\mathbb{R}x∈R 是一维时,高斯分布服从以下概率密度函数: P(x∣θ)=12πσ2exp(−(x−μ)22σ2)P(x|\theta)=\frac{1}{\sqrt{2\pi\...
  • sklearn之高斯混合模型

    2021-04-26 20:10:49
    什么是高斯分布? 高斯分布也叫正态分布,也就是常态分布,什么意思呢?比如说男性的身高,假如说有...一维正态分布公式:正态分布有两个参数,即期望(均数)μ和标准差σ,σ2为方差。 标准正态分布: ...
  • 机器学习(四)高斯混合模型

    千次阅读 2015-04-24 11:46:37
    开始学习高斯混合模型,需要先简单复习一下单高斯模型的参数估计方法,描述一个高斯模型其实就是要计算它的均值、协方差矩阵(一维空间为方差,二维以上称之为协方差矩阵): 假设有数据集X={x1,x2,x3……,xn},那么...

空空如也

空空如也

1 2 3 4 5 ... 13
收藏数 249
精华内容 99
关键字:

一维高斯混合模型