聊聊EM算法~

      EM算法(Expection Maximuzation)的中文名称叫做期望最大化，我的理解EM算法就是一种引入隐含变量的坐标向上法，它与其他优化方法的目的相同，就是求解一个数学模型的最优参数，不同之处在于EM算法交互迭代的对象是隐含变量与模型参数，一般隐含变量表现为数据的类别。期望说白了就是数据的权重平均，在EM算法中，可以理解为数据的类别，那么既然确定好了数据的类别，下一步就是想办法求取新数据构成分布的参数，一般采用最大似然法求取最优参数。剩下的就是最普通的迭代过程，先初始化参数，计算数据的概率分布，再对数据进行重新分类，然后重新计算参数，周而复始。幸运的是，EM算法的收敛性能够得到保证。

      曾经读过一篇特别好的博客，通过举例的方式对EM算法解释的很形象，现摘抄如下：

     “假设我们需要调查我们学校的男生和女生的身高分布。你怎么做啊？你说那么多人不可能一个一个去问吧，肯定是抽样了。假设你在校园里随便地活捉了100个男生和100个女生。他们共200个人（也就是200个身高的样本数据，为了方便表示，下面，我说“人”的意思就是对应的身高）都在教室里面了。那下一步怎么办啊？你开始喊：“男的左边，女的右边，其他的站中间！”。然后你就先统计抽样得到的100个男生的身高。假设他们的身高是服从高斯分布的。但是这个分布的均值u和方差∂²我们不知道，这两个参数就是我们要估计的。记作θ=[u, ∂]^T。

     再回到例子本身，如果没有“男的左边，女的右边，其他的站中间！”这个步骤，或者说我抽到这200个人中，某些男生和某些女生一见钟情，已经好上了，纠缠起来了。咱们也不想那么残忍，硬把他们拉扯开。那现在这200个人已经混到一起了，这时候，你从这200个人（的身高）里面随便给我指一个人（的身高），我都无法确定这个人（的身高）是男生（的身高）还是女生（的身高）。也就是说你不知道抽取的那200个人里面的每一个人到底是从男生的那个身高分布里面抽取的，还是女生的那个身高分布抽取的。用数学的语言就是，抽取得到的每个样本都不知道是从哪个分布抽取的。

        这个时候，对于每一个样本或者你抽取到的人，就有两个东西需要猜测或者估计的了，一是这个人是男的还是女的？二是男生和女生对应的身高的高斯分布的参数是多少？

        EM的意思是“Expectation Maximization”，在我们上面这个问题里面，我们是先随便猜一下男生（身高）的正态分布的参数：如均值和方差是多少。例如男生的均值是1米7，方差是0.1米（当然了，刚开始肯定没那么准），然后计算出每个人更可能属于第一个还是第二个正态分布中的（例如，这个人的身高是1米8，那很明显，他最大可能属于男生的那个分布），这个是属于Expectation一步。有了每个人的归属，或者说我们已经大概地按上面的方法将这200个人分为男生和女生两部分，我们就可以根据之前说的最大似然那样，通过这些被大概分为男生的n个人来重新估计第一个分布的参数，女生的那个分布同样方法重新估计。这个是Maximization。然后，当我们更新了这两个分布的时候，每一个属于这两个分布的概率又变了，那么我们就再需要调整E步……如此往复，直到参数基本不再发生变化为止。”

 这是我读过最好的EM算法例子，我对这个例子进行了仿真，代码如下所示：

% 函数功能：根据身高数据计算男性女性平均身高
% 1.对初始值很敏感，如果先验初始值不准确，最后必定收敛到不准确值
% 2.对初始数据分布比较敏感
% 3.初始值不能相同，否则将会陷入无限循环
function EM_algorithm()
    clear all;
    close all;
    clc;
    
    person_num = 600;

    % 男生数据
    u1 = 1.76;
    sigma1 = 0.1;
    x1 = normrnd(u1, sigma1, person_num / 2, 1); 
    y1 = gauss_function(u1, sigma1, x1);
%     plot(x1, ones(person_num / 2, 1), 'r*');
    plot(x1,y1, 'r*');
    
    % 女生数据
    u2 = 1.55;
    sigma2 = 0.1;
    x2 = normrnd(u2, sigma2, person_num / 2, 1); 
    y2 = gauss_function(u2, sigma2, x2);
    hold on;
