900字范文 > 机器学习十大经典算法：另辟蹊径EM算法+高斯混合模型图像像素分割实战——Nemo鱼图像

机器学习十大经典算法：另辟蹊径EM算法+高斯混合模型图像像素分割实战——Nemo鱼图像

时间：2019-11-19 10:37:15

前言

今天程序员节，据说发博客会有1024勋章，所以就来整理一下笔者模式识别课第一次大作业的第二道题——用EM算法来做Nemo鱼图像像素分割。网上看了很多关于EM算法的分析，大多都涉及很复杂的数学推导，包括要用凸函数的不等关系，要证单调性收敛性等，看了很久也没有理太清楚。在这里分享一种形象化理解EM算法的方法，来源于笔者的模式识别老师。她将EM算法看成是反复进行贝叶斯决策和最大似然估计的过程，当参数变化小于某个阈值时，算法收敛到局部最优值。并附笔者所写代码，供大家参考。贝叶斯理论与实战和最大似然估计在以下两篇博客中已经详细介绍。由于本文严格继承于贝叶斯实战，故往下阅读前，最好先读它。

贝叶斯理论：机器学习十大经典算法：深入浅出聊贝叶斯决策

贝叶斯实战：朴素贝叶斯图像分割实战——Nemo鱼图像分割

代码演示PPT及数据：下载

数据与任务

数据

图像fish.bmp与掩膜mask.mat，掩膜点乘图像，即可获得待分割区域ROI。小鱼区域主要有两种类型的区域（白色，橙色），以下就是要用EM算法+朴素贝叶斯把这两类区域分割出来——用不同的颜色表示不同类的区域。

训练数据sample.mat，它是一个二维的matlab数组，第一列为灰度值，第2-4列为RGB值。由于EM算法属于无监督学习，所以较贝叶斯实战中的训练数据而言，抹去了标签。训练数据中依旧隐含了两类区域的灰度值和RGB值，及每类中它们的分布规律。但哪个值属于哪类区域（隐变量），以及每类区域的灰度值和RGB值的分布是我们要用EM算法去估计的。

任务

任务1：用EM算法估计出每个训练数据的软标签，并结合软标签和训练数据，估计出两类区域灰度值的概率密度函数，最后用最小错误贝叶斯对fish.bmpROI灰度图像逐像素进行分割。

任务2：用EM算法估计出每个训练数据的软标签，并结合软标签和训练数据，估计出两类区域RGB值的概率密度函数，最后用最小错误贝叶斯对fish.bmpROI彩色图像逐像素进行分割。

copyright ©意疏：/sinat_35907936/article/details/109266603

理论与实现

EM算法

E-M即期望最大化，它是一种通过迭代来估计隐变量的利器，常用在无类别标签的高斯混合模型（GMM）中，GMM如下图所示，图源。这里并不打算去强推EM算法公式，而旨在用另一种方法理解它——EM=贝叶斯公式+最大似然估计，以达到在GMM模型中使用它的目的。

所谓高斯混合模型，笔者觉得就是先假定每个类别样本的产生都服从一个高斯分布，然后估计每类对应高斯分布的参数，最后，如果每个类别服从的高斯分布叠加（混合）起来的函数能很好拟合样本数据，说明就找到了一组好的参数。

从贝叶斯实战一文中可知，估计参数可以直接用最大似然估计。但是在训练数据没有标签时，对于GMM模型做分类时，除了要估计均值方差，还要类别估计标签（隐变量），这样使得很难找到似然函数极值点的解析解，此时就需要用EM算法，用迭代的方式解决。

换个角度看EM算法

E过程：估计隐变量，贝叶斯后验概率

对于分类（聚类）而言，隐变量就是指每个样本可能属于的类别标签，但此处的标签不是决策时非0即1的硬标签，而是衡量各类别对产生该样本贡献率的软标签。这由我们的目的决定，EM中不是为了决策，而是为了迭代下去，以产生更好的软标签。

由高斯函数的特点，我们可以知道三个函数在每个点上都是互相叠加的。换句话说，所有取值的产生都与三个分布有关系——三个分布都有产生该取值的可能，或者说都对产生该取值有贡献，只是贡献有大小之分。

通过贝叶斯理论，我们可以知道。在已知结果（采集到的样本是各类别的混合），并且知道结果有几个来源，要得到各来源对结果的贡献率，要用贝叶斯逆概公式。只有两个来源，即二分类时如下式：

