0和1

上一篇博客主要介绍了逻辑回归的理论知识，这篇博客咱们用Python机器学习包sklearn中的LogisticRegression做一个分类的实例。 数据还是学生样本，只有两个特征，分别是两门课的分数score1和score2，类标号y表示这个学生是否能被录取。先上分类效果图： 完整的Python代码如下： 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 # -*- coding: utf-8 -*- """ Created on Wed Nov 08 17:49:41 2017 @author: zhenlin """ import numpy as np import pandas as pd from sklearn.cross_validation import train_test_split from sklearn.linear_model import LogisticRegression import matplotlib.pyplot as plt from sklearn.metrics import precision_recall_curve from sklearn.metrics import classification_report # 1. 构造数据 sample_number = 200 # 第一个高斯分布参数 mean1 = [0, 4] # 两个维度上的均值 cov1 = [[5, 3], [3, 10]] # 两个维度的协方差矩阵，必须满足对称半正定 # 第二个高斯分布参数 mean2 = [7, 5] cov2 = [[7, 2], [2, 15]] # 从两个二元高斯分布中随机采样数据点 class1_x1, class1_x2 = np.random.multivariate_normal(mean1, cov1, sample_number).T # .T表示转置 class2_x1, class2_x2 = np.random.multivariate_normal(mean2, cov2, sample_number).T # 两个高斯分布对应两个类标号 data = [[class1_x1[i],class1_x2[i],0] for i in range(sample_number)]+[[class2_x1[i],class2_x2[i],1] for i in range(sample_number)] # 填充到pandas中 data = pd.DataFrame(data,columns=['score1','score2','result']) score_data = data[['score1','score2']] result_data = data['result'] # 2. 训练模型 average_precision = 0 # 平均准确度 iters = 10 # 交叉验证次数 for i in xrange(iters): # 数据划分，80%用于训练，20%用于预测 x_train, x_test, y_train, y_test = train_test_split(score_data, result_data, test_size = 0.2) # 构造默认逻辑回归模型 model = LogisticRegression() # 训练 model.fit(x_train, y_train) # 预测 predict_y = model.predict(x_test) # 计算测试集上的准确度 average_precision += np.mean(predict_y == y_test) average_precision /= iters # 3. 绘制分类面 - 法1 x1_min, x1_max = score_data['score1'].min() - .5, score_data['score1'].max() + .5 def generate_face(prob): y = -np.log(1.0 / prob - 1.0) n = 500 x1 = np.linspace(x1_min, x1_max, n) # w1x1+w2x2+b=y x2 = (-model.coef_[0][0] / float(model.coef_[0][1])) * x1 + (y - model.intercept_) / float(model.coef_[0][1]) return x1, x2 pos_data = data[data['result'] == 1] neg_data = data[data['result'] == 0] plt.scatter(x = pos_data['score1'], y = pos_data['score2'], color = 'black', marker = 'o') plt.scatter(x = neg_data['score1'], y = neg_data['score2'], color = 'red', marker = '*') face_04_x1, face_04_x2 = generate_face(0.4) face_05_x1, face_05_x2 = generate_face(0.5) face_06_x1, face_06_x2 = generate_face(0.6) plt.plot(face_04_x1, face_04_x2) plt.plot(face_05_x1, face_05_x2) plt.plot(face_06_x1, face_06_x2) plt.xlim(score_data['score1'].min(), score_data['score1'].max()) plt.ylim(score_data['score2'].min(), score_data['score2'].max()) plt.xlabel('score1') plt.ylabel('score2') plt.legend(['prob_threshold = 0.4', 'prob_threshold = 0.5', 'prob_threshold = 0.6'], loc='center left', bbox_to_anchor=(1, 0.865)) plt.show() # 4. 绘制分类面 - 法2 pos_data = data[data['result'] == 1] neg_data = data[data['result'] == 0] h = 0.02 s1_min, s1_max = score_data['score1'].min() - .5, score_data['score1'].max() + .5 s2_min, s2_max = score_data['score2'].min() - .5, score_data['score2'].max() + .5 # 生成s1在[s1_min, s1_max]，且s2在[s2_min, s2_max]的网格数据点 # meshgrid含义参见：http://blog.sciencenet.cn/blog-791749-675394.html s1, s2 = np.meshgrid(np.arange(s1_min, s1_max, h), np.arange(s2_min, s2_max, h)) # 把两个坐标的值按列拼在一起构成二维数据点 Z = model.predict(np.c_[s1.ravel(), s2.ravel()]) # 绘制边界和散点 Z = Z.reshape(s1.shape) # 坐标点是(s1[i], s2[i])，对应颜色是Z[i]，颜色主题使用plt.cm.Paired plt.pcolormesh(s1, s2, Z, cmap = plt.cm.Paired) plt.scatter(x = pos_data['score1'], y = pos_data['score2'], color = 'black', marker = 'o') plt.scatter(x = neg_data['score1'], y = neg_data['score2'], color = 'red', marker = '*') plt.xlim(s1.min(), s1.max()) plt.ylim(s2.min(), s2.max()) plt.xlabel('score1') plt.ylabel('score2') plt.show() # 5. 评估模型 # 对于测试数据，模型输出1的概率 answer = model.predict_proba(x_test)[:,1] # 计算不同概率阈值下的P和R precision, recall, thresholds = precision_recall_curve(y_test, answer) # prob > 0.5的报告为1 report = answer > 0.5 print(classification_report(y_test, report, target_names = ['neg', 'pos'])) print('average precision: %f'%average_precision) # 6. 绘制PRC曲线 # step阶跃图，在点(recall[i],precision[i])进行跳变 plt.step(recall, precision, color='b', alpha=0.2, where='post') # 对PRC下方填充颜色 plt.fill_between(recall, precision, step='post', alpha=0.2, color='b') plt.xlabel('Recall') plt.ylabel('Precision') plt.ylim([0.0, 1.05]) plt.xlim([0.0, 1.0]) plt.title('2-class Precision-Recall curve') plt.show() 下面将逐模块介绍代码细节，大神可以略过。 ...

