关闭

聊聊EM算法~

1506人阅读 评论(0) 收藏 举报
分类:

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

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


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

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

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

        EM的意思是“Expectation Maximization”,在我们上面这个问题里面,我们是先随便猜一下男生(身高)的正态分布的参数:如均值和方差是多少。例如男生的均值是17,方差是0.1米(当然了,刚开始肯定没那么准),然后计算出每个人更可能属于第一个还是第二个正态分布中的(例如,这个人的身高是18,那很明显,他最大可能属于男生的那个分布),这个是属于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

2
1

查看评论
* 以上用户言论只代表其个人观点,不代表CSDN网站的观点或立场
    个人资料
    • 访问:235225次
    • 积分:3164
    • 等级:
    • 排名:第11507名
    • 原创:42篇
    • 转载:7篇
    • 译文:0篇
    • 评论:469条
    博客专栏
    最新评论