由贝叶斯逆概公式，计算出每个类别后验概率（贡献率）的过程就是E过程，贡献率或者后验概率即软标签（soft_guess）。要用上述贝叶斯公式，需要知道先验概率P，类条件概率密度函数的PDF参数m和S。第一次迭代时，它们由我们给定初值，初值给定不是任意的，需要我们观察数据来给出，以保证EM算法尽快收敛。

M过程：更新先验概率和类条件概率密度函数参数，最大似然估计

迭代的目的是为了得到尽可能好分类结果——软标签，而要得到新的软标签，在训练数据给定的情况下，必须要更新先验概率，类条件概率密度函数的均值和方差。由于此时已经有标签，所以可以直接得到先验概率，并可以使用极大似然估计估计出类条件概率密度函数的均值和方差。

由以上分析，每个样本中，属于每个类别的个数为1* 贡献率，或者1* 软标签。那么，每个类别的样本总个数即为其软标签的和，如式（2）。其中N为训练样本个数。

先验概率：

贝叶斯实战中，可以知道最大似然估计均值方差或者均值协方差的解析解，有标签时非此即彼，soft_guess为0或者1，没有标签时。可能为此也可能为彼，soft_guess为某类别在该样本上的后验概率。这里只给出多维时均值协方差的表达式，一维时类似。

均值：累加所有样本值中属于某个类别的成分，再除以该类别的样本总数。

协方差：累加所有样本值中属于某个类别的协方差成分，再除以该类别的样本数减一，减一的目的是无偏估计。

E过程估计软标签，M过程估计先验概率和类条件概率密度函数。EM完成一次，即为一次迭代。如最大似然估计出的参数变化小于给定阈值，即可说EM算法收敛。上述几个式子，与EM算法推出来的是一致的，并且实践证明，它可以收敛。

copyright ©意疏：/sinat_35907936/article/details/109266603

EM实现流程

左边部分为训练过程，即E与M的迭代过程。算法收敛后，用得到的类概率密度函数和先验概率，对Nemo图像ROI区域逐像素进行最小错误贝叶斯决策，为不同的决策分配不同的颜色即可得到分割结果，如图所示。

copyright ©意疏：/sinat_35907936/article/details/109266603

算法实现

数据的加载与解析，并假定一个先验概率——两类各占一半。

import osfrom scipy import iofrom scipy.stats import normimport numpy as npimport PIL.Image as Imageimport matplotlib.pyplot as pltfrom scipy.stats import multivariate_normalplot_dir = 'EM_out'if os.path.exists(plot_dir) == 0:os.mkdir(plot_dir)# 数据加载与解析mask = io.loadmat('mask.mat')['Mask'] # 数据为一个字典，根据key提取数据sample = io.loadmat('sample.mat')['array_sample']src_image = Image.open('fish.bmp')RGB_img = np.array(src_image)Gray_img = np.array(src_image.convert('L'))# 通过mask，获取ROI区域Gray_ROI = (Gray_img * mask)/256RGB_mask = np.array([mask, mask, mask]).transpose(1, 2, 0)RGB_ROI = (RGB_img * RGB_mask)/255# 假设两类数据初始占比相同，即先验概率相同P_pre1 = 0.5P_pre2 = 0.5# 假设每个数据来自两类的初始概率相同，即软标签相同soft_guess1 = 0.5soft_guess2 = 0.5# 选择1维或是多维gray_status = Truegray_status = False

任务1：对灰度图像进行像素分割。经过试验，发现EM算法很容易过拟合，即将数据尽可能的分开，这样使得有些连续的区域分割出来不连续。所以笔者采用每迭代一次就进行一次Nemo与分割的测试，最后选择最好的一组分割结果。分割过程就是最小错误贝叶斯决策的过程。

根据上图中两种颜色，观察出一个初值来。笔者均值假定的0.5,0.8，方差假定的0.1，0.3。

