实战EM算法与图像分割

EM 算法是求参数极大似然估计的一种方法，它可以从非完整数据集中对参数进行估计，是一种非常简单实用的学习算法。这种方法可以广泛地应用于处理缺损数据、截尾数据以及带有噪声等所谓的不完全数据，可以具体来说，我们可以利用EM算法来填充样本中的缺失数据、发现隐藏变量的值、估计HMM中的参数、估计有限混合分布中的参数以及可以进行无监督聚类等等。

贴相关几个好文章：从最大似然到EM算法浅解

混合高斯模型（Mixtures of Gaussians）和EM算法

斯坦福大学机器学习——EM算法求解高斯混合模型

关于算法的原理推导上面文章介绍的很详细了，在这里再白话一下个人认识的关键几点：

（1）我们要估计的是什么？估计的是未知的参数，谁的参数？你认为你的数据符合的模型的参数，那么就需要认识到我们的数据符合什么模型。看到许多用的都是数据符合高斯（混合）模型，公式也就是正态分布的那个公式：

（2）关于EM算法的初始化参数问题，看了好多文章，感觉都没有怎么说在EM算法里面，参数初始化是多么的重要，直到自己在具体实现EM算法，发现结果怎么怎么不好的时候，查了许多原因才意到这个问题。虽然说EM算法是个鸡蛋问题，但是这个鸡蛋问题怎么感觉有一点特殊，在初始化的时候必须给一组像样的离的不是很远的蛋，最后在反复进化的过程中才能长成完全的鸡，否则的话你在初始化的时候给一个大石头，那么结果一般情况下不太容易进化成鸡。

好了，至于详细原理，网络资源众多就不在累述。然而原理众多但怎么用代码实现的却很少，这里以一个简单的例子（图像的分割，其实就是无监督聚类）来说明实现过程。在聚类时，设置聚类数，那么每一类数据认为是服从一个高斯分布的，那么我们通过EM算法来估计各个类中的高斯参数。关于灰度图像，在聚类里面其实就是一个一维特征的样本，每个像素点视为一个样本点，值得大小视为其具有的特征，这样就可以了，聚类完以后再把它转化为一个图像显示出来，matlab下代码如下：

clc
clear
close all
img = double(imread(‘lena.jpg‘));
%设置聚类数
cluster_num = 4;
%随机初始化参数--至关重要
%我这样选择可以更大程度的贴近真实均值
%从而效果才会更好
mu = (1:cluster_num)./(cluster_num + 1) * max(max(img));
sigma = mu;
pw = zeros(cluster_num,size(img,1)*size(img,2));
pc = rand(1,cluster_num);
pc = pc/sum(pc);%将类概率归一化
max_iter = 50;%以迭代次数来作为停止的条件
iter = 1;
while iter <= max_iter
    %----------E-------------------
    for i = 1:cluster_num
        %矩阵操作--速度快
        MU = repmat(mu(i),size(img,1)*size(img,2),1);
        %高斯模型
        temp = 1/sqrt(2*pi*sigma(i))*exp(-(img(:)-MU).^2/2/sigma(i));
        temp(temp<0.000001) = 0.000001;%防止出现0
%       pw(i,:) = log(pc(i)) + log(temp);
        pw(i,:) = pc(i) * temp;
    end
    pw = pw./(repmat(sum(pw),cluster_num,1));%归一
    %----------M---------------------
    %更新参数集
    for i = 1:cluster_num
         pc(i) = mean(pw(i,:));
         mu(i) = pw(i,:)*img(:)/sum(pw(i,:));
         sigma(i) = pw(i,:)*((img(:)-mu(i)).^2)/sum(pw(i,:));
    end
    %------------show-result---------------
    [~,label] = max(pw);
    %改大小
    label = reshape(label,size(img));
    imshow(label,[])
    title([‘iter = ‘,num2str(iter)]);
    pause(0.1);
    M(iter,:) = mu;
    S(iter,:) = sigma;
    iter = iter + 1;
end
%将均值与方差的迭代过程显示出来
figure
for i = 1:cluster_num
    plot(M(:,i));
    hold on
end
figure
for i = 1:cluster_num
    plot(S(:,i));
    hold on
end

均值随代数变化图

方差随代数变化图

可以看到这是设置4类的变化图，程序里面的参数初始化设置中，尤其是对均值Mu的设置很重要，设置的不好的话分割的效果就很常，也就是说参数估计的也很差，并且分割要求的类数越大（也就是需要估计的参数越多），对初始值的要求越高。