GMM

  许多概率模型有一系列可见变量$v$和一系列潜变量$h$，这时常常会涉及推断困难，就是指难以计算$p(h|v)$或其期望，而这样的操作在一些诸如最大似然学习的任务中往往是必需的。为此可以把精确推断问题描述为一个优化问题，借此推导出推断算法。为了构造这样一个优化问题，假设一个具有可见变量$v$和潜变量$h$的概率模型，按照最大似然估计，我们希望计算观察数据的对数概率$logp(v;\theta )$，则有
$$logp(v;\theta )=log\sum _{h}p(x,h;\theta )\text{\hspace{1em}(1)}$$
  有时候，边缘化消去$h$的操作很费时或难以计算，作为替代，我们可以计算一个$logp(v;\theta )$的证据下界$\mathcal{L}(v,\theta ,q)$
$$\begin{gathered} logp(v;\theta )=log\sum _{h}p(x,h;\theta )=log\sum _{h}q(h|v)\frac{p(v,h;\theta )}{q(h|v)} \\ \ge \sum _{h}q(h|v)log\frac{p(v,h;\theta )}{q(h|v)}={\mathbb{E}}_{h\sim q}[logp(h,v;\theta )-logq(h|v)] \\ =logp(v;\theta )-{\mathbb{E}}_{h\sim q}[logq(h|v)-logp(h,v;\theta )+logp(v;\theta )] \\ =logp(v;\theta )-{\mathbb{E}}_{h\sim q}log\frac{q(h|v)}{\frac{p(v,h;\theta )}{p(v;\theta )}}=logp(v;\theta )-{\mathbb{E}}_{h\sim q}log\frac{q(h|v)}{p(h|v)} \\ =logp(v;\theta )-D_{KL}(q(h|v)||p(h|v;\theta ))={\mathbb{E}}_{h\sim q}[logp(v,h)]+H(q)=\mathcal{L}(v,\theta ,q)\text{\hspace{1em}(2)} \end{gathered}$$
  其中$q(h|v)$表示给定一个数据的可见变量$v$，其潜变量为$h$的概率，对于一个选择的合适分布$q$来说，$\mathcal{L}$是容易计算的，对任意分布$q$的选择来说，$\mathcal{L}$提供了似然函数的一个下界。越好地近似$p(h|v)$的分布$q(h|v)$，其下界越紧，当$q(h|v)=p(h|v)$时，这个近似是完美的，也意味着$\mathcal{L}(v,\theta ,q)=logp(v;\theta )$，因此可以将推断问题看作找一个分布$q$使$\mathcal{L}$最大的过程。
  在潜变量模型中，EM算法是非常常见的训练方法，这是一种能够学到近似后验的算法。

