高斯混合模型 - 一抹烟霞

使用单高斯模型来建模有一些限制，例如，它一定只有一个众数，它一定对称的。举个例子，如果我们对下面的分布建立单高斯模型，会得到显然相差很多的模型：

于是，我们引入混合高斯模型（Gaussian Mixture Model，GMM）。高斯混合模型就是多个单高斯模型的和。它的表达能力十分强，任何分布都可以用GMM来表示。例如，在下面这个图中，彩色的线表示一个一个的单高斯模型，黑色的线是它们的和，一个高斯混合模型：

在小球检测的栗子中，我们试图对红色小球建立单高斯模型（红和绿这二元），会发现红色小球的观测值不是很符合所建立的模型，如下图：

此时，如果我们采取高斯混合模型（两个二元高斯分布），会发现效果好了很多，如下图：

下面，我们来详细地介绍一下高斯混合模型，高斯混合模型的数学形式如下：

$p(\boldsymbol{x})=\sum_{k=1}^Kw_kg_k(\boldsymbol{x}|\boldsymbol\mu_k,\Sigma_k)$

其中， $g_k$ 是均值为 $\boldsymbol\mu_k$ ，协方差矩阵为 $\Sigma_k$ 的单高斯模型， $w_k$ 是 $g_k$ 的权重系数（ $w_k>0, \sum_{k=1}^Kw_k=1$ ）， $K$ 是单高斯模型的个数。

前面我们提到了GMM的优点（能够表示任何分布），当然GMM也有缺点。其一，参数太多了，每一个单高斯模型都有均值、协方差矩阵和权重系数，另外还有单高斯模型的个数 $K$ 也是其中一个参数，这使得求解GMM十分复杂。其二，有时候所建立的高斯混合模型对观测值的建模效果很好，但是其他值可能效果不好，我们称这个缺点为过度拟合（Overfitting）。

2、求解高斯混合模型：EM算法

在求解单高斯分布的时候，我们用最大似然估计（MLE）的方法得到了理论上的最优解。当我们使用相同的方法试图求解高斯混合模型的时候，会卡在中间步骤上（具体来说，是单高斯分布求和出现在了对数函数里面）。索性我们可以用迭代的方法来求解GMM，具体来说，最大期望算法（Expectation Maximization algorith，EM）。

上一节提到，我们想要求解GMM如下：

$p(\boldsymbol{x})=\sum_{k=1}^Kw_kg_k(\boldsymbol{x}|\boldsymbol\mu_k,\Sigma_k)$

这其中需要求解的变量很多： $K$ 、 $w_k$ 、 $\boldsymbol\mu_k$ 、 $\Sigma_k$ 。为了简化问题，我们假定 $K$ 是预先设定好的，并且每个单高斯分布的权重相等，即假定：

$w_k = \frac{1}{K},k=1,2,\cdots,K$

这样，我们需要确定的就只剩下每个单高斯分布的参数（ $\boldsymbol\mu_k$ 、 $\Sigma_k$ ）了。

前面提到，这个模型是没有解析解（理论最优）的。取而代之，我们采取迭代的方法逼近最优解（大家可以回想一下牛顿法求近似求解方程，或者遗传算法）。

在牛顿法中，我们需要预先猜测一个初始解，同样在EM算法中，我们也需要预先猜测一个初始解（ $\boldsymbol\mu_1,\Sigma_1,\boldsymbol\mu_2,\Sigma_2,\cdots,\boldsymbol\mu_K,\Sigma_K$ ）。

在EM算法中，我们引入一个中间变量 $z_k^i = \frac{g_k(\boldsymbol{x_i}|\boldsymbol\mu_k,\Sigma_k)}{\sum_{i=1}^Kg_k(\boldsymbol{x_i}|\boldsymbol\mu_k,\Sigma_k)}(k=1,2,\cdots,K;\ i=1,2,\cdots,N)$ ，其中 $K$ 是单高斯分布的个数， $N$ 是观测值的数目。

直观地可以这样理解：在确定了所有单高斯分布之后，可以计算观测值 $\boldsymbol{x_i}$ 发生的概率（分母部分）和第 $k$ 个单高斯分布 $g_k(\boldsymbol\mu_k,\Sigma_k)$ 下观测值 $x_i$ 发生的概率（分子部分），这样 $z^i_k$ 就可以理解为第 $k$ 个高斯分布对观测值 $x_i$ 发生的概率的贡献了。如下表所示：

