机器学习笔记|机器学习笔记 - 20. EM算法实践（讲师（邹博））机器学习笔记-20.EM算法实践（

主要内容
文章图片
2018-12-12 20_25_51-【邹博_chinahadoop】机器学习升级版VII（七）.png 多维高斯混合分布聚类 EM算法的聚类效果或许比K均值聚类好一些。

文章图片
2018-12-12 20_27_30-【邹博_chinahadoop】机器学习升级版VII（七）.png 如图，对于二维数据形成概率密度曲线，或者说等值线：

文章图片
2018-12-12 20_28_09-【邹博_chinahadoop】机器学习升级版VII（七）.png 这个图也说明，身高一定符合高斯分布，不一定对。
下图表明，男性符合几个混合高斯分布，女性符合几个混合高斯分布

文章图片
2018-12-25 19_27_33-【邹博_chinahadoop】机器学习升级版VII（七）.png 问答

问：归一化的几种优劣之处？
答：比如做min-max，Scalar或标准版，如果数据服从均匀分布，可能做min-max好一些，但是如果数据服从高斯分布，可能标准化更好一些。
问：为什么鸢尾花没有隐变量？鸢尾花也可能有某个未知的特征决定它的分类，任何分布不是都可能有隐变量么？
答：是的。也许花萼长宽与花瓣长宽并不是鸢尾花最重要特征。只是没有提取
问：针对Kmeans多特征的情景，是不是用PCA处理以后，变成3维或2维，然后再用聚类的方式处理？
答：如果不需要做算法解释的话，这么做是合理的；但是如果需要做算法解释，建议不要用PCA，否则特征无法解释。
比如如果数据是200维的，那么就针对这200维做K-Means聚类。