# 一维时的EM# ----------------------------------------------------------------------------------------------------#if gray_status:# 观察图像，肉眼估计初始值gray1_m = 0.5gray1_s = 0.1gray2_m = 0.8gray2_s = 0.3# 绘制假定的PDFx = np.arange(0, 1, 1/1000)gray1_pdf = norm.pdf(x, gray1_m, gray1_s)gray2_pdf = norm.pdf(x, gray2_m, gray2_s)plt.figure(0)ax = plt.subplot(1, 1, 1)ax.plot(x, gray1_pdf, 'r', x, gray2_pdf, 'b')ax.set_title('supposed PDF')plt.figure(1)ax1 = plt.subplot(1, 1, 1)ax1.imshow(Gray_ROI, cmap='gray')ax1.set_title('gray ROI')plt.show()gray = np.zeros((len(sample), 5))gray_s_old = gray1_s + gray2_s# 迭代更新参数for epoch in range(30):for i in range(len(sample)):# 贝叶斯计算每个数据的后验，即得到软标签（E过程）soft_guess1 = (P_pre1*norm.pdf(sample[i][0], gray1_m, gray1_s))/(P_pre1*norm.pdf(sample[i][0], gray1_m, gray1_s) +P_pre2*norm.pdf(sample[i][0], gray2_m, gray2_s))soft_guess2 = 1 - soft_guess1gray[i][0] = sample[i][0]gray[i][1] = soft_guess1*1# 当前一个数据中类别1占的个数，1*后验，显然是小数gray[i][2] = soft_guess2*1gray[i][3] = soft_guess1*sample[i][0] # 对当前数据中属于类别1的部分，当前数据*后验gray[i][4] = soft_guess2*sample[i][0]# 根据软标签，再借助最大似然估计出类条件概率PDF参数——均值，标准差（M过程）gray1_num = sum(gray)[1] # 对每一个数据中类别1占的个数求和，就得到数据中类别1的总数gray2_num = sum(gray)[2]gray1_m = sum(gray)[3]/gray1_num # 对每一个数据中属于类别1的那部分求和，就得到类别1的x的和，用其除以类别1的个数就得到其均值gray2_m = sum(gray)[4]/gray2_numsum_s1 = 0.0sum_s2 = 0.0for i in range(len(gray)):sum_s1 = sum_s1 + gray[i][1]*(gray[i][0] - gray1_m)*(gray[i][0] - gray1_m)# 每个数据的波动中，属于类别1的部分sum_s2 = sum_s2 + gray[i][2]*(gray[i][0] - gray2_m)*(gray[i][0] - gray2_m)gray1_s = pow(sum_s1/gray1_num, 0.5) # 标准差gray2_s = pow(sum_s2/gray2_num, 0.5)# print(gray1_m, gray2_m, gray1_s, gray2_s)P_pre1 = gray1_num/(gray1_num + gray2_num) # 更新先验概率P_pre2 = 1 - P_pre1gray1_pdf = norm.pdf(x, gray1_m, gray1_s)gray2_pdf = norm.pdf(x, gray2_m, gray2_s)gray_s_d = abs(gray_s_old - gray2_s - gray1_s)gray_s_old = gray2_s + gray1_s# if gray_s_d < 0.0001: # 迭代停止条件，如果两次方差变化较小则停止迭代#break# 绘制更新参数后的pdfplt.figure(2)ax2 = plt.subplot(1, 1, 1)ax2.plot(x, gray1_pdf, 'r', x, gray2_pdf, 'b')ax2.set_title('epoch' + str(epoch + 1) + ' PDF')plt.savefig(plot_dir + '//' + 'PDF_' + str(epoch + 1) + '.jpg', dpi=100)plt.close()# plt.show()if epoch % 1 == 0: # 迭代2次进行一次分割测试gray_out = np.zeros_like(Gray_img)for i in range(len(Gray_ROI)):for j in range(len(Gray_ROI[0])):if Gray_ROI[i][j] == 0:continue# 贝叶斯公式分子比较，等价于最大后验elif P_pre1 * norm.pdf(Gray_ROI[i][j], gray1_m, gray1_s) > P_pre2 * norm.pdf(Gray_ROI[i][j],gray2_m, gray2_s):gray_out[i][j] = 100else:gray_out[i][j] = 255# 显示分割结果plt.figure(3)ax3 = plt.subplot(1, 1, 1)ax3.imshow(gray_out, cmap='gray')ax3.set_title('epoch' + str(epoch + 1) + 'gray segment')plt.savefig(plot_dir + '//' + 'Gray_segment_' + str(epoch + 1) + '.jpg', dpi=100)plt.close()# plt.show()

迭代过程概率分布变化情况，从图中可以看出，随着迭代进行，概率分布的交叉部分越来越少，这样会导致过拟合。