$\begin{tabular}{|c|c|c|c|c|c|c|c|} \hline & $\frac{1}{K}$ & $\frac{1}{K}$ & $\cdots$ & $\frac{1}{K}$ & $\cdots$ & $\frac{1}{K}$ \\ & $\boldsymbol\mu_1$ & $\boldsymbol\mu_2$ & $\cdots$ & $\boldsymbol\mu_k$ & $\cdots$ & $\boldsymbol\mu_K$ \\ & $\Sigma_1$ & $\Sigma_2$ & $\cdots$ & $\Sigma_k$ & $\cdots$ & $\Sigma_K$ \\ \hline $x_1$ & $z^1_1$ & $z^1_2$ & $\cdots$ & $z^1_k$ & $\cdots$ & $z^1_K$ \\ \hline $x_2$ & $z^2_1$ & $z^2_2$ & $\cdots$ & $z^2_k$ & $\cdots$ & $z^2_K$ \\ \hline $\vdots$ & $\vdots$ & $\vdots$ & & $\vdots$ & & $\vdots$ \\ \hline $x_i$ & $z^i_1$ & $z^i_2$ & $\cdots$ & $z^i_k$ & $\cdots$ & $z^i_K$ \\ \hline $\vdots$ & $\vdots$ & $\vdots$ & & $\vdots$ & & $\vdots$ \\ \hline $x_N$ & $z^N_1$ & $z^N_2$ & $\cdots$ & $z^N_k$ & $\cdots$ & $z^N_K$ \\ \hline & $z_1$ & $z_2$ & $\cdots$ & $z_k$ & $\cdots$ & $z_K$ \\ \hline \end{tabular}$

表中最后一行 $z_k=\sum_{i=1}^Nz^i_k$ 可以理解为第 $k$ 个单高斯分布对整个GMM的贡献。

例如，当GMM中 $K=2$ 时（即由两个单高斯分布组成）：

在这个图中，对于观测值 $\boldsymbol{x_i}$ ， $p_1/2$ 是第一个单高斯分布 $g_1$ 下 $\boldsymbol{x_i}$ 发生的概率， $p_2/2$ 是第二个单高斯分布 $g_2$ 下 $\boldsymbol{x_i}$ 发生的概率， $(p_1+p_2)/2$ 是在整个GMM下 $\boldsymbol{x_i}$ 发生的概率； $z_1^i=\frac{p_1}{p_1+p_2}$ 表示 $g_1$ 对 $\boldsymbol{x_i}$ 发生概率的贡献， $z_2^i=\frac{p_2}{p_1+p_2}$ 表示 $g_2$ 对 $\boldsymbol{x_i}$ 发生概率的贡献。

当我们算出了所有的 $z_k^i$ 之后，我们再反过来用它更新每一个单高斯分布。更新的公式为： $\boldsymbol\mu_k = \frac{1}{z_k}\sum_{i=1}^Nz_k^i\boldsymbol{x_i} <p>" /><br /><img src=$

大家可以想一想，对于第 $k$ 个高斯分布 $g_k$ ， $z_k^i/z_k$ 可以理解为第 $i$ 个观测值 $\boldsymbol{x_i}$ 的贡献，而 $\boldsymbol\mu_k$ 表示 $g_k$ 的平均值，用 $\frac{1}{z_k}\sum_{i=1}^Nz_k^i\boldsymbol{x_i}$ 来更新 $\boldsymbol\mu_k$ 就很好理解了，进而再更新协方差矩阵 $\Sigma_k$ 。

这样，两组值互相更新，当相邻两个步骤值的差小于事先设定一个阈值（threshold）时（收敛），我们停止循环，输出 $\boldsymbol\mu_k$ 和 $\Sigma_k$ 为所求。

值得一提的是，EM算法的结果和初值有很大关系，不同初值会得到不同的解。（想象一下GMM模型其实是多峰的，不同的初值可能最终收敛到不同的峰值，有的初值会收敛到全局最优，有的初值会收敛到局部最优）

实战：小球检测

具体任务在这个网址：小球检测

简单来说，给你一些图片作为训练集，让你对黄色小球进行建模，建模完成后，读入最左边的原始图片，要求能像最右边的图片一样标识出小球的位置：

（注：我们使用matlab完成所有代码）

第一步：建立颜色模型

前面几节我们提到，要建立高斯分布模型，首先要有很多观测值，所以，我们将在这一步里面采集很多观测值，再建立高斯模型。课程里面提供了19张图片作为训练集供采样：

我们编写代码mark.m，作用为读入一张一张图片，然后手动圈出图片中的黄色小球，将圈中部分的像素RGB值存入观测值。

imagepath = './train';
Samples = [];
for k=1:19
    I = imread(sprintf('%s/%03d.png',imagepath,k)); % read files into RGB
    R = I(:,:,1);
    G = I(:,:,2);
    B = I(:,:,3);
</span>%<span style="color: #000000;"> Collect samples 
disp(</span><span style="color: #800000;">''</span><span style="color: #000000;">);
disp(</span><span style="color: #800000;">'</span><span style="color: #800000;">INTRUCTION: Click along the boundary of the ball. Double-click when you get back to the initial point.</span><span style="color: #800000;">'</span><span style="color: #000000;">)
disp(</span><span style="color: #800000;">'</span><span style="color: #800000;">INTRUCTION: You can maximize the window size of the figure for precise clicks.</span><span style="color: #800000;">'</span><span style="color: #000000;">)
figure(</span><span style="color: #800080;">1</span>), mask =<span style="color: #000000;"> roipoly(I); 
figure(</span><span style="color: #800080;">2</span>), imshow(mask); title(<span style="color: #800000;">'</span><span style="color: #800000;">Mask</span><span style="color: #800000;">'</span><span style="color: #000000;">);