文章图片
2018-12-25 19_29_38-【邹博_chinahadoop】机器学习升级版VII（七）.png
高斯分布的公式：
f(x) = 1/((2*π) 0.5*σ) * e -(x-μ)2/(2*σ2)
其中：μ是均值，σ 2是方差
如果f(x)是多元的，则得到多元高斯分布的概率密度函数：f(x) = (2*π) -n/2*(Σ -1) n*e -(1/2)*(x-μ)T*ΣT*(x-μ)
此处Σ为协方差矩阵，首先它是一个nxn的对称方阵。
这个Σ矩阵，在做混合高斯模型的时候，就出问题了。
比如现在做一个二聚类：
GMM:
N(μ1,Σ1）以及N(μ2,Σ2）
μ1与μ2都是n元的，而Σ1与Σ2都是矩阵:

如果矩阵是单位矩阵，如：
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1，
则σ 乘以单位矩阵得到：σ · I
则图形为球面的，即图中的Spherical柱状图。
参数有1个
如果矩阵是对角矩阵，如：
σ12 0 0 0
0 σ22 0 0
0 0 σ32 0
0 0 0 σ42
则得到diag柱状图
参数有n个
如果Σ1 = Σ2，则形成tied，即相互关联的
即图中的tied柱状图
参数有nxn个，准确的说是nx(n+1)个参数
Σ1 与 Σ2没有任何关联，我们求正常的EM算法
理论上有2倍的nx(n+1)个参数
即图中的full柱状图

只要是做混合高斯模型，基本都会涉及这四个参数
问答

问：协方差矩阵为什么是个对称阵？
答：因为这是定义。协方差矩阵是对称阵。
问：怎么看出来是球形？
答：如果随机变量是三元的，协方差矩阵如果三个方差都相等，即主轴，副轴与短轴都相等，得到球形
问：这四种情况怎么来的？
答：就是参数的设置。不管是做EM, sklearn的设置，还有隐马尔科夫模型，如果说隐马尔科夫模型符合高斯分布，那么就是高斯隐马尔科夫模型，那个模型里面，不同的隐变量，如果是方差也有方差是否相等，方差是不是对角阵等情况
问：参数设置时，这几种情况怎么选？
答：如果不知道怎么选，我们选full，即参数有2倍的nx(n+1)个，如果知道是对角阵，选diagonal，当然都试一下也无妨。
问：EM算法是无监督学习么？
答：EM算法可以看成无监督学习，虽然它是一种算法，是描述how而不是what。比如EM，MLE（最大似然估计），SGD（随机梯度下降），L-BFGS（拟牛顿）都是讲的how，即解决what的具体的方法。

模型选择的准则

文章图片
2018-12-25 20_12_36-【邹博_chinahadoop】机器学习升级版VII（七）.png
AIC解释：
负的对数似然，就可以作为目标函数；
但是我们不希望过拟合，所以需要在损失函数的前提下，加一个模型的复杂程度，比如模型的维度作为一个复杂标准。
哪个模型的这个值小，哪个模型最优。
即2k就成为了正则项
BIC解释：
样本多可能带来模型复杂度变化，如果两个模型，一个样本多，一个样本少，在结果相同的情况下，样本少的模型，看起来要好一些。
所以乘以与样本个数有关的项，是有道理的。
即(lnn)k
BIC也可以认为贝叶斯信息准则
BIC看相对大小才有意义，绝对大小没有意义

文章图片
2018-12-25 20_27_03-【邹博_chinahadoop】机器学习升级版VII（七）.png
很显然，当参数选择full的时候，错误率几乎就是0，并且BIC最小。即选择full参数的时候，模型是最优的。
问题：为什么上图右下角有一小块红色？
答：因为红色方差大。
问答

问：上述的例子说明什么？
答：对于Σ1与Σ的选择，引入更多的参数是否值得
问：平时计算EM的Σ不都是full类型么？
答：是的
问：样本的个数n为什么越大，BIC就大呢？为什么和样本个数n有关系？
答：能达到相同效果的时候，如果样本比别人多，那么模型就没有别人好。

文章图片
2018-12-25 20_37_50-【邹博_chinahadoop】机器学习升级版VII（七）.png

文章图片
2018-12-25 20_38_39-【邹博_chinahadoop】机器学习升级版VII（七）.png
如图，上图中三分类效果远远比二分类要差，所以可以加入一些先验知识。
如图，如果模型中的参数θ是未知定值，则可以通过最大似然估计(MLE)以及期望最大化(EM)去求。
如果θ也是变化的，且符合概率分布，即P(θ|α)，这个是先验分布，
对于样本y，只要给出x就能算出y的分布，且是对于θ的概率分布，这个是似然分布
P(θ|x, y)，则是属于后验分布

文章图片
2019-01-03 19_41_14-【邹博_chinahadoop】机器学习升级版VII（七）.png

文章图片
2019-01-03 19_48_16-【邹博_chinahadoop】机器学习升级版VII（七）.png
接着进行计算，如果θ有无穷多个，那么哪一个θ是最大的，就是我们想要求的：

文章图片
2019-01-03 19_51_04-【邹博_chinahadoop】机器学习升级版VII（七）.png θ这个值如何去求？
如图，后验分布，可以认为与似然分布 * 先验分布成正比。

文章图片
2019-01-03 19_53_45-【邹博_chinahadoop】机器学习升级版VII（七）.png 【机器学习笔记|机器学习笔记 - 20. EM算法实践（讲师（邹博））】如果θ是Dirichlet(狄利克雷)分布，可以演化为：
Dirichlet分布（参数为α+x） = 多项分布 * Dirichlet分布（参数为α）
如果α采样的值，拍脑袋选择1，10，或100，
假定α=1，得到θ1，θ2，θ3，。。。，θ100
从而分别得到（x1, y1），(x2, y2), (x3, y3), 。。。, (x100, y100)
这些是我们看到的样本数据。
其中每一个θ都是根据α采样得到的，即每一个θ都是一个随机变量，构成了一个随机过程，或者说构成了Dirichlet过程

文章图片
2019-01-03 20_07_03-【邹博_chinahadoop】机器学习升级版VII（七）.png
α取1的时候是最特殊的，即此时为均匀分布
再回头来看这张图，我们使用了符合高斯混合分布的模型，但是我们希望对参数做一个影响，那么使用Dirichlet过程+高斯混合分布模型，就得到DPGMM的模型，此时分类是合理的。
如图左边是正常的高斯混合模型，分类分错了；
右边是使用了贝叶斯的高斯混合模型，是特定的Dirichlet过程+高斯混合分布模型得到的结果，即使分类选择3，但是得到的结果也是正确的
在sk-learn中， DPGMM的相关类为：BayesianGaussianMixture

文章图片
2018-12-25 20_38_39-【邹博_chinahadoop】机器学习升级版VII（七）.png
相关的代码如下：

# !/usr/bin/python # -*- coding:utf-8 -*-import numpy as np from sklearn.mixture import GaussianMixture, BayesianGaussianMixture import scipy as sp import matplotlib as mpl import matplotlib.colors import matplotlib.pyplot as plt from matplotlib.patches import Ellipsedef expand(a, b, rate=0.05): d = (b - a) * rate return a-d, b+dmatplotlib.rcParams['font.sans-serif'] = ['SimHei'] matplotlib.rcParams['axes.unicode_minus'] = Falseif __name__ == '__main__': np.random.seed(0) cov1 = np.diag((1, 2)) N1 = 500 N2 = 300 N = N1 + N2 x1 = np.random.multivariate_normal(mean=(3, 2), cov=cov1, size=N1) m = np.array(((1, 1), (1, 3))) x1 = x1.dot(m) x2 = np.random.multivariate_normal(mean=(-1, 10), cov=cov1, size=N2) x = np.vstack((x1, x2)) y = np.array([0]*N1 + [1]*N2) n_components = 3# 绘图使用 colors = '#A0FFA0', '#2090E0', '#FF8080' cm = mpl.colors.ListedColormap(colors) x1_min, x1_max = x[:, 0].min(), x[:, 0].max() x2_min, x2_max = x[:, 1].min(), x[:, 1].max() x1_min, x1_max = expand(x1_min, x1_max) x2_min, x2_max = expand(x2_min, x2_max) x1, x2 = np.mgrid[x1_min:x1_max:500j, x2_min:x2_max:500j] grid_test = np.stack((x1.flat, x2.flat), axis=1)plt.figure(figsize=(6, 6), facecolor='w') plt.suptitle('GMM/DPGMM比较', fontsize=15)ax = plt.subplot(211) gmm = GaussianMixture(n_components=n_components, covariance_type='full', random_state=0) gmm.fit(x) centers = gmm.means_ covs = gmm.covariances_ print('GMM均值 = \n', centers) print('GMM方差 = \n', covs) y_hat = gmm.predict(x)grid_hat = gmm.predict(grid_test) grid_hat = grid_hat.reshape(x1.shape) plt.pcolormesh(x1, x2, grid_hat, cmap=cm) plt.scatter(x[:, 0], x[:, 1], s=20, c=y, cmap=cm, marker='o', edgecolors='#202020')clrs = list('rgbmy') for i, (center, cov) in enumerate(zip(centers, covs)): value, vector = sp.linalg.eigh(cov) width, height = value[0], value[1] v = vector[0] / sp.linalg.norm(vector[0]) angle = 180* np.arctan(v[1] / v[0]) / np.pi e = Ellipse(xy=center, width=width, height=height, angle=angle, color=clrs[i], alpha=0.5, clip_box = ax.bbox) ax.add_artist(e)ax1_min, ax1_max, ax2_min, ax2_max = plt.axis() plt.xlim((x1_min, x1_max)) plt.ylim((x2_min, x2_max)) plt.title('GMM', fontsize=15) plt.grid(b=True, ls=':', color='#606060')# DPGMM dpgmm = BayesianGaussianMixture(n_components=n_components, covariance_type='full', max_iter=1000, n_init=5, weight_concentration_prior_type='dirichlet_process', weight_concentration_prior=0.1) dpgmm.fit(x) centers = dpgmm.means_ covs = dpgmm.covariances_ print('DPGMM均值 = \n', centers) print('DPGMM方差 = \n', covs) y_hat = dpgmm.predict(x) print(y_hat)ax = plt.subplot(212) grid_hat = dpgmm.predict(grid_test) grid_hat = grid_hat.reshape(x1.shape) plt.pcolormesh(x1, x2, grid_hat, cmap=cm) plt.scatter(x[:, 0], x[:, 1], s=20, c=y, cmap=cm, marker='o', edgecolors='#202020')for i, cc in enumerate(zip(centers, covs)): if i not in y_hat: continue center, cov = cc value, vector = sp.linalg.eigh(cov) width, height = value[0], value[1] v = vector[0] / sp.linalg.norm(vector[0]) angle = 180* np.arctan(v[1] / v[0]) / np.pi e = Ellipse(xy=center, width=width, height=height, angle=angle, color='m', alpha=0.5, clip_box = ax.bbox) ax.add_artist(e) plt.xlim((x1_min, x1_max)) plt.ylim((x2_min, x2_max)) plt.title('DPGMM', fontsize=15) plt.grid(b=True, ls=':', color='#606060') plt.tight_layout(2, rect=(0, 0, 1, 0.95)) plt.show()

得到的结果如图：

文章图片
2019-01-03 20_15_36_20.EM_20.png

问：看来要讲贝叶斯啊
答：是的，要为下次LDA做铺垫
问：DPGMM其实就是一个主题和样本的分布作为权重，乘以主题和样本的高斯混合分布？
答：一定程度可以这样解释。
问：P(,; |)这三个符号一般是指什么呢？
答：P(x,y) <=> P(y, x)，这个代表x与y的联合分布
P(y; x)与P(y|x)是等价的，即x是条件，y是x的因变量
但是如果代入θ，就不一样了：
但是P(y; θ)，属于频率学派；其中θ为参数，θ是未知的定值；
而P(y|θ)，属于贝叶斯学派，则样本θ是未知的随机变量
问：都是同样阿尔法，怎么能采样出多个θ呢？
答：高斯分布中，如果均值为170，方差为10，取样可以为226,175，168多个值；同理，同样的阿尔法，也可以采样出多个θ
问：模型和里面的实现可以组合么？比如这种混合高斯模型，能用随机或批量梯度下降达到目的么？
答：其实是达不到的。因为混合高斯分布，其目标函数是有隐变量的存在，所以没办法对其直接求取梯度，只能固定隐变量求梯度；固定梯度求隐变量，二者不断迭代，最后才得到EM
问：均值为0和不为0有什么区别？效果会变化么？
答：均值是否为0只是一个解释，因为不为0的时候，我们总是会将其调整为0附近的。比如事先减均值。

文章图片
2019-02-04 15_07_45-【邹博_chinahadoop】机器学习升级版VII（七）.png
求导的过程很简单：
?h(p)/?p = n*p n-1*(1-p) (N-n) - p n*(N-n)*(1-p) (N-n-1)
假定导数为0，即：
?h(p)/?p = n*p n-1*(1-p) (N-n) - p n*(N-n)*(1-p) (N-n-1) = 0
则等式两边除以p n*(1-p) (N-n),得到：
?h(p)/?p = n*p -1 - (N-n)*(1-p) -1 = n/p - (N-n)/(1-p) = 0，
可得：p = n/N
下面看二项分布与先验举例：

文章图片
2019-02-04 15_47_50-【邹博_chinahadoop】机器学习升级版VII（七）.png
可以观察到，修正公式的分子各加了一个5，而这个5是Dirichlet(狄利克雷)分布的超参数α。

文章图片
2019-02-04 15_52_03-【邹博_chinahadoop】机器学习升级版VII（七）.png

问答
问：是不是这个课程代码都敲会，加上一个项目经验就OK?
答：这个看需要。比如现在的情况，机器学习是研究整个这套方式的一个根基。不用基础这个词，是因为以为其很简单，所以用根基这个词。有了这个根基之后，大家再去往上做其他应用，不会感觉困难。比如用卷积网络，最后一层我们使用SoftMax的全连接，还是用SVM，本质上是换损失函数。然后我们解释模型是否有效，都是可以用上的。
基础是够的，但是如果大家没有深度学习的应用实践，或者只有一个项目经验，还是不够的。所以需要实际项目进行反复验证与活学活用，或者参加竞赛，比赛也可以。
问：强化学习和机器学习的关联大么？强化学习未来的应用前景如何？
答：强化学习可能是近两三年的爆发点，可能是大公司玩的。需要的算力，比数据要强。比如飞翔的小鸟，或者行走的人，需要对当前的动作进行反馈，然后根据反馈的结果，去更正动作，并不断学习。算力要求非常高，前景要求是有的，但目前只能进行简单游戏、博弈、对抗这种内容，没法成为最主力的算法应用。也许不对，但最主流的还是在有监督应用。
问：算力是什么？可以简单理解为硬件速度么？
答：这个理解没问题的。

EM算法代码下面代码有自己实现的高斯混合模型，以及通过sk-learn库的高斯混合模型类直接实现的两种方式。
自己实现的高斯混合模型，其实就是实现期望最大化，即EM算法，公式如下：

文章图片
2019-02-04 17_19_35-【邹博_chinahadoop】机器学习升级版VII（七）.png

import numpy as np from scipy.stats import multivariate_normal from sklearn.mixture import GaussianMixture from mpl_toolkits.mplot3d import Axes3D import matplotlib as mpl import matplotlib.pyplot as plt from sklearn.metrics.pairwise import pairwise_distances_argminmpl.rcParams['font.sans-serif'] = ['SimHei'] mpl.rcParams['axes.unicode_minus'] = Falseif __name__ == '__main__': #style = 'sklearn' style = 'myself' np.random.seed(0) mu1_fact = (0, 0, 0) cov1_fact = np.diag((1, 2, 3)) # 根据实际情况生成一个多元正态分布矩阵，np.random.multivariate_normal # 参数就是高斯分布所需的均值与方差 # 第一个参数： mean：mean是多维分布的均值维度为1； # 第二个参数：cov：协方差矩阵，注意：协方差矩阵必须是对称的且需为半正定矩阵； # 第三个参数：size：指定生成的正态分布矩阵的维度 data1 = np.random.multivariate_normal(mu1_fact, cov1_fact, 400) print('data1 shape: {0}'.format(data1.shape)) mu2_fact = (2, 2, 1) # 方差对称且正定(positive-semidefinite): (4, 1, 3), (1, 2, 1), (3, 1, 4) cov2_fact = np.array(((4, 1, 3), (1, 2, 1), (3, 1, 4))) data2 = np.random.multivariate_normal(mu2_fact, cov2_fact, 100) print('data2 shape: {0}'.format(data2.shape))data = https://www.it610.com/article/np.vstack((data1, data2)) print('data shape: {0}'.format(data.shape)) y = np.array([True] * 400 + [False] * 100)if style == 'sklearn': g = GaussianMixture(n_components=2, covariance_type='full', tol=1e-6, max_iter=1000) g.fit(data) print('类别概率:\t', g.weights_[0]) print('均值:\n', g.means_, '\n') print('方差:\n', g.covariances_, '\n') mu1, mu2 = g.means_ sigma1, sigma2 = g.covariances_ else: num_iter = 100 n, d = data.shape # 随机指定 # mu1 = np.random.standard_normal(d) # print mu1 # mu2 = np.random.standard_normal(d) # print mu2 mu1 = data.min(axis=0) mu2 = data.max(axis=0) # 创建d行d列的单位矩阵（对角线为1，其余为0） sigma1 = np.identity(d) sigma2 = np.identity(d) pi = 0.5 # EM for i in range(num_iter): # E Step # 通过初始化的均值与方差，做多元的正态分布 norm1 = multivariate_normal(mu1, sigma1) norm2 = multivariate_normal(mu2, sigma2) # 概率密度 * pi tau1 = pi * norm1.pdf(data) tau2 = (1 - pi) * norm2.pdf(data) gamma = tau1 / (tau1 + tau2)# M Step mu1 = np.dot(gamma, data) / np.sum(gamma) mu2 = np.dot((1 - gamma), data) / np.sum((1 - gamma)) sigma1 = np.dot(gamma * (data - mu1).T, data - mu1) / np.sum(gamma) sigma2 = np.dot((1 - gamma) * (data - mu2).T, data - mu2) / np.sum(1 - gamma) pi = np.sum(gamma) / n print(i, ":\t", mu1, mu2) print('类别概率:\t', pi) print('均值:\t', mu1, mu2) print('方差:\n', sigma1, '\n\n', sigma2, '\n')# 预测分类 # multivariate_normal获得多元正态分布 norm1 = multivariate_normal(mu1, sigma1) norm2 = multivariate_normal(mu2, sigma2) # pdf: Probability density function，连续性概率分布函数 tau1 = norm1.pdf(data) tau2 = norm2.pdf(data)fig = plt.figure(figsize=(10, 5), facecolor='w') ax = fig.add_subplot(121, projection='3d') ax.scatter(data[:, 0], data[:, 1], data[:, 2], c='b', s=30, marker='o', edgecolors='k', depthshade=True) ax.set_xlabel('X') ax.set_ylabel('Y') ax.set_zlabel('Z') ax.set_title('原始数据', fontsize=15) ax = fig.add_subplot(122, projection='3d') # 求取点距离 order = pairwise_distances_argmin([mu1_fact, mu2_fact], [mu1, mu2], metric='euclidean') # order = pairwise_distances_argmin([mu1_fact, mu2_fact], [mu1, mu2], metric='cosine')# 通过欧式距离，将点分为两类 print(order) if order[0] == 0: c1 = tau1 > tau2 else: c1 = tau1 < tau2 c2 = ~c1 # 机器学习计算准确率的常用做法 # 原理：真实值是y，预测值是c1，相等则为True，否则为False。True为1，False为0 # 求均值则为：预测准确的数目/总数目，这不就是准确率么 acc = np.mean(y == c1) print('准确率：%.2f%%' % (100*acc)) ax.scatter(data[c1, 0], data[c1, 1], data[c1, 2], c='r', s=30, marker='o', edgecolors='k', depthshade=True) ax.scatter(data[c2, 0], data[c2, 1], data[c2, 2], c='g', s=30, marker='^', edgecolors='k', depthshade=True) ax.set_xlabel('X') ax.set_ylabel('Y') ax.set_zlabel('Z') ax.set_title('EM算法分类', fontsize=15) plt.suptitle('EM算法的实现', fontsize=18) plt.subplots_adjust(top=0.90) # plt.tight_layout() plt.show()

得到的图形界面：

文章图片
emdraw.png
如果将data2，那100个数据，均值设置为(5, 5, 5)，则分类效果更明显，如图：

文章图片
emdraw.png
当均值为5, 5, 5时，且自己实现高斯分布，输出如下。
可以发现，迭代24次之后，均值就不再变化了，可以称为模型收敛了。
真实值为[0, 0, 0]与[5, 5, 5]，计算的均值为：
[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ]，有差别，但是靠谱。
当然，也可以增加样本，使得结果更接近真实值的情况。比如将400, 100替换为4000, 1000

data1 shape: (400, 3) data2 shape: (100, 3) data shape: (500, 3)0 :[-0.029927490.091468150.03351835] [5.43632719 5.10518101 5.44044355] 1 :[-0.05773430.027438370.00975419] [5.23010882 5.07679997 5.22085292] 2 :[-0.0647070.00613346 -0.0011657 ] [5.13849255 5.04871819 5.14698533] 3 :[-0.0673045-0.00120726 -0.00558504] [5.09950565.03026551 5.1158224 ] 4 :[-0.068465-0.00420981 -0.00743412] [5.08222193 5.02089901 5.10147628] 5 :[-0.06899479 -0.00554239 -0.00824704] [5.07428887 5.01640324 5.09475058] 6 :[-0.06923923 -0.00615431 -0.00861565] [5.07059253 5.01427675 5.09158405] 7 :[-0.06935293 -0.00643926 -0.00878579] [5.06886018 5.01327497 5.09009261] 8 :[-0.06940609 -0.00657271 -0.00886508] [5.06804645 5.01280356 5.0893904 ] 9 :[-0.06943103 -0.00663537 -0.00890221] [5.06766389 5.01258178 5.0890599 ] 10 :[-0.06944274 -0.00666482 -0.00891964] [5.06748396 5.01247744 5.08890438] 11 :[-0.06944825 -0.00667867 -0.00892783] [5.06739933 5.01242836 5.08883121] 12 :[-0.06945084 -0.00668519 -0.00893168] [5.06735951 5.01240527 5.08879678] 13 :[-0.06945206 -0.00668825 -0.00893349] [5.06734079 5.01239441 5.08878059] 14 :[-0.06945263 -0.00668969 -0.00893434] [5.06733197 5.01238935.08877297] 15 :[-0.0694529-0.00669037 -0.00893474] [5.06732783 5.01238695.08876938] 16 :[-0.06945302 -0.00669069 -0.00893493] [5.06732588 5.01238577 5.0887677 ] 17 :[-0.06945308 -0.00669084 -0.00893502] [5.06732496 5.01238523 5.0887669 ] 18 :[-0.06945311 -0.00669091 -0.00893506] [5.06732453 5.01238498 5.08876653] 19 :[-0.06945313 -0.00669095 -0.00893508] [5.06732433 5.01238487 5.08876635] 20 :[-0.06945313 -0.00669096 -0.00893509] [5.06732423 5.01238481 5.08876627] 21 :[-0.06945314 -0.00669097 -0.00893509] [5.06732419 5.01238478 5.08876623] 22 :[-0.06945314 -0.00669097 -0.00893509] [5.06732417 5.01238477 5.08876621] 23 :[-0.06945314 -0.00669097 -0.0089351 ] [5.06732416 5.01238477 5.08876621] 24 :[-0.06945314 -0.00669097 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 25 :[-0.06945314 -0.00669097 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 26 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 27 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 28 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 29 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 30 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 31 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 32 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 33 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 34 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 35 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 36 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 37 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 38 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 39 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 40 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 41 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 42 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 43 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 44 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 45 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 46 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 47 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 48 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 49 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 50 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 51 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 52 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 53 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 54 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 55 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 56 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 57 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 58 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 59 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 60 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 61 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 62 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 63 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 64 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 65 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 66 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 67 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 68 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 69 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 70 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 71 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 72 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 73 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 74 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 75 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 76 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 77 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 78 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 79 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 80 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 81 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 82 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 83 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 84 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 85 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 86 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 87 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 88 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 89 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 90 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 91 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 92 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 93 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 94 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 95 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 96 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 97 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 98 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 99 :[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 类别概率:0.7987220297951044 均值:[-0.06945314 -0.00669098 -0.0089351 ] [5.06732415 5.01238476 5.0887662 ] 方差: [[ 0.87148101 -0.056424940.03198856] [-0.056424942.09700921 -0.12547629] [ 0.03198856 -0.125476292.745459]] [[4.08142083 0.79087313 3.107469] [0.79087313 1.79995257 0.75954681] [3.1074690.75954681 4.04331614]] [0 1] 准确率：98.60%

问答
问：协方差一定是对称的么？
答：是的，协方差一定是对称的
问：np.identity是什么意思？
答：创建单位矩阵，即对角线为1，其余值为0

GMM代码实现对应业务为性别-身高-体重数据
通过高斯混合模型，预测身高与体重所属的性别

# !/usr/bin/python # -*- coding:utf-8 -*-import numpy as np from sklearn.mixture import GaussianMixture from sklearn.model_selection import train_test_split import matplotlib as mpl import matplotlib.colors import matplotlib.pyplot as pltmpl.rcParams['font.sans-serif'] = ['SimHei'] mpl.rcParams['axes.unicode_minus'] = False # from matplotlib.font_manager import FontProperties # font_set = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=15) # fontproperties=font_setdef expand(a, b): d = (b - a) * 0.05 return a-d, b+dif __name__ == '__main__': data = https://www.it610.com/article/np.loadtxt('./HeightWeight.csv', dtype=np.float, delimiter=',', skiprows=1) print(data.shape) y, x = np.split(data, [1, ], axis=1) x, x_test, y, y_test = train_test_split(x, y, train_size=0.6, random_state=0) gmm = GaussianMixture(n_components=2, covariance_type='full', random_state=0) x_min = np.min(x, axis=0) x_max = np.max(x, axis=0) gmm.fit(x) print('均值 = \n', gmm.means_) print('方差 = \n', gmm.covariances_) y_hat = gmm.predict(x) y_test_hat = gmm.predict(x_test) change = (gmm.means_[0][0] > gmm.means_[1][0]) if change: z = y_hat == 0 y_hat[z] = 1 y_hat[~z] = 0 z = y_test_hat == 0 y_test_hat[z] = 1 y_test_hat[~z] = 0 acc = np.mean(y_hat.ravel() == y.ravel()) acc_test = np.mean(y_test_hat.ravel() == y_test.ravel()) acc_str = '训练集准确率：%.2f%%' % (acc * 100) acc_test_str = '测试集准确率：%.2f%%' % (acc_test * 100) print(acc_str) print(acc_test_str)cm_light = mpl.colors.ListedColormap(['#FF8080', '#77E0A0']) cm_dark = mpl.colors.ListedColormap(['r', 'g']) x1_min, x1_max = x[:, 0].min(), x[:, 0].max() x2_min, x2_max = x[:, 1].min(), x[:, 1].max() x1_min, x1_max = expand(x1_min, x1_max) x2_min, x2_max = expand(x2_min, x2_max) x1, x2 = np.mgrid[x1_min:x1_max:500j, x2_min:x2_max:500j] grid_test = np.stack((x1.flat, x2.flat), axis=1) grid_hat = gmm.predict(grid_test) grid_hat = grid_hat.reshape(x1.shape) if change: z = grid_hat == 0 grid_hat[z] = 1 grid_hat[~z] = 0 plt.figure(figsize=(7, 6), facecolor='w') plt.pcolormesh(x1, x2, grid_hat, cmap=cm_light) plt.scatter(x[:, 0], x[:, 1], s=50, c=y.ravel(), marker='o', cmap=cm_dark, edgecolors='k') plt.scatter(x_test[:, 0], x_test[:, 1], s=60, c=y_test.ravel(), marker='^', cmap=cm_dark, edgecolors='k')p = gmm.predict_proba(grid_test) print(p) p = p[:, 0].reshape(x1.shape) CS = plt.contour(x1, x2, p, levels=(0.1, 0.5, 0.8), colors=list('rgb'), linewidths=2) plt.clabel(CS, fontsize=12, fmt='%.1f', inline=True) ax1_min, ax1_max, ax2_min, ax2_max = plt.axis() xx = 0.95*ax1_min + 0.05*ax1_max yy = 0.05*ax2_min + 0.95*ax2_max plt.text(xx, yy, acc_str, fontsize=12) yy = 0.1*ax2_min + 0.9*ax2_max plt.text(xx, yy, acc_test_str, fontsize=12) plt.xlim((x1_min, x1_max)) plt.ylim((x2_min, x2_max)) plt.xlabel('身高(cm)', fontsize=13) plt.ylabel('体重(kg)', fontsize=13) plt.title('EM算法估算GMM的参数', fontsize=15) plt.grid(b=True, ls=':', color='#606060') plt.tight_layout(2) plt.show()

其中HeightWeight.csv的数据如下，直接将其拷贝到文本文件，然后保存为文件名为HeightWeight.csv的文件即可

Sex,Height(cm),Weight(kg) 0,156,50 0,160,60 0,162,54 0,162,55 0,160.5,56 0,160,53 0,158,55 0,164,60 0,165,50 0,166,55 0,158,47.5 0,161,49 0,169,55 0,161,46 0,160,45 0,167,44 0,155,49 0,154,57 0,172,52 0,155,56 0,157,55 0,165,65 0,156,52 0,155,50 0,156,56 0,160,55 0,158,55 0,162,70 0,162,65 0,155,57 0,163,70 0,160,60 0,162,55 0,165,65 0,159,60 0,147,47 0,163,53 0,157,54 0,160,55 0,162,48 0,158,60 0,155,48 0,165,60 0,161,58 0,159,45 0,163,50 0,158,49 0,155,50 0,162,55 0,157,63 0,159,49 0,152,47 0,156,51 0,165,49 0,154,47 0,156,52 0,162,48 1,162,60 1,164,62 1,168,86 1,187,75 1,167,75 1,174,64 1,175,62 1,170,65 1,176,73 1,169,58 1,178,54 1,165,66 1,183,68 1,171,61 1,179,64 1,172,60 1,173,59 1,172,58 1,175,62 1,160,60 1,160,58 1,160,60 1,175,75 1,163,60 1,181,77 1,172,80 1,175,73 1,175,60 1,167,65 1,172,60 1,169,75 1,172,65 1,175,72 1,172,60 1,170,65 1,158,59 1,167,63 1,164,61 1,176,65 1,182,95 1,173,75 1,176,67 1,163,58 1,166,67 1,162,59 1,169,56 1,163,59 1,163,56 1,176,62 1,169,57 1,173,61 1,163,59 1,167,57 1,176,63 1,168,61 1,167,60 1,170,69

图形示例如下：

文章图片
2019-02-04 17_41_49-Start.png

问：from sklearn.metrics.pairwise import pairwise_distances_argmin，这个是干嘛的？
答：是用于计算任意的两个值里面，谁和谁是最小的。比如：order = pairwise_distances_argmin([mu1_fact, mu2_fact], [mu1, mu2], metric='euclidean')，返回的值是[0,1]，表明mu1与mu1_fact最近，mu2与mu2_fact最近。换句话说，我们做的顺序是做对了。

通过GMM实现鸢尾花分类通过高斯混合模型，对鸢尾花数据做分类

# !/usr/bin/python # -*- coding:utf-8 -*-import numpy as np import pandas as pd from sklearn.mixture import GaussianMixture import matplotlib as mpl import matplotlib.colors import matplotlib.pyplot as plt from sklearn.metrics.pairwise import pairwise_distances_argminmpl.rcParams['font.sans-serif'] = ['SimHei'] mpl.rcParams['axes.unicode_minus'] = Falseiris_feature = '花萼长度', '花萼宽度', '花瓣长度', '花瓣宽度'def expand(a, b, rate=0.05): d = (b - a) * rate return a-d, b+dif __name__ == '__main__': path = '..\9.Regression\iris.data' data = https://www.it610.com/article/pd.read_csv(path, header=None) x_prime = data[np.arange(4)] y = pd.Categorical(data[4]).codesn_components = 3 feature_pairs = [[0, 1], [0, 2], [0, 3], [1, 2], [1, 3], [2, 3]] plt.figure(figsize=(8, 6), facecolor='w') for k, pair in enumerate(feature_pairs, start=1): x = x_prime[pair] m = np.array([np.mean(x[y == i], axis=0) for i in range(3)])# 均值的实际值 print('实际均值 = \n', m)gmm = GaussianMixture(n_components=n_components, covariance_type='full', random_state=0) gmm.fit(x) print('预测均值 = \n', gmm.means_) print('预测方差 = \n', gmm.covariances_) y_hat = gmm.predict(x) print(y_hat) order = pairwise_distances_argmin(m, gmm.means_, axis=1, metric='euclidean') print(order) print('顺序：\t', order)n_sample = y.size n_types = 3 change = np.empty((n_types, n_sample), dtype=np.bool) for i in range(n_types): change[i] = y_hat == order[i] for i in range(n_types): y_hat[change[i]] = i acc = '准确率：%.2f%%' % (100*np.mean(y_hat == y)) print(acc)cm_light = mpl.colors.ListedColormap(['#FF8080', '#77E0A0', '#A0A0FF']) cm_dark = mpl.colors.ListedColormap(['r', 'g', '#6060FF']) x1_min, x2_min = x.min() x1_max, x2_max = x.max() x1_min, x1_max = expand(x1_min, x1_max) x2_min, x2_max = expand(x2_min, x2_max) x1, x2 = np.mgrid[x1_min:x1_max:200j, x2_min:x2_max:200j] grid_test = np.stack((x1.flat, x2.flat), axis=1) grid_hat = gmm.predict(grid_test)change = np.empty((n_types, grid_hat.size), dtype=np.bool) for i in range(n_types): change[i] = grid_hat == order[i] for i in range(n_types): grid_hat[change[i]] = igrid_hat = grid_hat.reshape(x1.shape) plt.subplot(2, 3, k) plt.pcolormesh(x1, x2, grid_hat, cmap=cm_light) plt.scatter(x[pair[0]], x[pair[1]], s=20, c=y, marker='o', cmap=cm_dark, edgecolors='k') xx = 0.95 * x1_min + 0.05 * x1_max yy = 0.1 * x2_min + 0.9 * x2_max plt.text(xx, yy, acc, fontsize=10) plt.xlim((x1_min, x1_max)) plt.ylim((x2_min, x2_max)) plt.xlabel(iris_feature[pair[0]], fontsize=11) plt.ylabel(iris_feature[pair[1]], fontsize=11) plt.grid(b=True, ls=':', color='#606060') plt.suptitle('EM算法无监督分类鸢尾花数据', fontsize=14) plt.tight_layout(1, rect=(0, 0, 1, 0.95)) plt.show()

图例如下：

文章图片
2019-02-04 18_17_05-Figure 1.png 绘制高斯混合模型的等值线

# !/usr/bin/python # -*- coding:utf-8 -*-import numpy as np from sklearn.mixture import GaussianMixture import scipy as sp import matplotlib as mpl import matplotlib.colors import matplotlib.pyplot as plt from matplotlib.patches import Ellipse import warningsdef expand(a, b, rate=0.05): d = (b - a) * rate return a-d, b+dif __name__ == '__main__': warnings.filterwarnings(action='ignore', category=RuntimeWarning) np.random.seed(0) cov1 = np.diag((1, 2)) N1 = 500 N2 = 300 N = N1 + N2 x1 = np.random.multivariate_normal(mean=(3, 2), cov=cov1, size=N1) m = np.array(((1, 1), (1, 3))) x1 = x1.dot(m) x2 = np.random.multivariate_normal(mean=(-1, 10), cov=cov1, size=N2) x = np.vstack((x1, x2)) y = np.array([0]*N1 + [1]*N2)gmm = GaussianMixture(n_components=2, covariance_type='full', random_state=0) gmm.fit(x) centers = gmm.means_ covs = gmm.covariances_ print('GMM均值 = \n', centers) print('GMM方差 = \n', covs) y_hat = gmm.predict(x)colors = '#A0FFA0', '#E080A0', levels = 10 cm = mpl.colors.ListedColormap(colors) x1_min, x1_max = x[:, 0].min(), x[:, 0].max() x2_min, x2_max = x[:, 1].min(), x[:, 1].max() x1_min, x1_max = expand(x1_min, x1_max) x2_min, x2_max = expand(x2_min, x2_max) x1, x2 = np.mgrid[x1_min:x1_max:500j, x2_min:x2_max:500j] grid_test = np.stack((x1.flat, x2.flat), axis=1) print(gmm.score_samples(grid_test)) grid_hat = -gmm.score_samples(grid_test) grid_hat = grid_hat.reshape(x1.shape) plt.figure(figsize=(7, 6), facecolor='w') ax = plt.subplot(111) cmesh = plt.pcolormesh(x1, x2, grid_hat, cmap=plt.cm.Spectral) plt.colorbar(cmesh, shrink=0.9) CS = plt.contour(x1, x2, grid_hat, levels=np.logspace(0, 2, num=levels, base=10), colors='w', linewidths=1) plt.clabel(CS, fontsize=9, inline=True, fmt='%.1f') plt.scatter(x[:, 0], x[:, 1], s=30, c=y, cmap=cm, marker='o', edgecolors='#202020')for i, cc in enumerate(zip(centers, covs)): center, cov = cc value, vector = sp.linalg.eigh(cov) width, height = value[0], value[1] v = vector[0] / sp.linalg.norm(vector[0]) angle = 180* np.arctan(v[1] / v[0]) / np.pi e = Ellipse(xy=center, width=width, height=height, angle=angle, color='m', alpha=0.5, clip_box = ax.bbox) ax.add_artist(e)plt.xlim((x1_min, x1_max)) plt.ylim((x2_min, x2_max)) mpl.rcParams['font.sans-serif'] = ['SimHei'] mpl.rcParams['axes.unicode_minus'] = False plt.title('GMM似然函数值', fontsize=15) plt.grid(b=True, ls=':', color='#606060') plt.tight_layout(2) plt.show()

图例如下：

文章图片
2019-02-04 18_29_19-Start.png

问答
问：DPGMM选的k是不是要尽量小？
答：不一定，与k值选择是否小没关系。
问：矩阵运算不是不能用交换律么？怎么直接交换了？
答：是对这个代码：sigma1 = np.dot(gamma * (data - mu1).T, data - mu1) / np.sum(gamma)，矩阵不能交换，但是标量值就可以进行交换。

EM算法内容完结