• E步：令${\theta }^{(0)}$表示在这一步开始时的参数值。对任何我们想要训练的（对所有的或者小批量数据均成立）索引为$i$的训练样本$v^{(i)}$，令$q(h|v^{(i)})=p(h|v^{(i)};\theta )$。通过这个定义，我们认为$q$在当前参数${\theta }^{(0)}$下定义，如果我们改变$\theta$，那么$p(h|v;\theta )$将会相应变化，但是$q(h|v)$还是不变并且等于$p(h|v;{\theta }^{(0)})$。
• M步：使用选择的优化算法完全地或部分地关于$\theta$最大化${\sum }_i\mathcal{L}(v^{(i)},\theta ,q)={\sum }_i{\sum }_{h}q(h|v^{(i)})log\frac{p(v^{(i)},h;\theta )}{q(h|v^{(i)})}$。
    对于高斯混合模型，其概率分布为
  $$p(v;\theta )=\sum _{k=1}^K{\pi }_{k}N(v|{\mu }_k,{\Sigma }_k)\text{\hspace{1em}(3)}$$
  其中${\sum }_{k=1}^K{\pi }_k=1$，这时的$v$是数据样本，潜变量 h = { h 1 , h 2 , ⋯   , h K } , h k = { 0 , 1 } h=\{h_1,h_2,\cdots,h_K\},h_k=\{0,1\} h={h1​,h2​,⋯,hK​},hk​={0,1}代表数据样本是否来自$K$个高斯分布中的第$k$个高斯分布，在E步中有
  $$q(h_k=1|v^{(i)})=\frac{{\pi }_{k}N(v^{(i)}|{\mu }_k,{\Sigma }_k)}{{\sum }_{k=1}^K{\pi }_{k}N(v^{(i)}|{\mu }_k,{\Sigma }_k)}=\frac{{\pi }_{k}N(v^{(i)}|{\mu }_k,{\Sigma }_k)}{p(v^{(i)};\theta )}\text{\hspace{1em}(4)}$$
    在M步中，为了求解模型参数${\pi }_k,{\mu }_k,{\Sigma }_k$，需要对各参数求偏导
  $$J_{\pi }=\sum _i\sum _{k}q(h_k|v^{(i)})log({\pi }_k)+\lambda (\sum _k{\pi }_k-1)\text{\hspace{1em}(5)}$$
  $$\frac{\partial J_{\pi }}{\partial {\pi }_k}=\sum _{i}q(h_k|v^{(i)})\frac{1}{{\pi }_k}+\lambda =0\text{\hspace{1em}(6)}$$
  设样本数为$m$，因为${\sum }_k{\pi }_k=1$，所以有
  $$\sum _i\sum _{k}q(h_k|v^{(i)})=-\lambda \sum _k{\pi }_k=-\lambda =\sum _{i}1=m\text{\hspace{1em}(7)}$$
  $${\pi }_k=\frac{{\sum }_{i}q(h_k|v^{(i)})}{m}\text{\hspace{1em}(8)}$$
  同理求导则有
  $${\mu }_k=\frac{{\sum }_{i}q(h_k|v^{(i)})v^{(i)}}{{\sum }_{i}q(h_k|v^{(i)})}\text{\hspace{1em}(9)}$$
  $${\Sigma }_k=\frac{{\sum }_{i}q(h_k|v^{(i)})(v^{(i)}-{\mu }_k)(v^{(i)}-{\mu }_k{)}^T}{{\sum }_{i}q(h_k|v^{(i)})}\text{\hspace{1em}(10)}$$

Matlab代码实现

  在代码实践中可能会出现一个问题，GMM迭代到某一次循环时，某类分布的方差会突然变的很小，导致一些本该分到其分布上的样本会因为概率太小而分到其它类，由于某类分布的样本较少甚至为零，计算过程中有一个除以每类样本数的步骤，这时可能会导致除以零的发生，从而产生非数$NaN$，程序会运行出错。为了缓解这个问题，在第一步时运用$kmeans$算法，对各类分布的均值方差参数进行初始化，使各类分布的均值均匀散布到数据样本的区域范围。