sample_ind </span>= find(mask &gt; <span style="color: #800080;">0</span>); % <span style="color: #0000ff;">select</span><span style="color: #000000;"> marked pixels
R </span>=<span style="color: #000000;"> R(sample_ind);
G </span>=<span style="color: #000000;"> G(sample_ind);
B </span>=<span style="color: #000000;"> B(sample_ind);

Samples </span>= [Samples; [R G B]]; %<span style="color: #000000;"> insert selected pixels into samples

disp(</span><span style="color: #800000;">'</span><span style="color: #800000;">INTRUCTION: Press any key to continue. (Ctrl+c to exit)</span><span style="color: #800000;">'</span><span style="color: #000000;">)
pause

end
save('Samples.mat', 'Samples'); % save the samples to file
figure,

scatter3(Samples(:,1),Samples(:,2),Samples(:,3),'.');

title('Pixel Color Distribubtion');

xlabel('Red');

ylabel('Green');

zlabel('Blue');

结果我们采集了6699个观测值，散点图如下：
我们选用多元高斯分布作为我们的模型，应用第四节中的结论，编写代码coefficient.m求出平均值 $\boldsymbol\mu$ 和协方差矩阵 $\Sigma$ 并存储到本地以供后续使用。

mu = mean(Samples); % mu
sig=zeros(3,3); % sigma

for i=1:N

data=double(Samples(i,:));

sig=sig+(data-mu)'*(data-mu)/N;

end
% save the coefficients to files

save('mu.mat', 'mu');

save('sig.mat', 'sig');

接下来，我们需要编写函数detecBall，以图像为参数，以划分好的黑白图像和小球中心为输出，如下，I是输入的图像，segI是划分好的黑白图像，loc是球的中心点。detecBall.m的具体代码如下：

function [segI, loc] = detectBall(I)
load('mu.mat');
load('sig.mat');
Id=double(I); % array in size of (row, col, 3)
row=size(Id,1);

col=size(Id,2);
% x_i - mu

for i=1:3

Id(:,:,i) = Id(:,:,i) - mu(i);

end
% reshape the image to a matrix in size of (rowcol, 3)

Id=reshape(Id,rowcol,3);
% calc possibility using gaussian distribution

% be careful of using * and .* in matrix multiply

Id = exp(-0.5* sum(Idinv(sig).Id, 2)) ./ (2*pi)^1.5 ./ det(sig)^0.5;
% reshape back, now each pixels is with the value of the possibility

Id=reshape(Id,row,col);
% set threshold

thr=8e-06;
% binary image about if each pixel 'is ball'

Id=Id>thr;
% find the biggest ball area

segI = false(size(Id));
CC = bwconncomp(Id);

numPixels = cellfun(@numel,CC.PixelIdxList);

[biggest,idx] = max(numPixels);

segI(CC.PixelIdxList{idx}) = true;

%figure, imshow(segI); hold on;
S = regionprops(CC,'Centroid');

loc = S(idx).Centroid;

%plot(loc(1), loc(2),'r+');

需要注意的一点是，求像素属于小球的概率的时候，尽量用矩阵运算，而不要用循环，否则会效率低下，请读者仔细揣摩计算技巧，灵活运用矩阵的 $.*$ 运算。

接下来，我们就可以测试小球检测的效果啦！

我们编写test.m测试训练集的19张图片：

imagepath = './train';
for k=1:19
    I = imread(sprintf('%s/%03d.png',imagepath,k));
    [segI, loc] = detectBall(I);
    figure, imshow(segI); hold on; 
    plot(loc(1), loc(2), '+b','MarkerSize',7); 
    disp('Press any key to continue. (Ctrl+c to exit)')
    pause
end

效果如下：

posted on 2019-08-23 21:59 一抹烟霞阅读(8729) 评论(1) 编辑收藏举报

刷新页面返回顶部

2、求解高斯混合模型：EM算法

在这个图中，对于观测值，是第一个单高斯分布下发生的概率，是第二个单高斯分布下发生的概率，是在整个GMM下发生的概率；表示对发生概率的贡献，表示对发生概率的贡献。

当我们算出了所有的之后，我们再反过来用它更新每一个单高斯分布。更新的公式为：

这样，两组值互相更新，当相邻两个步骤值的差小于事先设定一个阈值（threshold）时（收敛），我们停止循环，输出和为所求。

实战：小球检测

当我们算出了所有的 $z_k^i$ 之后，我们再反过来用它更新每一个单高斯分布。更新的公式为： $\boldsymbol\mu_k = \frac{1}{z_k}\sum_{i=1}^Nz_k^i\boldsymbol{x_i} <p>" /><br /><img src=$

这样，两组值互相更新，当相邻两个步骤值的差小于事先设定一个阈值（threshold）时（收敛），我们停止循环，输出 $\boldsymbol\mu_k$ 和 $\Sigma_k$ 为所求。