最近实验室在组织学习NG的机器学习视频，我也跟进了一下。讲到逻辑回归那一课，一直想不明白，逻辑回归到底是怎样分类的？逻辑回归的分类面在哪里？逻辑回归有像SVM的max margin那样清晰的推导过程吗？为什么需要Sigmoid函数？今天就让我们来一窥逻辑回归的始末。 假设有一堆学生样本，为了简单起见，只有两个特征，分别是两门课的分数score1和score2，类标号y表示这个学生是否能被录取。可视化如下图，黑点表示y=1即被录取，红点表示y=0即未被录取，可以看到黑点处在score1和score2都比较高的区域。我们的任务就是给定这些训练样本，需要确定一个分类面来划分这两类数据。 分类面的实质就是\(y=\mathbf{w^T x}+b\)，其中\(\mathbf{w}\)和\(\mathbf{x}\)都是向量，对应到本例中，展开为\(y=w_1x_1+w_2x_2+b\)。所以，寻找分类面的过程就是寻找倾斜度\(\mathbf{w}\)和截距\(b\)。 因为分类的结果是离散的，只有0或者1，可以用感知机来分类。即 但是怎样找这里的\(\mathbf{w}\)和\(b\)使得分类结果最好呢，我们需要定义一个优化的目标函数或者说损失函数。这里只能定义为分类错误的个数，只能一点点挪动超平面来找分类错误最少的超平面了，即只能用暴力枚举的方法来找\(\mathbf{w}\)和\(b\)。 另一种方法是令我们的分类面算出来的值和真实标号越接近越好，即最小化误差\((\mathbf{w^T}x^{(i)}+b-y^{(i)})^2\)，然后通过梯度下降求解\(\mathbf{w}\)和\(b\)。 但是这会有一个问题，对于上图数据，可以求到一个比较好的分类面，如下左图。但是如果数据中出现了如下右图右边的一个离群点或者噪声，为了最小化这个点的误差，即\((\mathbf{w^T}x^{(i)}+b-1)^2\)，导致分类面往右偏了，这一偏直接导致很多y=1的样本分类错误。所以这种最小化误差的方法对离群点很敏感。 假设分类超平面还是\(f(\mathbf{x})=\mathbf{w^T x}+b\)，我们希望的分类效果是这样的：\(f(\mathbf{x})=0\)是分类面；\(f(\mathbf{x})>0\)分类为1，且不管\(f(\mathbf{x})\)多大，都分为1；\(f(\mathbf{x})<0\)分类为0，且不管\(f(\mathbf{x})\)多小，都分为0。 因为类标号是离散的{0,1}，所以想到把\(f(\mathbf{x})\)映射到[0,1]之间，即\(g(f(\mathbf{x}))\)。为了满足上述条件，\(g(f(\mathbf{x}))\)需要满足： \(g(0)=0.5\)，即在分类面上无法判断类标号是0还是1 当\(f(\mathbf{x})>0\)时，\(g(f(\mathbf{x}))>0.5\) 当\(f(\mathbf{x})\rightarrow+\infty\)，\(g(f(\mathbf{x}))\rightarrow 1\)，且\(g'(f(\mathbf{x}))\rightarrow 0\)，即对于上右图右边的离群点，分类为1，且导数趋近于0，表示对其不敏感 当\(f(\mathbf{x})<0\)时，\(g(f(\mathbf{x}))<0.5\) 当\(f(\mathbf{x})\rightarrow-\infty\)，\(g(f(\mathbf{x}))\rightarrow 0\)，且\(g'(f(\mathbf{x}))\rightarrow 0\)，和第3点类似 满足上述性质的函数之一就是Sigmoid函数，其定义域为\([-\infty,+\infty]\)，值域为[0,1]，正好把原始的函数结果\(f(\mathbf{x})\)映射到了[0,1]，而概率的取值范围正好是[0,1]，所以Sigmoid函数正好可以作为分类为1的概率。 所以逻辑回归最终的形式就是： $$g(\mathbf{x})=\frac{1}{1+e^{-(\mathbf{w^T x}+b)}}$$分类面依然还是\(f(\mathbf{x})=\mathbf{w^T x}+b=0\)，因为\(f(\mathbf{x})=0\)时，\(g(\mathbf{x})=0.5\)，正好满足上述条件1。 Sigmoid函数的另一个优点是，它把原始函数值映射到了概率空间，这样就可以用极大似然的方法求解参数\(\mathbf{w}\)和\(b\)。下面用参数\(\mathbf{\theta}\)代表参数\(\mathbf{w}\)和\(b\)，用\(h_{\mathbf{\theta}}(\mathbf{x})\)代表\(g(f(\mathbf{x}))\)。则有： $$P(y=1|\mathbf{x};\mathbf{\theta})=h_{\mathbf{\theta}}(\mathbf{x})$$$$P(y=0|\mathbf{x};\mathbf{\theta})=1-h_{\mathbf{\theta}}(\mathbf{x})$$合并成一个式子就是： $$P(y|\mathbf{x};\mathbf{\theta})=h_{\mathbf{\theta}}(\mathbf{x})^y(1-h_{\mathbf{\theta}}(\mathbf{x}))^{1-y}$$由于所有样本独立同分布（I.I.D.），似然函数就是 $$L(\mathbf{\theta})=P(\mathbf{y}|X;\mathbf{\theta})=\prod\limits_{i}P(y^{(i)}|\mathbf{x}^{(i)};\mathbf{\theta})=\prod\limits_{i}h_{\mathbf{\theta}}(\mathbf{x}^{(i)})^{y^{(i)}}(1-h_{\mathbf{\theta}}(\mathbf{x}^{(i)}))^{1-y^{(i)}}$$最大化似然的含义就是，在给定样本\(X\)的情况下，我们想找一个参数\(\mathbf{\theta}\)，使得观测到类标号\(\mathbf{y}\)的概率最大。 最大化似然等价于最大化log似然，log展开之后就是： $$l(\mathbf{\theta})=logL(\mathbf{\theta})=\sum\limits_{i}y^{(i)}logh_{\mathbf{\theta}}(\mathbf{x}^{(i)})+(1-y^{(i)})log(1-h_{\mathbf{\theta}}(\mathbf{x}^{(i)}))$$而很巧的是，最大化log似然，其实等效于最小化log对数损失。对于单个样本，损失为： $$cost(h_{\theta}(\mathbf{x}),y) = \begin{cases}-log(h_{\theta}(\mathbf{x})) & \text {if y=1} \\ -log(1-h_{\theta}(\mathbf{x})) & \text{if y=0} \end{cases}$$即如果正确类标号是1，但算出来的\(h_{\theta}(\mathbf{x})\)很接近0的话，则损失\(-log(h_{\theta}(\mathbf{x}))\)就会很大。类标号为0的情况类似。把这两种情况合成一个式子就是： $$cost(h_{\theta}(\mathbf{x}),y) = -ylog(h_{\theta}(\mathbf{x})) – (1-y)log(1-h_{\theta}(\mathbf{x}))$$所有样本的损失之和就是： $$J(\mathbf{\theta})=cost(h_{\theta}(X),\mathbf{y}) = \sum\limits_{i}-y^{(i)}logh_{\mathbf{\theta}}(\mathbf{x}^{(i)})-(1-y^{(i)})log(1-h_{\mathbf{\theta}}(\mathbf{x}^{(i)}))$$所以最大化对数似然\(\max l(\mathbf{\theta})\)和最小化对数损失\(\min J(\mathbf{\theta})\)是等效的！最优化求解的方法用梯度下降和牛顿法都可以。 和Sigmoid很像的函数还有很多，比如y=arctan(x)也可以，不过Sigmoid有一个很好的特点，即它的导数可以由自身算出来，\(f'(x)=f(x)(1-f(x))\)。 传统的逻辑回归只能处理二分类问题，那么怎样将其扩展到解决多分类问题呢？有两种方法，一种方法是建立k个二元分类器，比如类标号有A,B,C，则建立3个二元分类器，分别是1）A和非A；2）B和非B；3）C和非C。训练每个2元分类器时，都对所有的数据重新标号为0或1，这样共需要训练k个二元分类器。 还有一种方法是直接将二元逻辑回归推广到多元回归，即Softmax回归，有关Softmax回归的内容，请参考此博客，非常详细。简单来说，二元逻辑回归的假设函数是：多元Softmax回归的假设函数在形式上有所不同：其中是模型的参数。请注意这一项对概率分布进行归一化，使得所有概率之和为1。 Softmax回归的损失函数如下，其实就是logistic回归损失函数的推广： 二元逻辑回归是多元Softmax回归在k=2时的特例，令k=2并利用Softmax回归参数冗余的特点，可以得到一般形式的二元逻辑回归假设函数，具体可以看原博客。 面对一个多元分类问题，是选择Softmax回归还是k个二元分类器呢，这取决于你的类别之间是否互斥，如果互斥，可以用Softmax回归，否则，请用k个二元分类器。 这就是逻辑回归的理论知识，下一篇博客将介绍逻辑回归在Python中的应用。 参考： http://www.cnblogs.com/sparkwen/p/3441197.html https://tech.meituan.com/intro_to_logistic_regression.html http://blog.csdn.net/bitcarmanlee/article/details/51165444 http://ufldl.stanford.edu/wiki/index.php/Softmax%E5%9B%9E%E5%BD%92

上一篇博客简单介绍了Linux平台下不同性能分析工具的特点，最后是Google出品的gperftools胜出。今天我们将举一个用gperftools对链接库进行性能分析的例子。 假设我们创建了一个TestGperftools项目，其中包括一个动态链接库模块complex和主程序main，build文件夹用于存放动态链接库和可执行程序，文件结构如下所示。 1 2 3 4 5 6 7 8 9 10 11 12 13 czl@ubuntu:~/czl/TestGperftools$ tree . ├── build ├── complex │ ├── complex.cpp │ ├── complex.h │ └── Makefile ├── main │ ├── main.cpp │ └── Makefile └── Makefile 3 directories, 6 files 其中complex借用了这篇博客的第三个公式计算圆周率π，主要用来模拟一个很耗时的操作。该模块会生成一个libcomplex.so动态链接库供外部调用。其complex.cpp函数如下： ...

这半年一直在研究pLink 2的加速策略，了解了相关的性能分析工具，现记录如下。 对软件进行加速的技术路线一般为：首先对代码进行性能分析（ performance analysis，也称为 profiling），然后针对性能瓶颈提出优化方案，最后在大数据集上评测加速效果。所以进行性能优化的前提就是准确地测量代码中各个模块的时间消耗。听起来很简单，不就是测量每行代码的运行时间吗，直接用time_t t=clock();不就好了，但是事情并没有那么简单。如果只进行粗粒度的性能分析，比如测量几个大的模块的运行时间，clock()还比较准确，但是如果测量的是运算量比较小的函数调用，而且有大量的小函数调用，clock()就不太准确了。 比如下面的一段代码，我最开始的性能分析方法是在fun1()~fun3()前后添加time_t t=clock()，然后作差求和的。但是3个fun()加起来的时间居然不等于整个while循环的时间，有将近50%的时间不翼而飞了！ 1 2 3 4 5 6 7 8 9 10 11 while (true) { if (fun1()) { for (int i = 0; i < k; ++k) { if (flag1) { fun2(); } } } else { fun3(); } } 一种可能是while循环以及内部的for循环本身占用了较多的时间，但是好像不太可能呀。还有一种可能是clock()测量有误，time_t只能精确到秒，clock_t只能精确到毫秒，如果一次fun*()的时间太短，导致一次测量几乎为0，那么多次的while和for循环调用累加的时间也几乎为0，导致实际测量到的fun*()时间远小于真实时间。所以自己用代码进行性能分析可能会有较大的误差，最好借助已有的性能分析工具。 ...

今天阅读《C++ Primer, 5e》的第二章，介绍C++的基本内置类型，觉得有一些平时工作容易出错的知识点，现摘录如下： 1 unsigned char c = -1; // 假设char占8比特，c的值为255 当我们赋给无符号类型一个超出它表示范围的值时，结果是初始值对无符号类型表示数值总数取模后的余数。例如，8比特大小的unsigned char可以表示0至255区间内的值，如果我们赋了一个区间以外的值，则实际的结果是该值对256取模后所得的余数。因此，把-1赋给8比特大小的unsigned char所得的结果是255。 1 signed char c2 = 256; // 假设char占8比特，c2的值是未定义的 当我们赋给带符号类型一个超出它表示范围的值时，结果是未定义的（undefined）。此时，程序可能继续工作、可能崩溃，也可能生成垃圾数据。 1 2 3 4 unsigned u = 10; int i = -42; std::cout << i + i << std::endl; // 输出-84 std::cout << u + i << std::endl; // 如果int占32位，输出4294967264 在第一个输出表达式里，两个（负）整数相加并得到了期望的结果。在第二个输出表达式里，相加前首先把整数-42转换成无符号数。把负数转换成无符号数类似于直接给无符号数赋一个负数，结果等于这个负数加上无符号数的模。unsigned (int)的取值范围是0~\(2^{32}-1\)，所以总数有\(2^{32}\)个数，-42%\(2^{32}\)=-42+\(2^{32}\)，u+i=10+(-42+\(2^{32}\))=4294967264。 1 2 3 unsigned u1 = 42, u2 = 10; std::cout << u1 – u2 << std::endl; // 正确：输出32 std::cout << u2 – u1 << std::endl; // 正确：不过，结果是取模后的值 当从无符号数中减去一个值时，不管这个值是不是无符号数，我们都必须确保结果不能是一个负值。 ...

随机矩阵是这样一类方阵，其元素为非负实数，且行和或列和为1。如果行和为1，则称为行随机矩阵；如果列和为1，则称为列随机矩阵；如果行和和列和都为1，则称为双随机矩阵。 前面我们介绍的谷歌矩阵和HMM中的转移矩阵都属于随机矩阵，所以随机矩阵也称为概率矩阵、转移矩阵、或马尔可夫矩阵。 随机矩阵有一个性质，就是其所有特征值的绝对值小于等于1，且其最大特征值为1。下面通过两种方法证明这个结论。 首先，随机矩阵A肯定有特征值1，即 $$\begin{equation}A\vec 1=1\times\vec 1\end{equation}$$其中的单位向量\(\vec 1=(\frac{1}{n},…,\frac{1}{n})^T\)，因为A的行和为1，所以上述等式成立。即1是A的特征值。 反证法 假设存在大于1的特征值\(\lambda\)，则有\(A\vec x=\lambda\vec x\)。令\(x_k\)是\(\vec x\)中最大的元素。又因为A的元素非负，且行和为1，所以\(\lambda\vec x\)中的每个元素都是\(\vec x\)中元素的凸组合，所以\(\lambda\vec x\)中的每个元素都小于等于\(x_k\)。 $$\begin{equation}a_{i1}x_1+a_{i2}x_2+…+a_{in}x_n=\lambda x_i\leq x_k\end{equation}$$但是如果\(\lambda>1\)，则\(\lambda x_k>x_k\)，和(2)式矛盾，所以\(\lambda\leq 1\)。又因为(1)式，所以A的最大特征值为1。 常规证法 设对称随机矩阵A的特征值\(\lambda\)对应的特征向量为\(x\)（为了简便，以下省略向量符号），则有\(Ax=\lambda x\)，即\(x^TAx=\lambda x^Tx\)，欲证明\(|\lambda|\leq 1\)，只需证明 $$\begin{equation}\lambda=\frac{< x, Ax >}{< x, x >}\leq 1\end{equation}$$根据定义有： $$\begin{equation}< x, Ax >=\sum_{i=1}^na_{ii}x_i^2+2\sum_{i < j, i\sim j}a_{ij}x_ix_j\end{equation}$$对于\(i < j, i\sim j\)，有： $$\begin{equation}a_{ij}(x_i-x_j)^2=a_{ij}x_i^2-2a_{ij}x_ix_j+a_{ij}x_j^2\end{equation}$$两边求和并移项得到： $$ \begin{equation} \begin{array} \displaystyle{2\sum_{i < j}}a_{ij}x_ix_j & = & \displaystyle{\sum_{i < j}a_{ij}x_i^2+\sum_{i < j}a_{ij}x_j^2-\sum_{i < j}a_{ij}(x_i-x_j)^2}\\ & = & \displaystyle{\sum_{i < j}a_{ij}x_i^2+\sum_{i < j}a_{ji}x_j^2-\sum_{i < j}a_{ij}(x_i-x_j)^2}\\ & = & \displaystyle{\sum_{i < j}a_{ij}x_i^2+\sum_{i > j}a_{ij}x_i^2-\sum_{i < j}a_{ij}(x_i-x_j)^2}\\ & = & \displaystyle{\sum_i(\sum_{j\neq i}a_{ij}x_i^2)-\sum_{i < j}a_{ij}(x_i-x_j)^2}\\ & = & \displaystyle{\sum_i(x_i^2(1-a_{ii}))-\sum_{i < j}a_{ij}(x_i-x_j)^2} \end{array} \end{equation} $$第2、3个等号都是因为A是对称矩阵，所以可以把\(a_{ij}\)替换为\(a_{ji}\)，然后互换\(i,j\)下标。最后一个等号是因为A的行和为1。 ...

马尔可夫聚类算法（The Markov Cluster Algorithm, MCL）是一种快速可扩展的基于图的聚类算法。它的基本思想为：在一个稀疏图G中，如果某个区域A是稠密的（是一个聚类），则在A中随机游走k步，还在A内的概率很大，也就是说，A内的k步路径（k-length path）很多。所以我们可以在图中随机游走k步，如果某个区域连通的概率很大，则该区域是一个聚类。随机游走的下一步只和当前所处节点有关，也就是说这是一个马尔可夫的随机游走过程。 我们用一个例子来演示马尔可夫聚类算法的过程。 上图是一个很小的网络，我们用肉眼大概能看出有三个聚类，分别是左边的{1,6,7,10}，中间的{2,3,5}和右边的{4,8,9,11,12}。我们用MCL看看结果如何。 为了随机游走，我们常用邻接矩阵来表示图，如果i,j有边，则N[i][j]=1，否则N[i][j]=0。又随机游走可能有自回路，所以加上单位矩阵I，得到矩阵N+I。 MCL有两个关键的步骤，分别是Expansion和Inflation。 Expansion就是不断对矩阵进行幂次运算，相当于随机游走。假设随机游走了2步，则得到如下图的关联矩阵\((N+I)^2\)，第1行第10列为4，说明1到10的2-length path有4条：1→6→10，1→7→10，1→1→10，1→10→10。随机游走k步之后，\((N+I)^k[i][j]\)越大，说明\(i\)和\(j\)之间的连通性越强。 $$\begin{equation}Expand(M)=M^k\end{equation}$$ Inflation是为了增强更强的连接，减弱更弱的连接，只有这样才能得到边界比较明确的聚类。MCL的做法是对元素做幂次运算，然后按列归一化，公式为： $$\begin{equation}(\Gamma_rM)_{pq}=\frac{(M_{pq})^r}{\sum_{i=1}^k(M_{iq})^r}\end{equation}$$参数经验值是\(k=r=2\)。不断做Expansion和Inflation操作，直到算法收敛，得到若干个聚类。中间过程请点此查看，下图为最终结果。 从图中可以看出，和1有边的只剩下6,7,10了，所以得到聚类{1,6,7,10}，同理能得到聚类{2,3,5}和{4,8,9,11,12} ，和我们肉眼得到的结果是一致的。 MCL算法的原理很简单，得到的聚类效果也不错。下面总结一下MCL的算法过程： 给定无向图G，Expansion和Inflation的参数\(k\)和\(r\) 生成G的邻接矩阵\(N\) 添加自回路，得到矩阵\(N+I\) 循环对\(N+I\)做Expansion和Inflation操作，即计算公式(1)和(2)，直到收敛 根据最终得到的矩阵，进行划分聚类 此算法是我在上《生物信息学中的算法设计》课上是学到的，当时觉得这个算法真是神奇，如此简单，但又如此有效，实在高明。查阅文献得知，此为Stijn van Dongen的博士论文，本博客的图片均来自其博士论文，想深入了解图聚类算法，请下载他的论文。

上一回介绍了HMM的解码问题，今天我们介绍HMM的学习问题和识别问题，先来看学习问题。 正如上一回结束时所说，HMM的学习问题是：仅已知观测序列\(\vec y\)，要估计出模型参数组\(\vec\lambda=(\mu,A,B)\)，其中\(\mu\)为初始概率分布向量，\(A\)为转移概率矩阵，\(B\)为发射概率矩阵。 算法设计 求解HMM的参数学习问题，就是求解如下的最优化问题： $$\begin{equation} P(\vec Y = \vec y|\hat \lambda)=\max\limits_{\vec \lambda} P(\vec Y = \vec y|\vec \lambda)\end{equation}$$也就是找一个参数\(\vec \lambda\)，使得模型在该参数下最有可能产生当前的观测\(\vec y\)。如果使用极大似然法求解，对于似然函数\(P(\vec Y=\vec y|\vec \lambda)=\sum\limits_{i_1,…,i_T}\mu_{i_1}b_{i_1y_1}a_{i_1i_2}…a_{i_{T-1}i_T}b_{i_Ty_T}\)而言，这个最大值问题的计算量过大，在实际中是不可能被采用的。为此，人们构造了一个递推算法，使其能相当合理地给出模型参数\(\vec \lambda\)的粗略估计。其核心思想是：并不要求备选\(\vec\lambda\)使得\(P(\vec Y=\vec y|\vec \lambda)\)达到最大或局部极大，而只要求使\(P(\vec Y=\vec y|\vec \lambda)\)相当大，从而使计算变为实际可能。 EM算法 为此，我们定义一个描述模型“趋势”的量\(Q(\vec\lambda^*|\vec\lambda)\)代替似然函数\(P(\vec Y=\vec y|\vec\lambda)\)，其定义为： $$\begin{equation} Q(\vec\lambda^*|\vec\lambda)=\sum\limits_{\vec x}P(\vec x,\vec y|\vec\lambda)\ln P(\vec x,\vec y|\vec\lambda^*)\end{equation}$$利用在\(0 < x < 1\)时，不等式\(\ln x\leq x-1\)成立，可以证明： $$\begin{equation} Q(\vec\lambda^*|\vec\lambda)-Q(\vec\lambda|\vec\lambda)\leq P(\vec Y=\vec y|\vec\lambda^*)-P(\vec Y=\vec y|\vec\lambda)\end{equation}$$由此可见，对于固定的\(\vec\lambda\)，只要\(Q(\vec\lambda^*|\vec\lambda)>Q(\vec\lambda|\vec\lambda)\)，就有\(P(\vec Y=\vec y|\vec\lambda^*)>P(\vec Y=\vec y|\vec\lambda)\)。于是想把模型\(\vec\lambda_m\)修改为更好的模型\(\vec\lambda_{m+1}\)，只需找\(\vec\lambda_{m+1}\)使得： $$\begin{equation}Q(\vec\lambda_{m+1}|\vec\lambda_m)=\sup_{\vec\lambda}Q(\vec\lambda|\vec\lambda_m)\end{equation}$$即只要把\(Q(\vec\lambda|\vec\lambda_m)\)关于\(\vec\lambda\)的最大值处取成\(\vec\lambda_{m+1}\)，就有\(P(\vec Y=\vec y|\vec\lambda_{m+1})>P(\vec Y=\vec y|\vec\lambda_m)\)。 这样得到的模型序列\(\{\vec\lambda_m\}\)能保证\(P(\vec Y=\vec y|\vec\lambda_m)\)关于\(m\)是严格递增的，虽然在这里还不能在理论上证明\(P(\vec Y=\vec y|\vec\lambda_m)\)收敛到\(\max_{\vec\lambda}P(\vec Y=\vec y|\vec\lambda)\)，但是当\(m\)充分大时，\(\vec\lambda_m\)也还能提供在实际中较为满意的粗略近似。 ...

隐马尔可夫模型（Hidden Markov Model, HMM）是统计模型，它用来描述一个含有隐含未知参数的马尔可夫过程。其难点是从可观察的参数中确定该过程的隐含参数。然后利用这些参数来作进一步的分析，例如模式识别。 先举一个简单的例子以直观地理解HMM的实质——韦小宝的骰子。 假设韦小宝有两个骰子，一个正常的骰子A，A以1/6的概率均等的出现每个点；一个不正常的骰子B，B出现5,6点数的概率为0.3，出现其他点数的概率为0.1。显然投掷B更容易出现大的点数。每次试验第一次投掷时，韦小宝会以0.4的概率出千（即投掷B）。但是在一次试验中，韦小宝不太可能一直出千，所以骰子会在A、B之间转换，比如这次投了B，下次可能会以0.1的概率投A。A、B之间的转移概率如下图。 某一次试验，我们观察到韦小宝掷出的骰子序列为\(O=(1,3,4,5,5,6,6,3,2,6)\)，请问韦小宝什么时候出千了。这个问题就可以通过HMM求解。 HMM有2个状态： 观测状态。我们观察到的骰子序列称为观测状态\(\mathbf{Y}=\{y_1,y_2,…,y_T\}\) 隐状态。隐含在每个观测状态里面的是隐状态\(\mathbf{X}=\{x_1,x_2,…,x_T\}\) T是时间，也可以认为是观测的次数。HMM有3个参数： 初始分布\(\mathbf{\mu}=(\mu_i)\)，\(\mu_i=Pr(x_1=i)\)，即第一次观测时，每个隐状态出现的概率 转移概率矩阵\(A=(a_{ij})\)，\(a_{ij}=Pr(x_{t+1}=j|x_t=i)\)，即t时刻的隐状态为i，t+1时刻转移到隐状态j的概率 发射概率矩阵\(B=(b_{il})\)，\(b_{il}=Pr(y_t=l|x_t=i)\)，即t时候隐状态为i的情况下，观测到状态为l的概率 参数\(\mathbf{\lambda=\{\mu,A,B\}}\)称为HMM的模型参数。具体到上面的例子，我们有初始分布和转移概率为： 发射概率为： 观测状态为\(\mathbf{Y}=(1,3,4,5,5,6,6,3,2,6)\)，问题就是求解出隐状态\(\mathbf{X}\)，此问题被称为HMM的解码问题，可以由著名的维特比算法（Viterbi algorithm）解决。 解码问题是要求出使得观测状态\(Y\)出现概率最大的隐状态\(X\)，假设有N个隐状态（本例为2），共有T个时刻（本例为10），则每个时刻有N个取值可能，则共有\(N^T\)条可能的隐状态链（本例为\(2^{10}\)）。我们需要求出每一条隐状态链下T个发射概率的乘积，然后取最大值，这是指数时间复杂度的（\(O(N^T)\)）。 但是Viterbi算法是一个动态规划算法，只需多项式时间即可解决该问题。该算法的原理很好理解，假设我们求得到\(s_{i2}\)的最大概率路径为下图中的红线\(s_{11}\rightarrow s_{22}\rightarrow … s_{i2}\)，则在求经过\(s_{i2}\)到\(s_{(i+1)1}\)的最大概率路径时，不需要再测试\(s_{13}\rightarrow s_{21}\rightarrow s_{i2}\rightarrow s_{(i+1)1}\)这条路径（下图蓝线），因为显然已经知道红线概率大于蓝线概率了。图中还有很多类似蓝线的路径都可以不用计算了，大大提高了求解速度。 因为计算第\(i+1\)时刻的累积概率只和第\(i\)时刻的概率有关，每次至多计算\(N*N\)个概率乘积（可以从\(i\)时刻的\(N\)个状态到达\(i+1\)时刻的某个状态，\(i+1\)时刻共有\(N\)个状态），最多计算T次（共T个时刻），所以时间复杂度降到了\(O(N^2T)\)。 下面我们形式化的描述Viterbi算法。 假设\(\delta_t(i)\)为\(t\)时刻取到隐状态\(i\)，且1~t的观测状态都符合观测值\(Y\)的各个路径的最大概率，即 $$ \begin{equation}\delta_t(i)=\underset{i_1,…,i_{t-1}}{\max}Pr(X_t=i,X_{t-1}=i_{t-1},…,X_1=i_1,Y_t=y_t,…,Y_1=y_1|\mathbf{\lambda})\end{equation} $$联系上图，可认为\(\delta_t(i)\)为红线。则递推公式为： $$ \begin{equation}\delta_{t+1}(i)=b_{iy_{t+1}}\underset{j}{\max}(\delta_t(j)a_{ji})\end{equation} $$由\(j\)到\(i\)的转移概率，再乘上\(i\)发射\(y_{t+1}\)的概率。 在初始时刻\(t=1\)，有： $$ \begin{equation}\delta_1(i)=\mu_ib_{iy_1}\end{equation} $$最后的全局最大概率为\(\underset{j}{\max}\delta_T(j)\)。为了得到完整路径，我们保留每一隐状态取得最大概率时的上一隐状态，即： $$ \begin{equation}\psi_{t+1}(i)=j^*\end{equation} $$其中\(j^*\)要满足 $$ \begin{equation}\delta_{t+1}(i)=b_{iy_{t+1}}\delta_t(j^*)a_{j^*i}\end{equation} $$最后使用如下回溯法得到所有最佳隐状态： $$\begin{equation}X_T=i^*\in\{i:\delta_T(i)=\underset{j}{\max}\delta_T(j)\}\end{equation}$$$$\begin{equation}X_t=\psi_{t+1}(X_{t+1})\end{equation}$$下面我们利用Viterbi算法来求解韦小宝的骰子这个例子。 \(t=1\)时，\(y_1=1\)，有\(\delta_1(A)=0.6*1/6=0.1\)，\(\delta_1(B)=0.4*0.1=0.04\)。 \(t=2\)时，\(y_2=3\)，有： 隐状态为A：a）A->A有\(\delta_2(A)=(1/6)*0.1*0.8=1.33*10^{-2}\)；b）B->A有\(\delta_2(A)=(1/6)*0.04*0.1=6.6*10^{-4}\)。所以A->A，\(\psi_2(A)=A\)。 隐状态为B：a）A->B有\(\delta_2(B)=0.1*0.1*0.2=2*10^{-3}\)；b）B->B有\(\delta_2(B)=0.1*0.04*0.9=3.6*10^{-3}\)。所以B->B，\(\psi_2(B)=B\)。 如此计算下去，可以得到如下表： \(t=10\)时最大概率为\(\delta_{10}(B)\)，经过回溯得到最佳隐状态为： 所以HMM很神奇吧，可以抓住韦小宝从第5次开始就一直在出千，而且出千之后，掷出的点数大部分为5和6。 Viterbi算法还可用于解决语音识别或者拼音输入法。我们知道中文的一个拼音可以对应多个汉字，连续的一段拼音就能组成成千上万种可能的句子，哪一个句子才是最佳候选呢？我们可以把每个拼音当成观测状态，同音的汉字当成可能的隐状态。通过背景语料库统计得到每个汉字出现在词首的概率、汉字之间的转移概率和汉字与拼音之间的发射概率，这样我们就能得到模型参数，然后利用Viterbi算法求解出一个最佳的隐状态序列，这样就能完成一个简易的拼音输入法。 HMM在实际中主要有3个方面的应用，分别是： 从一段观测序列\(\mathbf{Y}\)及已知模型\(\mathbf{\lambda=(\mu,A,B)}\)出发，估计出隐状态\(\mathbf{X}\)的最佳值，称为解码问题，这是状态估计问题。这篇博客讨论的就是这个问题。 从一段观测序列\(\mathbf{Y}\)出发，估计模型参数组\(\mathbf{\lambda=(\mu,A,B)}\)，称为学习问题，就是参数估计问题。 对于一个特定的观测链\(\mathbf{Y}\)，已知它可能是由已经学习好的若干模型之一所得的观测，要决定此观测究竟是得自其中哪一个模型，这称为识别问题，就是分类问题。 关于HMM的学习问题和识别问题，请听下回分解。

我们知道常规的快速排序算法是一个不稳定的算法，也就是两个相等的数排序之后的顺序可能和在原序列中的顺序不同。这是因为当选定一个枢轴（pivot），要把其他数分到小于pivot和大于pivot的两边的时候，不同实现的分法不一样。 下面我实现了一种稳定版快速排序算法，在Partition函数中保持了原序列中所有元素的相对顺序，只把pivot放到了它的正确位置。具体方法是三遍扫描原序列：1）第一遍先把小于pivot的元素按先后顺序放到tmp里，然后把pivot放到它的正确位置tmp[k]；2）第二遍把大于pivot的元素按先后顺序追加在tmp里，这样除了pivot以前的其他元素，都保持了和原序列中一样的顺序；3）第三遍把tmp赋值回原数组A。 当排序算法稳定之后，就可以借此统计逆序数了，文件Q5.txt中共包含100000个不同的整数，每行一个数。我们可以使用稳定版快速排序算法对其排序，并统计出其中的逆序数个数。 具体的Python 3实现如下： 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 # -*- coding: utf-8 -*- """ Created on Tue Oct 6 00:21:37 2015 @author: bitjoy """ import time inversions = 0 def Partition(A, p, r): global inversions tmp = [0] * (r-p+1) pivot = A[p] k = 0 for i in range(p+1, r+1): # first if A[i] < pivot: tmp[k] = A[i] inversions = inversions + i – k – p k = k + 1 tmp[k] = pivot ans = k + p k = k + 1 for i in range(p+1, r+1): # second if A[i] > pivot: tmp[k] = A[i] k = k + 1 k = 0 for i in range(p, r+1): # third A[i] = tmp[k] k = k + 1 return ans def QuickSortAndCount(A, p, r): if p < r: q = Partition(A, p, r) QuickSortAndCount(A, p, q-1) QuickSortAndCount(A, q + 1, r) if __name__ == "__main__": Q5 = open('Q5.txt', encoding = 'utf-8') data = [ int(x) for x in Q5 ] Q5.close() start = time.clock() QuickSortAndCount(data, 0, len(data) -1 ) end = time.clock() print("number of inversions:%d\ntime:%f s"%(inversions,end-start)) 虽然这种快排的时间复杂度还是O(nlgn)，但是在Partition函数中扫描了3次数组，并且借用了辅助数组tmp，不再是in-place排序算法，所以排序用时会比常规快排或者归并排序要慢。 ...