%     plot(x2, 2 * ones(person_num / 2, 1), 'bo');
    plot(x2, y2, 'bo');
    keyboard
    
    % EM算法
    flag        = 1;
    iter        = 1;
    sigma_joint = 0.1;
    u0_man      = 1.5;
    u0_woman    = 1.4;
    persons     = [x1; x2];
    theold      = 0.001;
    while(flag == 1)
        % 确定样本类别，E步骤
        z = [];
        for i = 1 : length(persons)
            p_man = gauss_function(u0_man, sigma_joint, persons(i));
            p_woman = gauss_function(u0_woman, sigma_joint, persons(i));
            if p_man > p_woman
                z = [z, 1];
            else
                z = [z, -1];
            end
        end
        
        % 提取样本类别序列
        man_index = find(z == 1);
        woman_index = find(z == -1);
        
        % 最大似然估计，M步骤
%         u0_man = fminsearch(@(k)StdMonochrome(k,G),[1.5, 2]);
%         u0_woman = fminsearch(@(k)StdMonochrome(k,G),[1.5, 2]);

        % 记录之前的参数值
        u0_man_old = u0_man;
        u0_woman_old = u0_woman;
        
        disp(['u0_man = ', num2str(u0_man_old)]);
        disp(['u0_woman = ',num2str(u0_woman_old)]);
        disp(['number of man = ', num2str(length(man_index))]);
        disp(['number of woman = ',num2str(length(woman_index))]);

        % 针对于高斯分布这种情况，最大似然估计解的解析式就是均值
        u0_man = mean(persons(man_index), 1);
        u0_woman = mean(persons(woman_index), 1);
        
        if abs(u0_man - u0_man_old) < theold && abs(u0_woman_old - u0_woman) < theold
            flag = 0;
            disp('success');
            disp(['u0_man = ', num2str(u0_man)]);
            disp(['u0_woman = ',num2str(u0_woman)]);
            disp(['iter = ', num2str(iter)]);
        end
        
        iter = iter + 1;
    end
end

% 高斯核函数，经常写
function[p] = gauss_function(u, sigma, data)
%     p = (1/sqrt(2*pi)).*exp(-(data - u).* (data - u)/2 * sigma^2);
    p = exp(-(data - u).* (data - u)/sigma^2);
end

其中，数据分布的具体参数有两个，代码中我直接用的是参数经过最大似然估计后的解析解。当然，也可以直接调用Matlab中的函数进行最大似然估计，不过这个例子就没太大必要了。

抽样数据的真实分布如下所示：

运行代码，得到的迭代数据如下所示，算法迭代9步最终收敛：

u0_man = 1.5
u0_woman = 1.4
number of man = 557
number of woman = 43
u0_man = 1.6714
u0_woman = 1.4049
number of man = 464
number of woman = 136
u0_man = 1.7067
u0_woman = 1.4669
number of man = 400
number of woman = 200
u0_man = 1.7295
u0_woman = 1.498
number of man = 355
number of woman = 245
u0_man = 1.7458
u0_woman = 1.5169
number of man = 326
number of woman = 274
u0_man = 1.7567
u0_woman = 1.5282
number of man = 307
number of woman = 293
u0_man = 1.7641
u0_woman = 1.5352
number of man = 296
number of woman = 304
u0_man = 1.7685
u0_woman = 1.5392
number of man = 292
number of woman = 308
u0_man = 1.7701
u0_woman = 1.5407
number of man = 291
number of woman = 309
success
u0_man = 1.7705
u0_woman = 1.541
iter = 9

其结果与真实值（男：1.76米，女：1.55米）较为接近。但是，经过几个实验得知，EM算法对初始值的敏感程度较高，同时，假设有两类数据，那么这两类数据的类间距离应该满足一定条件，如果我们将分布调整为（男：1.56米，女：1.55米），那么EM算法将无法收敛。大家可以试一试看。

这段时间一直比较忙，忙于写代码实现文献&改进算法，好久没写博客了，如有问题，请代价不吝指教！

参考文献：http://blog.csdn.net/zouxy09