0和1

上一回介绍了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排序算法，所以排序用时会比常规快排或者归并排序要慢。 ...

哈弗曼编码是一个很经典的压缩算法，压缩率能达到50%，甚至更低。它的基本原理包括四个步骤： 统计文件中每个字符出现的频率。 构建一个哈弗曼树。建树的过程是不断的合并频率最小的两个节点，父亲节点的频率为两个孩子节点的频率之和。如此循环直到合并成一个根节点。叶子节点为不同的字符及其频率。 生成哈弗曼编码。从树根开始对树进行编码，比如进入左孩子的边标记为0，进入右孩子的边标记为1，这里的0和1都是二进制位。这样之后，每个叶子节点都有一个唯一的二进制编码，这就是哈弗曼编码。频率越低的字符哈弗曼编码越长，频率越高的字符哈弗曼编码越短，这样就能起到压缩的效果。 第二遍扫描文件，把字符转换为对应的哈弗曼编码，保存成压缩文件。 解压缩的过程就是解析二进制位，然后查找哈弗曼树，每找到一个叶子节点，就解析出一个字符，直到解析完所有二进制位。下面详细解释我的C++实现。 首先定义一个哈弗曼编码类，对外只提供压缩Compress和解压缩Decompress两个接口。值得注意的是有一个Node结构体，用于构成哈弗曼树的节点。此外count_node的key是字符频率，value是所在节点，且是multimap类型的，所以count_node会自动按字符频率有小到大排序，在构建哈弗曼树时，每次只需要取count_node的前两个节点进行合并即可。 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 class HuffmanCode { public: HuffmanCode(); void Compress(string src, string dest); void Decompress(string src, string dest); virtual ~HuffmanCode(); private: void CountLetter(string src); void ConstructHuffmanTree(); void GenerateHuffmanCode(); void WriteHuffmanCode(ofstream &os); void Compressing(string src, string dest); void InsertIntoHuffmanTree(char letter, string &code, int &k); void ConstructHuffmanTreeFromFile(ifstream &is); void Decompressing(ifstream &is, ofstream &os); map<char, int> letter_count; typedef struct Node { int id; bool is_leaf; char letter; int parent, lchild, rchild; Node() { } Node(int i, bool il, char lt, int p, int lc, int rc) : id(i), is_leaf(il), letter(lt), parent(p), lchild(lc), rchild(rc) { } }; multimap<int, Node> count_node; vector<Node> huffman_tree; map<char, vector<char>> letter_hcode; // hufman code for each letter }; 压缩函数Compress串起压缩的整个流程，包括统计字符频率、构建哈弗曼树、生成哈弗曼编码以及最后将原始文件转换成哈弗曼编码的二进制文件。 ...

之前写了七篇博客详细介绍了搜索引擎的工作原理。彼时的搜索引擎主要讲查询和网页的相关性匹配，是动态的、在线的、实时的。相关性匹配有一个问题，网页很容易作弊，比如可以在一个网页中写满诸如“免费”、“美容”之类的垃圾关键词，进而提升查询相关性。但是用户在查询时，一定希望返回的网页比较权威可信，比如同样搜索“苹果电脑”，排名第一的应该是Apple的官网，而不应该是中关村在线之类的第三方网站。 权威性是一个静态的（或者说变化较慢的）衡量网页重要性的指标。但是应该怎样度量权威性呢，HITS算法使用authority来度量，即指向自身的网页数量越多，则自身的authority值越大。谷歌的PageRank算法是用PageRank值来衡量权威性的。HITS和PageRank一个比较大的区别是HITS和查询有关，而PageRank和查询无关，所以PageRank可以离线计算。下面主要介绍PageRank算法。 PageRank’s thesis is that a webpage is important if it is pointed to by other important pages. 我先不加解释的给出PageRank的公式，然后带领大家一步步推导出这个公式。 $$\pi^T=\pi^T(\alpha S+(1-\alpha)E)$$我们首先明确目标：PageRank计算的是网页的静态权威度（PR值），也就是如果给定了一个网络结构，则每个网页的PR值就可以通过PageRank算法计算出。假设网页\(P_i\)的PR值为\(r(P_i)\)，则\(r(P_i)\)等于所有指向\(P_i\)的网页的PR值之和，即 $$\begin{equation}r(P_i)=\sum\limits_{P_j\in B_{P_i}}\frac{r(P_j)}{|P_j|}\end{equation}$$其中\(B_{P_i}\)为指向\(P_i\)的网页集合，\(|P_j|\)为\(P_j\)的出边的数量。这个式子很好理解，包括两方面内容：1）\(\sum\limits_{P_j\in B_{P_i}}\)表示如果指向\(P_i\)的网页数量越多，说明网页\(P_i\)越重要；2）\(\frac{r(P_j)}{|P_j|}\)表示如果\(P_j\)指向的页面数量越少，但有一个指向了\(P_i\)，说明网页\(P_i\)越重要（如果一个大牛写了很多推荐信（\(|P_j|\)大），则这些推荐信的效力就下降了，如果大牛只给你写了推荐信（\(|P_j|=1\)），则这封推荐信的效力一定很高）。 (1)式有一个问题，初始给定一个网络结构时，并不知道\(r(P_i), r(P_j)\)，如何计算呢？Brin和Page利用递归的思想求解，初始假设所有网页的PR值相等，都为\(\frac{1}{n}\)，其中\(n\)为网络中网页的数量。则第\(k+1\)轮的PR计算公式为： $$\begin{equation}r_{k+1}(P_i)=\sum\limits_{P_j\in B_{P_i}}\frac{r_k(P_j)}{|P_j|}\end{equation}$$初始对所有网页\(P_i\)有\(r_0(P_i)=\frac{1}{n}\)，迭代\(k\)步之后，可以计算出所有网页的PR值，然后按PR值从大到小排序，就可以知道每个网页的重要性了。 对于上图的小网络，我们可以计算出其每一步的PR值： 可以看到经过2次迭代之后，节点4的PR值最大，从图中也可以看出，节点4的出入边较多，它可能比较重要。 注意到对于(2)式，当\(i,j\)之间有边时，\(\frac{1}{|P_j|}\)相当于对\(P_j\)出度的归一化，设矩阵\(H\)为图的邻接矩阵的行归一化矩阵，对于上图，为 设行向量\(\pi^{(k)T}\)为第\(k\)轮迭代时所有网页的PR值，则式(2)可以转换为如下的矩阵形式： $$\begin{equation}\pi^{(k+1)T}=\pi^{(k)T}H\end{equation}$$初始有\(\pi^{(0)T}=\frac{1}{n}e^T\)，\(e^T\)为全1的行向量。我们可以从(3)式观测出几点信息： (3)式的每一轮计算涉及到向量和矩阵的乘法，复杂度为\(O(n^2)\)，\(n\)为矩阵\(H\)的大小 \(H\)是一个稀疏矩阵，因为大部分网页只和很少的网页有链接关系，所以上述向量和矩阵的乘法复杂度还可以降低 \(H\)有点像马尔科夫链中的随机转移矩阵，但又不完全是，因为如果有dangling nodes，则这一行就是全0，所以\(H\)被称为substochastic matrix 上图中的节点3就是一个dangling node，它只有入边，没有出边，也就是说，每一轮迭代，PR值只会流入3号节点，不会从3号节点流出，久而久之，3就像一个水槽(sink)一样，吸走了大部分的PR，导致PR值虚高。 所以问题随之而来，怎样保证(3)式一定能够收敛到一个平稳概率分布\(\pi^T\)，\(\pi^T\)和\(\pi^{(0)T}\)有关吗，怎样解决dangling nodes问题，等等。此时需要引入一点马尔科夫链理论的知识。 在马尔科夫理论呢中，如果一个矩阵\(P\)是随机的（stochastic）、不可约的（irreducible）和非周期的（aperiodic），则对于任意的起始向量，都能收敛到一个唯一的平稳正向量。所以如果PageRank矩阵\(H\)满足上述三个条件，则可以用幂法（Power Method）找到一个平稳概率分布\(\pi^T\)。幂法是用来计算最大特征值的特征向量。因为\(H\)的最大特征值为1，所以可以用幂法找到稳态时（\(\pi^T=\pi^TH\)）的概率分布\(\pi^T\)。 下面我们就将矩阵\(H\)调整为随机的（stochastic）、不可约的（irreducible）和非周期的（aperiodic）。 行随机矩阵是指行和为1的非负矩阵。如果图中含有dangling nodes，则\(H\)不是随机的，比如上面的例子，第二行为全0。所以第一个调整是对于所有dangling nodes，都加上一个随机跳转向量\(e^T/n\)，含义就是如果进入死胡同（dangling nodes），则随机跳转到网络中的任意一个网页。定义向量\(a\)： $$\begin{equation}a_i=\begin{cases}1\quad\text{if page}~i\text{ is a dangling node}\\0\quad\text{otherwise}\end{cases}\end{equation}$$则新的Google矩阵为： $$\begin{equation}S=H+a\frac{1}{n}e^T\end{equation}$$新矩阵\(S\)就是一个行随机矩阵了。对于上图的例子，有 为了保证矩阵\(S\)满足不可约性（irreducible）和非周期性（aperiodic），必须使\(S\)对应的图是强连通的且每个节点有自回路。所以再次调整为： $$\begin{equation}G=\alpha S+(1-\alpha)\frac{1}{n}ee^T\end{equation}$$令 $$\begin{equation}E=\frac{1}{n}ee^T\end{equation}$$则得到本博客开头的Google矩阵公式： $$\begin{equation}G=\alpha S+(1-\alpha)E\end{equation}$$\(E\)即为随机平均游走矩阵。矩阵\(G\)也很好解释，大家上网的时候以\(\alpha\)的概率沿着某个网页里面的链接一步步深入进去（\(S\)），当沿着链接走累的时候，以\(1-\alpha\)的概率在地址栏输入一个新地址，随机跳走了（\(E\)）。 此时的矩阵\(G\)满足随机性（stochastic）、不可约性（irreducible）和非周期性（aperiodic），所以可以根据幂法（Power Method）找到一个平稳概率分布\(\pi^T\)，\(\pi^T_i\)就衡量了网页\(P_i\)的重要性或者权威性。 ...

$$\begin{equation}Pr(|\hat{p}-p|\geq 5\%)\leq 5\%\end{equation}$$上一回我们讲到当\(p\)本身很小的时候，容易被5%（绝对误差）给淹没掉，导致结果的不可信。我们可以引入相对误差，把(1)式转换为如下的不等式 $$\begin{equation}Pr(|\hat{p}-p|\geq\delta p)\leq\epsilon\end{equation}$$同理，我们可以用 $$\begin{equation}\hat{p}=\frac{x_1+x_2+…+x_n}{n}\end{equation}$$代替\(\hat{p}\)（建议先看上一篇博客），转换为 $$\begin{equation}Pr(|X-np|\geq\delta np)\end{equation}$$类似的，\(X=x_1+x_2+…+x_n\)，\(E(X)=\mu=np\)，所以(4)式等价为 $$\begin{equation}Pr(|X-\mu|\geq\delta\mu)\end{equation}$$这个时候，因为不等号右边和均值\(\mu\)有关，不能再用切比雪夫不等式了，我们需要另外一个武器：Chernoff bound。它有两种形式： $$\begin{equation}Pr(X\geq (1+\delta)\mu)\leq[\frac{e^\delta}{(1+\delta)^{1+\delta}}]^\mu\leq e^{-\frac{\mu}{3}\delta^2}\quad\forall\delta>0\end{equation}$$$$\begin{equation}Pr(X\leq (1-\delta)\mu)\leq[\frac{e^{-\delta}}{(1-\delta)^{1-\delta}}]^\mu\leq e^{-\frac{\mu}{2}\delta^2}\quad\forall 0<\delta<1\end{equation}$$Chernoff bound的证明需要用到马尔可夫不等式，有一点技巧。以上两种形式可以统一成 $$\begin{equation}Pr(|X-\mu|\geq\delta\mu)\leq 2e^{-\frac{\mu}{3}\delta^2}\end{equation}$$也是一个很漂亮的不等式。 利用Chernoff bound求解(5)式： $$\begin{equation}Pr(|X-\mu|\geq\delta\mu)\leq 2e^{-\frac{\mu}{3}\delta^2}\\=2e^{-\frac{np}{3}\delta^2}\leq\epsilon\end{equation}$$解得 $$\begin{equation}n\geq\left\lceil\frac{3ln\frac{2}{\epsilon}}{p\delta^2}\right\rceil\end{equation}$$这个结果看起来就很复杂了。也就是说，如果要设计调查问卷使满足(2)式的精度，抽样的样本数必须满足(10)式。从(10)式可知，当要求的精度越高（即\(\delta\)和\(\epsilon\)越小），所需的样本数越大。并且结果还和真实值\(p\)有关。

每年春晚过后，央视又要吹嘘说今年春晚收视率创新高了，但是我们总感觉央视在骗我们，因为我是越长大越不看春晚了[笑cry]，所以收视率到底是怎么统计出来的，央视的说法是否靠谱呢？ 最近的美国大选真是热闹，很多机构都会发放一些调查问卷，然后统计出希拉里或者唐纳德的民众支持率是多少，但是我并没有收到调查问卷，凭什么就得出了民众支持率了，意思是把我排除在民众之外咯？所以引出这样一个问题，调查问卷是否可信，即调查问卷的有效性。 其实，央视统计收视率并不要问全中国14亿人口有多少人看了春晚，他只需要从14亿人口里面随机抽\(n\)个人，问一下这\(n\)个人里有多少人看了春晚，然后把看的人数除以总数就大概估计出全国的收视率了。同理调查民众支持率也是一样，只需要随机调查\(n\)个人的意向，把支持希拉里的人数除以总数就大概得到了希拉里的支持率。 但是你要问了，通过抽样调查出来的收视率和支持率靠谱吗，需要随机抽样多少人才能得到一个比较好的全局近似解呢？今天我们就来解决这个问题。 假设我们随机抽样了\(n\)个人，分别是\(x_1,x_2,…,x_n\)。如果第\(i\)个人看了春晚，则\(x_i=1\)，否则\(x_i=0\)。那么通过这\(n\)个人的收视情况，我们可以估计出一个收视率 $$\begin{equation}\hat{p}=\frac{x_1+x_2+…+x_n}{n}\end{equation}$$假设全国的真实收视率是\(p\)，那么平均到每一个人，他看了春晚的概率就是\(p\)，也即\(Pr(x_i=1)=p\)，所以有 $$\begin{equation}E(x_i)=p\quad E(x_i^2)=p\quad Var(x_i)=p(1-p)\end{equation}$$我们的目的就是希望通过\(n\)个人估计出来的\(\hat{p}\)和\(p\)越接近越好。换句话说，我们希望\(\hat{p}\)和\(p\)相差大于5%的概率要小于5%。再换句话说就是有至少95%的概率，\(\hat{p}\)和\(p\)相差在5%以内，即\(\hat{p}\)和\(p\)很接近。注意这里的两个5%都是可以换成任意你想要的精度。用数学语言表示就是，\(n\)至少为多少时，以下不等式可以被满足。 $$\begin{equation}Pr(|\hat{p}-p|\geq 5\%)\leq 5\%\end{equation}$$把(1)式代入(3)式，用\(\frac{1}{20}\)代替5%，得到等价形式： $$\begin{equation}Pr(|(\frac{x_1+x_2+…+x_n}{n})-p|\geq\frac{1}{20})\\ \Longleftrightarrow~Pr(|X-np|\geq\frac{n}{20})\end{equation}$$其中\(X=x_1+x_2+…+x_n\)。根据期望的线性可加性，有 $$\begin{equation}E(X)=E(x_1+x_2+…+x_n)=E(x_1)+E(x_2)+…+E(x_n)=np\end{equation}$$所以(4)又等价于 $$\begin{equation}Pr(|X-E(X)|\geq\frac{n}{20})\end{equation}$$我们需要利用著名的切比雪夫不等式来求解上式，切比雪夫不等式如下： $$\begin{equation}Pr(|X-E(X)|\geq~c)\leq\frac{Var(X)}{c^2}\end{equation}$$切比雪夫不等式可以直接由马尔可夫不等式得到，马尔可夫不等式的证明也不难，略过。 利用切比雪夫不等式求解(6)式 $$\begin{equation}Pr(|X-E(X)|\geq\frac{n}{20})\leq\frac{Var(X)}{n^2}*400\\ =\frac{n*Var(x_i)}{n^2}*400\\ =\frac{p(1-p)}{n}*400\\ \leq\frac{1/4}{n}*400=\frac{100}{n} \end{equation}$$第一个等号是因为\(n\)个变量是独立同分布的，所以方差也有类似于(5)式的线性性质。最后一个不等号是因为\(p(1-p)\)是一个开口向下的抛物线，在\(p=1/2\)时取到极值\(1/4\)。 回到最初的不等式(3)，则(8)式要满足\(\frac{100}{n}\leq 5\%\)，解得\(n\geq 2000\)。注意到求出的\(n\)和总体人数是无关的，也就是说，虽然全中国有十几亿人口，但是央视只要随机抽样调查2000个人的收视情况，就能以比较高的概率准确估计出全国的收视率。 这个结论还是很漂亮的，但是这种方法有两个限制条件： 采样满足独立同分布，即这\(n\)个人是独立同分布的，不能针对某一特定人群调查 (3)式的5%是一个绝对误差，当\(p\)本身很小的时候，容易被5%淹没 对于第1个问题，稍微好处理一点，抽样的时候尽量随机一点。对于第2个问题，比较好的解决办法是引入相对误差，即把(3)式转换为如下的不等式 $$\begin{equation}Pr(|\hat{p}-p|\geq\delta p)\leq\epsilon\end{equation}$$(9)式的求解就比较复杂了，得出的结论也没有上面那么简单，具体的求解方法请听下回分解。

你是否想过如下问题：怎样向色盲证明两只袜子的颜色是不一样的？怎样证明两个图是不同构的？怎样证明一个数是二次非剩余的？ 咋听起来觉得很有意思吧，色盲是区分不了颜色的，怎么能让他相信两只袜子的颜色不一样呢。图同构问题目前既没有被证明属于P，也没有被证明属于NP-Complete。二次非剩余问题也没有被证明属于NP。 这些听起来很“难”的问题，却可以通过交互式证明进行证明，下面先通过“向色盲证明两只袜子的颜色不同”这个有趣的例子一窥交互式证明的强大。 向色盲证明两只袜子的颜色不同 P有一只红袜子和黄袜子，她的一个色盲朋友V不相信P的袜子颜色不同，P如何才能让V相信这是真的呢？一个简单的办法如下： P把两只袜子给V，V每只手拿了一只袜子 P转过身背对V V抛一枚硬币，如果头面朝上，则保持两只袜子不动，否则交换左右手的袜子 P转过身，V问P是否交换过袜子 如果P回答错误，则V不相信；否则，重复100次实验，如果P都回答正确，则V相信这两只袜子是不同颜色的 如果两只袜子的颜色确实不一样，则P通过区分两只袜子的颜色能正确回答V有没有交换过袜子。但是如果两只袜子颜色一样，则不管V有没有交换过，P都无法分辨这两只袜子，所以只好猜V有没有交换，而猜对的概率只有1/2，重复100次，都猜对的概率只有\((1/2)^{100}\)，这是一个非常小的数，可以认为几乎不会发生，即出错的概率极低。 这就是交互式证明的一个例子，上述证明有三个特点：1）交互过程，整个证明需要P和V进行交互才能完成；2）具有随机性，即V需要抛一枚硬币，来决定是否交换袜子；3）零知识，虽然V最终相信了这两只袜子是不同颜色的，但V还是不知道这两只袜子是什么颜色的。 下面我们给出交互式证明的形式化定义。 交互式证明（Interactive Proofs, IP） 令\(k\)是\(N\rightarrow N\)的一个函数，我们称语言\(L\)属于\(IP[k]\)，如果存在一个\(k(|x|)\)多项式时间概率图灵机TM \(V\)，使得： $$ \begin{equation} \begin{cases} x\in L \Longrightarrow\exists P\quad Pr[V \text{ accepts }x, V(x)=1]\geq 2/3 \\ x\notin L \Longrightarrow\forall P\quad Pr[V \text{ accepts }x, V(x)=1]\leq 1/3 \end{cases} \end{equation} $$定义 $$IP=\underset{c}{\bigcup} IP[n^c]$$上述定义的第一条称为“完备性”（Completeness）：如果\(x\in L\)，则存在一个证明者P（prover），使得验证者V（verfier）能以多项式时间接受\(x\)，且接受的概率大于2/3；第二条称为“可靠性”（Soundness），如果\(x\notin L\)，则对于所有证明者P，V接受\(x\)的概率都不会超过1/3。 对应到上面的例子，完备性：当两只袜子的颜色确实不同时，V接受的概率为1>2/3；可靠性：当两只袜子的颜色相同时，重复100次实验，V接受的概率只有\((1/2)^{100}<1/3\)。 IP这个复杂性类就是所有IP[k]的并集。在IP中，P的能力是无穷的，但它不一定是诚实的；V能力较弱，只能进行多项式时间的计算。 下面我们给出另外两个交互式证明的协议。 图不同构（Graph Non-isomorphism, GNI）的交互式证明 如果图\(G_1\)和\(G_2\)可以通过对顶点进行恰当标记来将它们转换为同一个图，则称\(G_1\)和\(G_2\)同构，记为\(G_1\cong G_2\)。换句话说，如果\(G_1\)和\(G_2\)同构，则在\(G_1\)的所有顶点标签上存在一个置换\(\pi\)使得\(\pi (G_1)=G_2\)，其中\(\pi (G_1)\)是将\(\pi\)作用到\(G_1\)的各个顶点上之后得到的图。下图就是两个同构图，右边给出了置换\(\pi\)。 图同构的补集为图不同构（Graph Non-isomorphism, GNI），即判定给定的两个图是否不同构。下面是GNI的一个交互式证明过程。 给定两个图\(G_1\)和\(G_2\)，证明者P想要向验证者V证明\(G_1\ncong G_2\)。 V：随机选一个\(i\in \{1,2\}\)，对\(G_i\)做一个随机的置换，得到新图\(H\)，则有\(H\cong G_i\)，将\(H\)发送给P P：发送\(j\)给V V：如果\(i\neq j\)，则拒绝；否则重复100次实验，都有\(i==j\)，则相信\(G_1\ncong G_2\) 完备性：如果\(G_1\ncong G_2\)，则\(H\)只和\(G_1, G_2\)中的一个图同构，P因为能力无穷，一定能找出和\(H\)同构的图\(G_j\)，且满足\(j==i\)。 ...

\(\LaTeX\)的强大我就不赘述了，这里简单介绍一下怎样在Windows快速配置一个完美好用的\(\LaTeX\)环境。 我大学刚接触\(\LaTeX\)的时候，使用的是CTeX，CTeX是一个大礼包，整合了编译器编辑器等，但是由于久不更新，很多宏包和语法都变了，而且CTeX附带的WinEdt是商业软件，30天之后需要收费，我又不想用盗版，所以就打算自己配置\(\LaTeX\)环境。 目前使用的是MiKTeX+Texmaker的完美组合！MiKTeX是\(\LaTeX\)编译器，Texmaker是\(\LaTeX\)编辑器。两者都是开源软件。 MiKTeX非常棒的地方在于“MiKTeX has the ability to install needed packages automatically (on-the-fly)”，就是说，你用MiKTeX时，不需要担心某个宏包是否存在，你只管用就是了，MiKTeX会在你第一次用到某个宏包时，自动从网上下载，非常方便。正因为这样，MiKTeX的安装包很小，只有175MB。当然，因为是on-the-fly的，所以必须联网使用，而且MiKTeX只有Windows版本。 MiKTeX自带了一个TeXworks编辑器的，但是这软件用户体验并不好。我以前一直都用WinEdt，很好用，但是它是商业软件，我又不想盗版（说到底是没钱…），所以换了Texmaker。Texmaker可以媲美WinEdt，软件布局合理，各种快捷键用起来也很方便。不过在上手之前要简单配置一下。 如果是写英文文章，点击“快速构建”左边的箭头（或者F1快捷键），就能一键编译并刷新pdf视图。但是默认的快速构建使用的引擎是PdfLaTeX，如果你是中文用户，使用了xeCJK宏包，则必须使用XeLaTeX引擎编译，所以依次点击“选项->配置Texmaker->（左边）快速构建”，选择快速构建命令为”XeLaTeX + View PDF”。 构建好的PDF默认是以弹窗的形式展现的，我们可以设置让代码和PDF并排显示，这样方便在PDF和源代码之间切换，配置如下： Texmaker自带了一个PDF阅读器，当然你也可以使用外部阅读器，比如非常棒的Sumatra PDF，只需填入Sumatra PDF的路径跟上%.pdf，并选中External Viewer。 Texmaker还有一个很好用的功能是“正向/反向搜索”。反向搜索是点击PDF某个位置，会跳到tex源代码对应位置，快捷方式是ctrl+click。正向搜索是点击tex源代码某个位置，会跳到PDF对应的位置，默认快捷方式ctrl+space，但是这个快捷方式好像用不了，可以自行配置成其他快捷方式，比如ctrl+1，我当时是打开下图的快捷方式窗口才发现这个问题的。 正反向搜索都可以通过鼠标右键菜单实现，但是快捷键还是更方便的。最重要的一点是，源文件*.tex所在路径不能有中文！！！要不然正反向搜索不能用，这点很重要，我当时郁闷了好久。 另外还可以配置一下编辑器的字体，勾选”Backup documents every 10 min”之类的。 OK，大功告成，这种三段式的界面、F1快速构建以及正向/反向搜索，用起来真是太顺手了，Just Enjoy \(\LaTeX\)~ 下面是我常用的\(\LaTeX\)中文模板： LaTeXDemo.pdf LaTeXDemo.tex

SVM回顾 支持向量机（SVM）的一大特点是最大化间距（max margin）。对于如上图的二分类问题，虽然有很多线可以将左右两部分分开，但是只有中间的红线效果是最好的，因为它的可活动范围（margin）是最大的，从直观上来说很好理解。 对于线性二分类问题，假设分类面为 $$\begin{equation} u=\vec w \cdot \vec x-b \end{equation}$$则margin为 $$\begin{equation} m=\frac{1}{||w||_2} \end{equation}$$根据max margin规则和约束条件，得到如下优化问题，我们要求的就是参数\(\vec w\)和\(b\)： $$\begin{equation} \min\limits_{\vec w,b}\frac{1}{2}||\vec w||^2 \quad\text{subject to}\quad y_i(\vec w\cdot \vec x_i-b) \geq 1, \forall i,\end{equation}$$对于正样本，类标号\(y_i\)为+1，反之则为-1。根据拉格朗日对偶，(3)可以转换为如下的二次规划（QP）问题，其中\(\alpha_i\)为拉格朗日乘子。 $$\begin{equation} \min\limits_{\vec \alpha}\Psi(\vec\alpha)=\min\limits_{\vec \alpha}\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^Ny_iy_j(\vec x_i\cdot\vec x_j)\alpha_i\alpha_j-\sum_{i=1}^N\alpha_i,\end{equation}$$其中N为样本数量。上式还需满足如下两个约束条件： $$\begin{equation} \alpha_i\geq 0, \forall i,\end{equation}$$$$\begin{equation} \sum_{i=1}^Ny_i\alpha_i=0.\end{equation}$$一旦求解出所有的拉格朗日乘子，则我们可以通过如下的公式得到分类面参数\(\vec w\)和\(b\)。 $$\begin{equation}\vec w=\sum_{i=1}^Ny_i\alpha_i\vec x_i,\quad b=\vec w\cdot\vec x_k-y_k\quad\text{for some}\quad\alpha_k>0.\end{equation}$$当然并不是所有的数据都可以完美的线性划分，可能有少量数据就是混在对方阵营，这时可以通过引入松弛变量\(\xi_i\)得到软间隔形式的SVM： $$\begin{equation} \min\limits_{\vec w,b,\vec\xi}\frac{1}{2}||\vec w||^2+C\sum_{i=1}^N\xi_i \quad\text{subject to}\quad y_i(\vec w\cdot \vec x_i-b) \geq 1-\xi_i, \forall i,\end{equation}$$其中的\(\xi_i\)为松弛变量，能假装把错的样本分对，\(C\)对max margin和margin failures的trades off。对于这个新的优化问题，约束变成了一个box constraint： $$\begin{equation}0\leq\alpha_i \leq C,\forall i.\end{equation}$$而松弛变量\(\xi_i\)不再出现在对偶公式中了。 ...