classdef GMM
    properties
        numOfGMs;                         %高斯分布的个数
        pDim;                             %可见变量的空间维数
        pMix;                             %每个分布的混合比例
        mu;                               %均值参数
        sigma;                            %方差参数
        trainData_x;                      %训练数据
        trainData_q_h;                    %训练数据的潜变量编码分布
        N_sample;                         %训练数据样本总数
        batchSize;                        %批次大小
        IterMax;                          %最大迭代次数
        iter_train;                       %当前迭代次数
        isInited;                         %是否已经执行初始化
    end
    methods
        function obj = GMM(pDim,numOfGMs) %类库构造函数
            obj.numOfGMs = numOfGMs;      %初始化类对象参数
            obj.pDim = pDim;
            obj.pMix = 1 / numOfGMs * ones(1,numOfGMs);
            obj.mu = zeros(numOfGMs,pDim) + rand(numOfGMs,pDim) * 0.1;
            obj.sigma = zeros(pDim,pDim,numOfGMs);
            for k = 1 : obj.numOfGMs
                obj.sigma(:,:,k) = eye(obj.pDim);
            end
            obj.batchSize = 100;
            obj.iter_train = 0;
            obj.IterMax = 100;
            obj.isInited = 0;
        end
        function obj = train(obj,trainData_x,Iter)%训练
            if(nargin == 2)
                Iter = obj.IterMax;
            end
            if obj.isInited == 0
                obj = obj.trainInit(trainData_x);
            end
            obj.N_sample = size(obj.trainData_x,1);
            batch_size = obj.batchSize;
            trainData_q_h_batch = zeros(batch_size,obj.numOfGMs);
            diagMat = eye(batch_size);
            diagMat0 = find(diagMat);
            findDiag = diagMat0;
            while obj.iter_train < Iter
                obj.iter_train = obj.iter_train + 1;
                trainData_x_batch = obj.trainData_x(randperm(obj.N_sample , batch_size) , :);
                %E步，求解q(h|v)
                for k = 1 : obj.numOfGMs
                    trainData_q_h_batch(:,k) = obj.pMix(k) * mvnpdf(trainData_x_batch,obj.mu(k,:),obj.sigma(:,:,k));
                end
                trainData_q_h_batch = trainData_q_h_batch ./ sum(trainData_q_h_batch,2);
                %M步，求解模型参数
                obj.pMix = mean(trainData_q_h_batch);
                sumTrainDataQH = sum(trainData_q_h_batch)';
                obj.mu = trainData_q_h_batch' * trainData_x_batch ./ sumTrainDataQH;
                for k = 1 : obj.numOfGMs
                    diagMat(findDiag) = trainData_q_h_batch(:,k);
                    obj.sigma(:,:,k) = (trainData_x_batch - obj.mu(k,:))' * diagMat * (trainData_x_batch - obj.mu(k,:)) / sumTrainDataQH(k);
                end
            end
        end
        
        function [obj] = predict(obj)%预测重构
            batch_size = obj.N_sample;
            indexRandBatch = randperm(obj.N_sample , batch_size);
            testData_x_batch = obj.trainData_x(indexRandBatch , :);
            %求解q(h|v)
            trainData_q_h_batch = zeros(batch_size,obj.numOfGMs);
            for k = 1 : obj.numOfGMs
                trainData_q_h_batch(:,k) = obj.pMix(k) * mvnpdf(testData_x_batch,obj.mu(k,:),obj.sigma(:,:,k));
            end
            trainData_q_h_batch = trainData_q_h_batch ./ sum(trainData_q_h_batch,2);
            obj.trainData_q_h(indexRandBatch,:) = trainData_q_h_batch;
        end
        
        function [obj]=trainInit(obj,trainData_x)%使用kmeans进行初始化
            obj.isInited = 1;
            obj.trainData_x = trainData_x;
            labelMat = k_means(obj.trainData_x,obj.numOfGMs,100);
            labelCenter = obj.trainData_x' * labelMat ./ sum(labelMat);
            obj.mu = labelCenter';
            for k = 1 : obj.numOfGMs
                obj.sigma(:,:,k) = cov(obj.trainData_x(labelMat(:,k)==1,:));
            end
        end
    end
end

function yy = k_means(x0 , k , iter)
% kmeans修改输出参数y为yy，进行GMM的初始化
[n , p] = size(x0);
x = repmat(x0 , 1 , k);
%% 随机分配初始类中心
c = repmat( reshape( x( randperm(n , k) , 1:p)' , 1 , k * p) , n , 1);
%% 求和矩阵
sum_mat = kron(eye(k) , ones(p , 1));
%% 迭代计算
for i = 1 : iter
    distance = ( x - c ).^2 * sum_mat;
    yy = ( distance==min(distance , [] , 2) );
    c0 = x' * yy ./ repmat(sum(yy) , k * p , 1);
    c = repmat( reshape(c0(1 : p , :) , 1 , k * p) , n , 1);
end
end