EM算法

2019-04-02 本文已影响1人 Vophan

含有隐参数的混合模型的参数估计方法

不是所有问题都可以求得解析解（最简单的似然估计可以）。

迭代的算法。

分为E步，M步。
$\theta^{t+1} = \arg\max\int_Z\log P(x,z|\theta)P(z|x,\theta^t)$
E step: 得到Q函数

M step: 求使Q函数最大的 $\theta$

高斯混合模型GMM

高斯：高斯分布

多个高斯分布叠加而成
$p(x) = \sum_{k=1}^K\alpha_kN(\mu_k,\sigma_k),\quad \sum_{k=1}^K\alpha_k=1$
x：观察变量

z：隐变量，表示的对应的样本x属于哪一个高斯分布。（离散的随机变量）

(x,z) :完全数据

为什么用EM算法：

因为在高维的高斯分布问题上，用MLE得不到解析解，首先log里面是连加，很难算出解析解，而且高维的高斯分布式子本身就很难计算。

EM算法的导出：

事实上，EM算法是通过迭代逐步近似极大化 $L(\theta)$ 的，我们希望新的估计值可以让 $L(\theta)$ 增加，即逐步到达最大值。

为此，两者的差：
$L(\theta) - L(\theta_{i}) = \log(\sum_zP(Y|Z,\theta)P(Z|\theta)) - log(P(Y|\theta_i))$
然后利用jenson不等式，得到差的下界：
$L(\theta) - L(\theta_i) = \log(\sum_zP(Y|Z,\theta)\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Y|Z,\theta)}) - log(P(Y|\theta_i))\\ \geq\sum_ZP(Z|Y,\theta_i)\log(\sum_zP(Y|Z,\theta)\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Y|Z,\theta)}) - log(P(Y|\theta_i))\\ =\sum_Z P(Z|Y,\theta_i)\log(\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Y|Z,\theta)P(Y|\theta_i)})$
下界：
$B(\theta, \theta_i) = L(\theta_i) + \sum_Z P(Z|Y,\theta_i)\log(\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Y|Z,\theta)P(Y|\theta_i)})$
为什么要求下界呢？

通过下界，我们可以发现我们将前面提到的log中的连加变成了连乘，而且下界增大的方向一定可以使 $L(\theta)$ 增大。

然后，我们省去在对 $\theta$ 求偏导时为常数的项：
$\theta_{i+1} = \arg \max (\sum_Z P(Z|Y,\theta_i)\log(P(Y,Z|\theta)))\\ =\arg \max_{\theta} Q(\theta, \theta_i)$
由于是近似逼近的方法，所以EM算法并不能保证取到全局最优解。

result

代码：

import numpy as np
import math

class EmAlgorithm:
    """
    Implements of EM algorithm
    """
    def __init__(self, data, models_num):

        self.data = np.array(data)
        self.data_num = len(data)
        self.models_num = models_num
        self.mu = np.array([[1,2]]).T
        self.theta = np.array([[1,1]]).T
        self.alpha = np.array([0.1, 0.9],dtype=float).reshape(1,2)
        self.gama = np.zeros((self.models_num, self.data_num))
        self.max_iters = 500

    def compute_gama(self):
        
        self.gauss = self._compute_Guass(self.data, self.mu, self.theta)
        print(self.gauss,"\n",self.alpha.T)
        self.gama =  self.gauss*self.alpha.T / np.sum(self.gauss*self.alpha.T,axis=0)

    def _compute_Guass(self, x, mu, theta):

        x = np.tile(x, (self.models_num,1))
        mu = np.tile(mu, (1,self.data_num))
        theta = np.tile(theta, (1,self.data_num))
        return 1/np.sqrt(2*np.pi)*theta*np.exp((x-mu)**2/2*theta**2)

    def update(self):

        self.mu = (np.sum(self.gama*self.data, axis=1) / np.sum(self.gama, axis=1)).T.reshape(self.models_num,1)

        self.theta = (np.sum((((np.tile(self.data, (2,1)).T - self.mu.T)**2)*self.gama.T).T,axis=1) / np.sum(self.gama, axis=1)).T.reshape(self.models_num,1)

        self.alpha = ((np.sum(self.gama, axis=1) / self.data_num).T.reshape(self.models_num,1)).T

    def train(self):
        
        for i in range(self.max_iters):
            print(i)
            self.compute_gama()
            self.update()
        print("----------------------------")
        print("alpha:", self.alpha)
        print("mu:", self.mu)
        print("theta:", self.theta)

没时间排版了，有空在做吧。

EM算法

猜你喜欢

热点阅读