分割结果动图展示：

copyright ©意疏：/sinat_35907936/article/details/109266603

任务2：对彩色图像进行像素分割。协方差矩阵为3X3的实对称矩阵，对角线为方差。分割过程是最小错误贝叶斯决策的过程。

# 三维时的EM# -------------------------------------------------------------------------------------------------------#else:# 观察图像，肉眼估计初始值RGB1_m = np.array([0.5, 0.5, 0.5])RGB2_m = np.array([0.8, 0.8, 0.8])RGB1_cov = np.array([[0.1, 0.05, 0.04],[0.05, 0.1, 0.02],[0.04, 0.02, 0.1]])RGB2_cov = np.array([[0.1, 0.05, 0.04],[0.05, 0.1, 0.02],[0.04, 0.02, 0.1]])RGB = np.zeros((len(sample), 11))# 显示彩色ROIplt.figure(3)cx = plt.subplot(1, 1, 1)cx.set_title('RGB ROI')cx.imshow(RGB_ROI)plt.show()# 迭代更新参数for epoch in range(20):for i in range(len(sample)):# 贝叶斯计算每个数据的后验，即得到软标签soft_guess1 = P_pre1*multivariate_normal.pdf(sample[i][1:4], RGB1_m, RGB1_cov)/(P_pre1*multivariate_normal.pdf(sample[i][1:4], RGB1_m, RGB1_cov) + P_pre2*multivariate_normal.pdf(sample[i][1:4], RGB2_m, RGB2_cov))soft_guess2 = 1 - soft_guess1RGB[i][0:3] = sample[i][1:4]RGB[i][3] = soft_guess1*1RGB[i][4] = soft_guess2*1RGB[i][5:8] = soft_guess1*sample[i][1:4]RGB[i][8:11] = soft_guess2*sample[i][1:4]# print(RGB[0])# 根据软标签，再借助最大似然估计出类条件概率PDF参数——均值，标准差RGB1_num = sum(RGB)[3]RGB2_num = sum(RGB)[4]RGB1_m = sum(RGB)[5:8]/RGB1_numRGB2_m = sum(RGB)[8:11]/RGB2_num# print(RGB1_num+RGB2_num, RGB1_m, RGB2_m)cov_sum1 = np.zeros((3, 3))cov_sum2 = np.zeros((3, 3))for i in range(len(RGB)):# print(np.dot((RGB[i][0:3]-RGB1_m).reshape(3, 1), (RGB[i][0:3]-RGB1_m).reshape(1, 3)))cov_sum1 = cov_sum1 + RGB[i][3]*np.dot((RGB[i][0:3]-RGB1_m).reshape(3, 1), (RGB[i][0:3]-RGB1_m).reshape(1, 3))cov_sum2 = cov_sum2 + RGB[i][4]*np.dot((RGB[i][0:3]-RGB2_m).reshape(3, 1), (RGB[i][0:3]-RGB2_m).reshape(1, 3))RGB1_cov = cov_sum1/(RGB1_num-1) # 无偏估计除以N-1RGB2_cov = cov_sum2/(RGB2_num-1)P_pre1 = RGB1_num/(RGB1_num + RGB2_num)P_pre2 = 1 - P_pre1print(RGB1_cov, P_pre1)# 用贝叶斯对彩色图像进行分割RGB_out = np.zeros_like(RGB_ROI)for i in range(len(RGB_ROI)):for j in range(len(RGB_ROI[0])):if np.sum(RGB_ROI[i][j]) == 0:continue# 贝叶斯公式分子比较elif P_pre1 * multivariate_normal.pdf(RGB_ROI[i][j], RGB1_m, RGB1_cov) > P_pre2 * multivariate_normal.pdf(RGB_ROI[i][j], RGB2_m, RGB2_cov):RGB_out[i][j] = [255, 0, 0]else:RGB_out[i][j] = [0, 255, 0]# print(RGB_ROI.shape)# 显示彩色分割结果plt.figure(4)ax3 = plt.subplot(1, 1, 1)ax3.imshow(RGB_out)ax3.set_title('epoch' + str(epoch + 1) + ' RGB segment')plt.savefig(plot_dir + '//' + 'RGB_segment_' + str(epoch + 1) + '.jpg', dpi=100)plt.close()

分割结果动图展示：