统计学习方法

第一版

第二版

第 1 章 统计学习及监督学习概论

统计学习的主要特点是

  1. 统计学习以计算机及网络为平台,是建立在计算机及网络之上的;
  2. 统计学习以数据为研究对象,是数据驱动的学科;
  3. 统计学习的目的是对数据进行预测与分析;
  4. 统计学习以方法为中心,统计学习方法构建模型并应用模型进行预测与分析;
  5. 统计学习是概率论、统计学、信息论、计算理论、最优化理论及计算机科学等多个领域的交叉学科,并且在发展中逐步形成独自的理论体系与方法论。

假设空间(hypothesis space)
H={f(x;θ)θRD}orF={PP(YX;θ),θRD}\mathcal H = \{ f(x;\theta) | \theta \in \mathbb{R}^D\} \\ or \quad \mathcal F = \{P|P(Y|X;\theta),\theta \in \mathbb{R}^D\}
其中f(x;θ)f(x; \theta)是参数为θ\theta 的函数(决策函数),也称为模型(Model),参数向量θ\theta取值与DD维欧式空间RD\mathbb{R}^D,也称为参数空间(parameter space),DD 为参数的数量(维度)

模型的假设空间(hypothesis space)包含所有可能的条件概率分布或决策函数

特征空间(feature space)
每个具体的输入是一个实例(instance),通常由特征向量(feature vector)表示。这
时,所有特征向量存在的空间称为特征空间(feature space)。特征空间的每一维对应于
一个特征。

输入空间中的一个输入向量x=(x1,x2)x = (x_1,x_2),在多项式模型中特征向量是(x12,x1x2,x22,...x_1^2,x_1x_2,x_2^2,...)
一般说的线性模型,指的是特征向量的线性组合,而不是指输入向量,所以说模型都是定义在特征空间上的

统计学习的三要素

  1. 模型的假设空间(hypothesis space),简称:模型(model)。假设空间即我们对模型形式的先验假设,最终我们求得的模型必定符合我们对模型形式的先验假设。
  2. 模型选择的准则(evaluation criterion),简称:策略(strategy)或者学习准则。即我们用什么标准来评价一个模型的好坏。策略决定了我们从假设空间中选择模型的偏好。
  3. 模型学习的算法(algorithm),简称:算法(algorithm)。优化算法指的是通过什么样的方式调整我们的模型结构或模型超参数取值,使得模型的目标函数取值不断降低。优化算法决定了我们用什么样的步骤在假设空间中寻找合适的模型。

以线性回归(Linear Regression)为例:
模型: f(x;w,b)=wTx+bf(x;w,b) = w^Tx +b
策略(strategy)或者学习准则: 平方损失函数 L(y,y^)=(yf(x,θ))2\mathcal L(y,\hat{y}) = (y-f(x,\theta))^2
算法:解析解 analytical solution(闭式解 closed-form solution)和数值解 numerical solution,如:closed-form 的最小二乘的解以及梯度下降法

机器学习的定义

graph LR; F(["未知的目标函数(理想中完美的函数):𝑓: 𝒙⟶𝑦"])-->D["训练样本D:{(𝒙¹,𝑦¹),...,(𝒙ⁿ,𝑦ⁿ)}"]; D-->A{{"算法"}} H{{"假设空间"}}-->A A-->G["模型 g≈f"]

使用训练数据来计算接近目标 𝑓 的假设(hypothesis )g (来自:Machine Learning Foundations(机器学习基石)- the learning problem,25 页

监督学习
监督学习(supervised learning)是指从标注数据中学习预测模型的机器学习问题。本质是学习输入到输出的映射的统计规律

输入变量与输出变量均为连续变量的预测问题称为回归问题
输出变量为有限个离散变量的预测问题称为分类问题
输入变量与输出变量均为变量序列的预测问题称为标注问题(分类问题的推广,如:隐马尔可夫模型 HMM,条件随机场 CRF)。

监督学习的模型可以是概率模型或非概率模型,由条件概率分布P(YX)P(Y|X)决策函数(decision function)Y=f(X)Y=f(X)表示,随具体学习方法而定。对具体的输入进行相应的输出预测时,写作P(yx)P(y|x)Y=f(x)Y=f(x)
y=arg maxyP(yx)y =\displaystyle\argmax_{y} P(y|x)

联合概率分布
监督学习假设输入与输出的随机变量 X 和 Y 遵循联合概率分布P(X,Y)P(X,Y)P(X,Y)P(X,Y)表示分布函数,或分布密度函数。注意,在学习过程中,假定这一联合概率分布存在,但对学习系统来说,联合概率分布的具体定义是未知的。训练数据与测试数据被看作是依联合概率分布P(X,Y)P(X,Y)独立同分布产生的
统计学习假设数据存在一定的统计规律,XXYY具有联合概率分布的假设就是监督学习关于数据的基本假设。

非监督学习
非监督学习(unsupervised learning)是指从无标注数据中学习预测模型的机器学习问题。本质是学习数据中的统计规律或潜在结构

非监督学习的模型可以表示为函数z=g(x)z = g(x)或者条件概率分布P(zx)P(z|x) (输出zz可以是聚类或者降维
z=arg maxzP(zx)z =\displaystyle\argmax_{z} P(z|x)
以及 条件概率分布P(xz)P(x|z) (用来做概率密度估计,比如 GMM 中P(xz)P(x|z)属于高斯分布,如果假设知道数据来自哪个高斯分布,即知道zz,我们可以用极大似然估计来估计相关参数)。

核密度估计 Kernel Density Estimation. - 应用密度估计检测离群值(outlier)的LocalOutlierFactor

概率模型(probabilistic model)与非概率模型(non-probabilistic model)或者确定性模型(deterministic model)

概率模型(probabilistic model)- 条件概率分布 P(y|x)和 非概率模型(non-probabilistic model) - 函数 y=f(x)可以相互转化,条件概率分布最大化后得到函数,函数归一化后得到条件概率分布。所以概率模型与非概率模型的区别不在于输入输出之间的映射关系,而在于模型的内部结构:概率模型一定可以表示为联合概率分布的形式,而非概率模型则不一定存在这样的联合概率分布。

概率模型的代表是概率图模型(probabilistic graphical model)参考文献[13]^{参考文献[1-3]},联合概率分布可以根据图的结构分解为因子乘积的形式,可以用最基本的加法规则和乘法规则进行概率推理:
P(x)=yP(x,y)P(x,y)=P(x)P(yx)P(x) = \sum_yP(x,y) \\ P(x,y) = P(x)P(y|x)

参数化模型(parametric model)和非参数化模型(non-parametric model)

参数化模型假设模型参数的维度固定,模型可以由有限维参数完全刻画,不随数据点的变化而变化。(如:感知机、GMM、logistic regression、朴素贝叶斯、k 均值聚类、潜在语义分析、概率潜在语义分析、潜在狄利克雷分配)
非参数化模型假设模型参数的唯独不固定或者说无穷大,随着训练数据量的增加而不断增大。(如:决策树、支持向量机、AdaBoost、k 近邻)

非参数化模型意味着决策树没有假设空间分布和分类器结构?

在线学习(online learning)和批量学习(batch learning)

在线学习每次接受一个样本,预测后学习模型,并不断重复该操作。
批量学习一次接受所有数据,学习模型之后进行预测。

在线学习比批量学习更难,因为每次模型更新中可利用的数据有限。

贝叶斯学习(Bayesian learning)/ 贝叶斯推理(Bayesian inference)
Bayes  Rule:P(XY)posterior=P(YX)likelihoodP(X)priorP(Y)evidence=P(YX)likelihoodP(X)priorxP(YX)P(X)evidence\mathrm{Bayes \; Rule:} \\ \underbrace{P(X|Y)}_{\mathrm{posterior}} = \frac{\overbrace{P(Y|X)}^{\mathrm{likelihood}}\overbrace{P(X)}^{\mathrm{prior}}}{\underbrace{P(Y)}_{\mathrm{evidence}}} = \frac{\overbrace{P(Y|X)}^{\mathrm{likelihood}}\overbrace{P(X)}^{\mathrm{prior}}}{\underbrace{\sum_{x}P(Y|X)P(X)}_{\mathrm{evidence}}}

核技巧(kernel trick)/ 核方法(kernel method)

核方法是一类把低维空间的非线性可分问题,转化为高维空间的线性可分问题的方法。
核技巧是一种利用核函数直接计算 ϕ(x),ϕ(z)\lang \phi(x),\phi(z) \rang ,以避开分别计算 ϕ(x)\phi(x)ϕ(z)\phi(z) ,从而加速核方法计算的技巧。

核函数Kernel function
X\mathcal X 是输入空间(即 xiXx_i \in \mathcal XX\mathcal XRn\mathbb R^n 的子集或离散集合 ),又设 H\mathcal H 为特征空间(​ 希尔伯特空间附加知识:各种空间介绍^{附加知识:各种空间介绍}),如果存在一个从 X\mathcal XH\mathcal H 的映射

ϕ(x):XH\phi(x) : \mathcal X \to \mathcal H

使得对所有 x,zXx,z \in \mathcal X ,函数 K(x,z)K(x,z) 满足条件

K(x,z)=ϕ(x).ϕ(z)=ϕ(x),ϕ(z)K(x,z) = \phi(x).\phi(z) = \lang \phi(x),\phi(z) \rang

则称 K(x,z)K(x,z) 为核函数。其中 ϕ(x)\phi(x) 为映射函数, ϕ(x),ϕ(z)\lang \phi(x),\phi(z) \rang 为内积。

核技巧的想法是,在学习和预测中只定义核函数 K(x,z)K(x,z) ,而不显式地定义映射函数 ϕ\phi。通常直接计算K(x,z)K(x,z)比较容易,而通过ϕ(x)\phi(x)ϕ(z)\phi(z)计算K(x,z)K(x,z)并不容易。

注意:ϕ\phi是输入空间Rn\mathbb{R}^n到特征空间H\mathcal H的映射,特征空间H\mathcal H一般是高维的,甚至是无穷维的。所以ϕ\phi不好计算,甚至会带来维度灾难又称维度诅咒(Curse of Dimensionality)附加知识:维度诅咒^{附加知识:维度诅咒}

附加知识

正则化

正则化符合奥卡姆剃刀(Occam's razor)原理。

参考:L1L2 正则化和凸优化

模型选择

参考:模型选择

生成模型和判别模型

参考:生成模型和判别模型

各种空间介绍

线性空间就是定义了加法和数乘的空间(空间里的一个元素就可以由其他元素线性表示)。


度量空间就是定义了距离的空间(曼哈顿距离,欧氏距离,闵可夫斯基距离,马氏距离,切比雪夫距离)。
定义距离时,有三条公理必须遵守:

  1. 非负性、同一性:dist(xi,xj)0dist(x_i,x_j) \geq 0(非负性),dist(xi,xj)=0dist(x_i,x_j) = 0当且仅当xi=xjx_i=x_j(同一性)
  2. 对称性:dist(xi,xj)=dist(xj,xi)dist(x_i,x_j) = dist(x_j,x_i)
  3. 三角不等式(也叫直递性):dist(xi,xj)dist(xi,xk)+dist(xk,xj)dist(x_i,x_j) \leq dist(x_i,x_k) + dist(x_k,x_j)
    希尔伯特空间(Hilbert)

    文字解释:【两点之间距离不为负;两个点只有在 空间 上重合才可能距离为零;a 到 b 的距离等于 b 到 a 的距离;a 到 c 的距离加上 c 到 b 的距离大于等于 a 直接到 b 的距离;】


赋范空间就是定义了范数的空间。
x 的范数||x||就是 x 的长度。那么这里的长度和上一节中说的距离到底有什么区别呢。距离的概念是针对两个元素来说的,例如 d(x,y)指的是 x 与 y 两个元素之间的距离,而范数是针对一个元素来说的,每一个元素都对应一个范数,可以将范数理解为一个元素到零点的距离(这只是一种理解,并不是定义),也就是它自己的长度。
定义:
称 映射.:RnR||.|| : \mathbb{R}^n \to \mathbb{R}Rn\mathbb{R}^n 上的范数,当且仅当:

  1. 非负性: xRn,x0\forall x \in \mathbb{R}^n ,||x|| \geq 0 ,x=0||x|| = 0当且仅当x=0x=0
  2. 数乘:xRn,aRn,ax=a.x\forall x \in \mathbb{R}^n ,a \in \mathbb{R}^n, ||ax|| = |a|.||x||
  3. 三角不等式: x,yRn,x+yx+y\forall x,y \in \mathbb{R}^n ,||x+y|| \leq ||x|| + ||y||

如果我们定义了范数,可以在这基础上定义距离:dist(x,y)=||x-y||。根据范数的三条性质,我们可以证明我们这样定义的距离也满足距离的定义,聪明的你可以自己证明一下(对称性的证明,提一个-1 出来,一加绝对值就是 1 了)。

也就是说范数其实是一个更加具体的概念,有了范数一定能利用范数定义距离,但是有距离不能定义范数

也许你会问,你不是说理解范数就是一个元素到零点的距离吗,那定义范数为||x||=dist(x,0) 不就行了吗。这样的话,对于范数的第二条性质就不一定会满足,||ax||=dist(ax,0),而 dist(ax,0)不一定等于|a|dist(x,0),具体等不等于还要看你的距离是怎么定义的。

例如:Lp范数
欧式距离对应 L2 范数
曼哈顿距离对应 L1 范数
切比雪夫距离对应 L∞ 范数
Lp范数:当 p>=1 时,向量的 Lp范数是凸的。(这也是为什么一般不用 L0 范数的原因之一)


线性赋范空间就是定义了加法、数乘和范数的空间。


巴拿赫空间就是完备的赋范线性空间。(Banach space)
完备的空间的定义:如果一个空间是完备的,那么该空间中的任何一个柯西序列都收敛在该空间之内。

首先来说一下柯西序列是什么,柯西序列就是随着序数增加,值之间的距离越来越小的序列。换一种说法是,柯西序列可以在去掉有限个值之后,使任意两个值之间的距离\underline{\mathrm{距离}}都小于任意给定正常数(其实这就是定义了一个极限而已)。

那么任意一个柯西序列都收敛在该空间内是什么意思呢,举个例子你就明白了。

设定义在有理数空间 Q 上的序列:xn=[2n]nx_n = \frac{[\sqrt{2}n]}{n},其中[x]表示 x 取整数部分。
对于这个数列来说,每一个元素的分子分母都是整数,所以每一个xnx_n都在有理数空间 Q 上,那这个序列的极限呢,稍有常识的人都能看出,这个序列的极限是2\sqrt{2},而这并不是一个有理数,所以这个柯西序列的极限不在该空间里面,也就是说有理数空间 Q 是不完备的。

所以完备的意义我们可以这样理解,那就是在一个空间上我们定义了极限,但是不论你怎么取极限,它的极限的值都不会跑出这个空间,那么这个空间就是完备空间

另外,不知道你有没有发现,上面在解释什么是柯西序列的时候,有一个词我加了下划线,那就是距离,也就说说在定义完备空间之前,要先有距离的概念。所以完备空间,其实也是完备度量空间

所以,巴拿赫空间满足几条特性呢:距离、范数、完备。


内积空间就是定义了内积的空间。Inner product space
有时也称准希尔伯特空间。
内积就是我们所说的点乘、标积,它的定义方式也不是唯一的,但如同距离范数的定义一样,内积的定义也要满足某些条件,不能随便定义。

定义映射.,.:V×VF\lang .,. \rang : V \times V \to \mathbb{F}, 其中VV是向量,F\mathbb{F}是标量
x,y,zV,sFx,y,z \in V ,s \in \mathbb{F},那么内积满足

  1. 第一个参数中的线性:
    sx,y=sx,yx+y,z=x,z+y,z0,x=0\lang sx,y \rang = s\lang x,y \rang \\ \lang x+y,z \rang = \lang x,z \rang + \lang y,z \rang \\ \lang 0,x \rang = 0

  2. 共轭对称:x,y=y,x\lang x,y \rang = \overline{\lang y,x \rang }

  3. 正定性:x,x>0if  x0\lang x,x \rang > 0 \quad\mathrm{if}\; x \neq 0

  4. 正半定性或非负定性:x,x,x0\forall{x}, \lang x,x \rang \geq 0

  5. 确定性:x,x=0必然有x=0\lang x,x \rang = 0 必然有 x=0

3,4,5 可以跟上面定义范数和距离一样写成一个

例子-欧几里得向量空间:
x,yRn,x,y=xTy=_i=1nxiyix,y \in \mathbb{R}^n , \lang x,y \rang = x^Ty=\sum\_{i=1}^n{x_iy_i}

只有定义了内积,才会有夹角的概念,才会有正交的概念,另外内积也可以定义范数,也就是说内积是比范数更具体的一个概念。


欧式空间就是定义了内积的有限维实线性空间。


希尔伯特空间就是完备的内积空间。(Hilbert space)
希尔伯特空间中的元素一般是函数,因为一个函数可以视为一个无穷维的向量。

graph LR; LS(("Linear Space"))-->NLS(("Normed Linear Space")); NLS-->BS(("Banach Space")) NLS-->IPS(("Inner Product Space")) IPS-->HS(("Hilbert Space")) IPS-->ES(("Euclid Space"))

参考:一片文章带你理解再生核希尔伯特空间(RKHS)以及各种空间

维度诅咒

维度诅咒通常是指在涉及到向量的计算的问题中,随着维数的增加,计算量呈指数倍增长的一种现象。高维度有更大的特征空间,需要更多的数据才可以进行较准确的估计。

若特征是二值的,则每增加一个特征,所需数据量都在以 2 的指数级进行增长,更何况很多特征不只是二值的。

几何角度 1:

background Layer 1 0.5

上图表示一个多维空间(以二维为例),设正方形边长为 1,则其内切圆半径为r=0.5r=0.5,则正方形面积为 1,内切圆面积为π(0.5)2\pi(0.5)^2 。若将此变为三维情况下,正方体体积为 1,内切球体积为43π(0.5)3\frac{4}{3}\pi(0.5)^3

因此球体的体积可以表示为V(d)=πd/2Γ(d2+1)0.5d=k(0.5)dV(d) = \frac{\pi^{d/2}}{\varGamma(\frac{d}{2}+1)}0.5^d = k(0.5)^d(d 为维度),则 limdk(0.5)d=0\lim_{d \to \infty}k(0.5)^d = 0,其内切超球体的体积为 0。由此可知,高维情况下,数据大都分布在四角(正方形内,内切圆外),稀疏性太大,不好分类。

维度越大,超球体体积越小。说明落在超球体内的样本越少,因为超球体是超立方体的内切球。不在球内,那只能在角落!

几何角度 2:

background Layer 1

上图也表示一个多维空间(以二维为例),则其中图形的体积有如下关系:外圆半径r=1r=1,内圆半径为rεr−\varepsilon 。同样在高维情况下,外圆体积为V外圆=k.1d=kV_{外圆} = k.1^d = k,中间的圆环体积为V圆环=kk(1ε)dV_{圆环} = k - k(1-\varepsilon)^d,则:
limdV圆环V外圆=limdkk(1ε)dk=limd(1(1ε)d)=1\lim_{d \to \infty}\frac{V_{圆环}}{V_{外圆}} = \lim_{d \to \infty}\frac{ k - k(1-\varepsilon)^d}{k} = \lim_{d \to \infty}(1-(1-\varepsilon)^d) = 1

高维情况下,无论ε\varepsilon多小,只要 d 足够大,圆环几乎占据了整个外圆,内圆体积趋向于 0,导致数据稀疏

参考:
The Curse of Dimensionality in classification
机器学习-白板推导系列(五)-降维(Dimensionality Reduction)

不等式(Inequality)

所有不等式 以及所有概率(Probabilistic)不等式

由 KL divergence 就能证明
DKL(PQ)i=1npilogpiqi0.{\displaystyle D_{\mathrm {KL} }(P\|Q)\equiv \sum _{i=1}^{n}p_{i}\log {\frac {p_{i}}{q_{i}}}\geq 0.}

概率不等式 Probabilistic inequalities

参考:初等数学学习笔记

参考文献

[1-1] Hastie T,Tibshirani R,Friedman J. The Elements of Statistical Learning: DataMining,Inference,and Prediction. Springer. 2001(中译本:统计学习基础——数据挖掘、推理与预测。范明,柴玉梅,昝红英等译。北京:电子工业出版社,2004)

[1-2] Bishop M. Pattern Recognition and Machine Learning. Springer,2006

[1-3] Probabilistic Graphical Models: Principles and Techniques by Daphne Koller, Nir Friedman from The MIT Press

[1-4] Deep Learning (Ian Goodfellow, Yoshua Bengio, Aaron Courville)

[1-5] Tom M Michelle. Machine Learning. McGraw-Hill Companies,Inc. 1997(中译本:机器学习。北京:机械工业出版社,2003)

[1-6] Bayesian Reasoning and Machine Learning by David Barber 2007–2020 ,other version

[1-7] Reinforcement Learning:An Introduction (second edition 2020) by Richard S. Sutton and Andrew G. Barto ,other version

[1-8] 周志华,机器学习,清华大学出版社 (手推笔记 以及 公式推导解析)

[1-9] Lecture Notes in MACHINE LEARNING Dr V N Krishnachandran

第 2 章 感知机

判别模型

感知机Perceptron神经网络支持向量机的基础。最早在 1957 年由 Rosenblatt 提出参考文献[21]^{参考文献[2-1]}。Novikoff参考文献[22]^{参考文献[2-2]},Minsky 与 Papert参考文献[23]^{参考文献[2-3]}等人对感知机进行了一系列理论研究。感知机的扩展学习方法包括口袋算法(pocket algorithm)参考文献[24]^{参考文献[2-4]}、表决感知机(voted perceptron)参考文献[25]^{参考文献[2-5]}、带边缘感知机(perceptron with margin)参考文献[26]^{参考文献[2-6]}等。
Brief History of Machine Learning

要求:数据集线性可分(linearly separable data set)

感知机是一种线性分类模型,属于判别模型。感知机模型的假设空间是定义在特征空 间中的所有线性分类模型(linear classification model)或线性分类器(linear classifier),即 函数集合{ff(x)wx+b}\{f|f(x)=w·x+b\}

超平面 S:w.x+b=0w.x+b = 0,其中ww是 S 的法向量,bb是 S 的截距,超平面 S 称为分离超平面(separating hyperplane)

函数间隔:y(w.x+b)y(w.x + b)
几何间隔:1ww.x+b\frac{1}{||w||}|w.x + b| (在上面的 loss function 中没有考虑1w\frac{1}{||w||})

  1. 初始化参数(随机法):w0,b0w_0,b_0
  2. 选取数据(xi,yi)(x_i,y_i)
  3. 如果(xi,yi)(x_i,y_i)是误分类点,也就是yi(w.xi+b)0y_i(w.x_i + b) \leq 0,则对w,bw,b进行更新
    (xi,yi)点处梯度为:wL(w,b)=yixibL(w,b)=yi更新wwk+1wk+ηyixi更新bbk+1bk+ηyi其中学习率η(0,1]在(x_i,y_i)点处梯度为:\\ \nabla_wL(w,b) = -y_ix_i \\ \nabla_bL(w,b) = -y_i\\ 更新w:w_{k+1} \gets w_{k}+\eta y_ix_i \\ 更新b:b_{k+1} \gets b_{k}+\eta y_i \\其中学习率\eta \in (0,1]
  4. 循环 2-3,直到训练集中没有误分类点。

Novikoff 定理:
设训练集T={(x1,y1),...,(xN,yN)}T = \{(x_1,y_1),...,(x_N,y_N)\}是线性可分的,

  1. 设完美超平面w^opt.x^=0,w^opt=1\hat{w}_{opt}.\hat{x} = 0 , ||\hat{w}_{opt}||=1 将训练集完全正确分开(简化起见 w^opt.x^=wopt.x+b\hat{w}_{opt}.\hat{x} = w_{opt}.x +b),存在γ>0\gamma >0 ,对所有点有yi(w^opt.xi^)γy_i(\hat{w}_{opt}.\hat{x_i}) \geq \gamma

  2. R=max1iNxi^R = \max_{1\leq i\leq N}||\hat{x_i}||,则算法会在有限步 k 满足不等式k(Rγ)2k \leq (\frac{R}{\gamma})^2

证明(注意:带 hat 的表示扩充向量):

  1. 因为数据线性可分,对于所有点yi(w^opt.xi^)>0y_i(\hat{w}_{opt}.\hat{x_i}) > 0,所以存在
    γ=miniyi(w^opt.xi^)yi(w^opt.xi^)(2-1)\gamma = \min_i{y_i(\hat{w}_{opt}.\hat{x_i})} \leq {y_i(\hat{w}_{opt}.\hat{x_i})} \label{2-1}\tag{2-1}
    所以这里的γ\gamma代表了所有点离完美超平面的最小距离;

  2. 为了方便计算 设 扩充向量w^=(wT,b)T\hat{w} = (w^T,b)^T, 有
    w^k=w^k1+ηyixi^(2-2)\hat{w}_{k} = \hat{w}_{k-1}+\eta y_i\hat{x_i} \label{2-2}\tag{2-2}

  3. 推导不等式
    w^k.w^optkηγ(2-3)\hat{w}_{k}.\hat{w}_{opt} \geq k\eta\gamma \label{2-3}\tag{2-3}

(2-1)\eqref{2-1}(2-2)\eqref{2-2}
w^k.w^opt=w^k1.w^opt+ηyiw^opt.xi^w^k1.w^opt+ηγw^k2.w^opt+2ηγkηγ\hat{w}_{k}.\hat{w}_{opt} = \hat{w}_{k-1}.\hat{w}_{opt} + \eta{y_i}\hat{w}_{opt}.\hat{x_i} \\ \geq \hat{w}_{k-1}.\hat{w}_{opt} + \eta\gamma \\ \geq \hat{w}_{k-2}.\hat{w}_{opt} + 2\eta\gamma \\ \geq k\eta\gamma

  1. 推导不等式
    w^k2kη2R2(2-4)||\hat{w}_{k}||^2 \leq k\eta^2R^2 \label{2-4}\tag{2-4}
    (2-2)\eqref{2-2}
    w^k2=w^k1+ηyixi^2=w^k12+2ηyiw^k1.x^i+η2x^i2||\hat{w}_{k}||^2=||\hat{w}_{k-1}+\eta y_i\hat{x_i}||^2 = ||\hat{w}_{k-1}||^2 + 2\eta{y_i}\hat{w}_{k-1}.\hat{x}_{i} + \eta^2||\hat{x}_{i}||^2
    假设 k 次完全分对,那么 k-1 次有误分类点,则yiw^k1.x^i0{y_i}\hat{w}_{k-1}.\hat{x}_{i} \leq 0
    所以
    w^k2=w^k12+2ηyiw^k1.x^i+η2x^i2w^k12+η2x^i2w^k12+η2R2w^k22+2η2R2...kη2R2||\hat{w}_{k}||^2 =||\hat{w}_{k-1}||^2 + 2\eta{y_i}\hat{w}_{k-1}.\hat{x}_{i} + \eta^2||\hat{x}_{i}||^2 \\ \leq ||\hat{w}_{k-1}||^2 + \eta^2||\hat{x}_{i}||^2 \\ \leq ||\hat{w}_{k-1}||^2 + \eta^2R^2 \\ \leq ||\hat{w}_{k-2}||^2 + 2\eta^2R^2 \leq ... \\ \leq k\eta^2R^2

  2. (2-3)\eqref{2-3}(2-4)\eqref{2-4}

kηγw^k.w^optw^k.w^opt=1柯西-施瓦茨 (Cauchy–Schwarz) 不等式kηR  k2γ2kR2k(Rγ)2k\eta\gamma \leq \underbrace{\hat{w}_{k}.\hat{w}_{opt} \leq ||\hat{w}_{k}||.\underbrace{||\hat{w}_{opt}||}_{=1} }_{\text{柯西-施瓦茨 (Cauchy–Schwarz) 不等式}} \leq \sqrt{k} \eta R \\ \; \\ \Rightarrow k^2\gamma^2 \leq kR^2 \\ \Rightarrow k \leq (\frac{R}{\gamma})^2

也就是说 k 是有上界的。

书中还介绍了原形式的对偶形式,也就是等价形式(SVM 中 7.2.2 节 127 页也是等价的意思,区别于拉格朗日对偶),这两个地方的等价都是经过基本推导,求出 w 参数,然后对原问题进行了替换。

参考文献

[2-1] Rosenblatt, F. (1958). The perceptron: A probabilistic model for information storage and organization in the brain. Psychological Review, 65(6), 386–408.

[2-2] Novikoff, A. B. (1962). On convergence proofs on perceptrons. Symposium on the Mathematical Theory of Automata, 12, 615-622. Polytechnic Institute of Brooklyn.

[2-3] Minsky M L and Papert S A 1969 Perceptrons (Cambridge, MA: MIT Press)

[2-4] Gallant, S. I. (1990). Perceptron-based learning algorithms. IEEE Transactions on Neural Networks, vol. 1, no. 2, pp. 179-191.

[2-5] Freund, Y. and Schapire, R. E. 1998. Large margin classification using the perceptron algorithm. In Proceedings of the 11th Annual Conference on Computational Learning Theory (COLT' 98). ACM Press.

[2-6] Li YY,Zaragoza H,Herbrich R,Shawe-Taylor J,Kandola J. The Perceptron algorithmwith uneven margins. In: Proceedings of the 19th International Conference on MachineLearning. 2002,379–386

[2-7] Widrow, B., Lehr, M.A., "30 years of Adaptive Neural Networks: Perceptron, Madaline, and Backpropagation," Proc. IEEE, vol 78, no 9, pp. 1415-1442, (1990)。

[2-8] Cristianini N,Shawe-Taylor J. An Introduction to Support Vector Machines and OtherKernelbased Learning Methods. Cambridge University Press,2000

第 3 章 k 近邻法

判别模型

k 近邻法(k-nearest neighbor,k-NN)1968 年由 Cover 和 Hart 提出,是一种基本分类与回归方法。本书只讨论分类问题中的 k 近邻法。
k 值的选择、距离度量及分类决策规则是 k 近邻法的三个基本要素。
最后讲述 k 近邻法的一个实现方法——kd 树,介绍构造 kd 树和搜索 kd 树的算法

k 近邻法的三个基本要素
k 值的选择:超参数,可以使用交叉验证法来选取最优 k 值
距离度量:L2L_2距离/欧氏距离,LpL_p距离/Minkowski 距离
分类决策规则:多数表决(0-1 损失也就是指示函数)

附加知识

距离度量

Distance

sklearn.neighbors.DistanceMetric

Distance computations(scipy.spatial.distance)

24 种距离度量小结

先了解度量空间和赋范空间

实值向量空间的度量:

实值向量空间的度量(scipy):

整数值向量空间的度量:

布尔值向量空间的度量:

经纬度距离:

其它:

参考文献

[3-1] Cover T,Hart P. Nearest neighbor pattern classification. IEEE Transactions onInformation Theory,1967

[3-2] Hastie T,Tibshirani R,Friedman J. The Elements of Statistical Learning: DataMining,Inference,and Prediction,2001(中译本:统计学习基础——数据挖掘、推理与预测。范明,柴玉梅,昝红英等译。北京:电子工业出版社,2004)

[3-3] Friedman J. Flexible metric nearest neighbor classification. Technical Report,1994

[3-4] Weinberger KQ,Blitzer J,Saul LK. Distance metric learning for large margin nearestneighbor classification. In: Proceedings of the NIPS. 2005

[3-5] Samet H. The Design and Analysis of Spatial Data Structures. Reading,MA: Addison-Wesley,1990

第 4 章 朴素贝叶斯法

朴素贝叶斯(Naïve Bayes)法是基于贝叶斯定理特征条件独立假设(Naive 天真的)的分类方法。
对于给定的训练数据集,首先基于特征条件独立假设学习输入/输出的联合概率分布;然后基于此模型,对给定的输入 x,利用贝叶斯定理求出后验概率最大的输出 y。
朴素贝叶斯法实现简单,学习与预测的效率都很高,是一种常用的方法。并且支持 online learning(有 partial_fit 方法)。

朴素贝叶斯法是典型的生成学习方法。生成方法由训练数据学习联合概率分布 P(X,Y),然后求得后验概率分布 P(Y|X)。具体来说,利用训练数据学习 P(X|Y)和 P(Y)的估计,得到联合概率分布:P(X,Y)= P(Y)P(X|Y) ;概率估计方法可以是极大似然估计或贝叶斯估计。

贝叶斯定理(Bayes' theorem)
P(AB)=P(BA)P(A)P(B)P(A|B) = \frac{P(B|A)P(A)}{P(B)}

特征条件独立假设
条件独立
(ABC)    P(AB,C)=P(AC)(ABC)    P(A,BC)=P(AC)P(BC)(A\perp B|C) \iff P(A|B,C) = P(A|C) \\ (A\perp B|C) \iff P(A,B|C) = P(A|C)P(B|C)

特征条件独立假设就是已知 y 的情况下,x 中每个特征相互独立。

数据集T={(x1,y1),...,(xN,yN)}T = \{(x_1,y_1),...,(x_N,y_N)\}KK为类别个数,其中xix_i是 n 维向量xi=(xi(1),...,xi(n))Tx_i = (x_i^{(1)},...,x_i^{(n)})^T

附加知识

参数估计

参数估计(Parameter Estimation) 有点估计(point estimation)和区间估计(interval estimation)两种

点估计法:

共轭先验(Conjugate prior:如果先验分布 prior 和后验分布 posterior 属于同一分布簇,则 prior 称为 likehood 的共轭先验
likehood 为高斯分布,prior 为高斯分布,则 posterior 也为高斯分布。
likehood 为伯努利分布(二项式分布),prior 为 beta 分布,则 posterior 也为 beta 分布。
likehood 为多项式分布,prior 为 Dirichlet 分布(beta 分布的一个扩展),则 posterior 也为 Dirichlet(狄利克雷)分布。beta 分布可以看作是 dirichlet 分布的特殊情况。

最小二乘估计(Least squares estimation, LSE)

矩估计(Method of moments estimators)

区间估计法:
区间估计最流行的形式是置信区间 confidence intervals (一种频率论方法)和 可信区间 credible intervals(一种贝叶斯方法),此外还有预测区间(Prediction interval)等

采样法: 贝叶斯估计,近似推断
马尔可夫链蒙特卡罗法 Markov chain Monte Carlo, MCMC

参考文献

[4-1] Mitchell TM. Chapter 1: Generative and discriminative classifiers: Naïve Bayes andlogistic regression. In: Machine Learning. Draft,2005.

[4-2] Hastie T,Tibshirani R,Friedman J. The Elements of Statistical Learning. DataMining,Inference,and Prediction. Springer-Verlag,2001(中译本:统计学习基础——数据挖掘、推理与预测。范明,柴玉梅,昝红英等译。北京:电子工业出版社,2004)

[4-3] Bishop C. Pattern Recognition and Machine Learning,Springer,2006

第 5 章 决策树

判别模型

决策树(decision tree)是一种基本的分类与回归方法,具有良好的可解释性(可视化),通常包括 3 个步骤:特征选择、决策树的生成和决策树的修剪

特征选择
特征选择在于选取对训练数据具有分类能力的特征。(sklearn 中可以返回 feature_importances_特征重要性,属性越重要,特征空间划分的面积越大)

也就是计算每个特征的(信息增益,基尼指数)来选择特征(作为根节点)进行特征空间划分,注意:划分后再次计算每个特征的(信息增益,基尼指数),除非该特征所在的空间就只有一类了(也就是该特征不可分了,那么就直接生成叶子节点);

特征选择的准则:

决策树的生成
常见算法(参见:Decision tree learning以及Tree algorithms):

决策树的修剪 Decision tree pruning
修剪是机器学习和搜索算法中的一种数据压缩技术,它通过删除树中对分类实例不重要和冗余的部分来减小决策树的大小。剪枝降低了最终分类器的复杂度,从而通过减少过拟合来提高预测精度。

统计学习方法三要素:

Overview of Decision Trees

Decision Tree

Decision Trees

附加知识

信息论(Information Theory

Entropy, Relative Entropy, Cross Entropy

熵(Entropy

在信息论中,熵用来衡量一个随机事件的不确定性。也叫香农熵 Shannon's(人名) entropy。
H(X)=Ep(x)[I(X)]=Ep(x)[logp(x)]=i=1np(xi)logp(xi)=Xp(x)logp(x)dxH(X) = E_{p(x)}[I(X)] = E_{p(x)}[-\log {p(x)}] \\= -\sum_{i=1}^n {p(x_i)} \log {p(x_i)} \\= -\int_{X} {p(x)} \log {p(x)} dx
其中I(X)=logp(x)I(X) = -\log {p(x)} 称为自信息Self Information),是一个随机事件所包含的信息量。一个随机事件发生的概率越高,其自信息越低。如果一个事件必然发生,其自信息为 0。
在自信息的定义中,对数的底可以使用 2、自然常数 𝑒 或是 10。当底为 2 时,自信息的单位为 bit;当底为 𝑒 时,自信息的单位为 nat。

熵越高,则随机变量的信息越多(不确定性越大,系统越复杂);熵越低,则随机变量的信息越少。


求最大熵:假设概率分布

X 1 2 ... n
p(x) p₁ p₂ ... pⁿ

maxH(p)=maxi=1npilogpis.t.i=1npi=1\max H(p) = \max -\sum_{i=1}^n p_i \log p_i \\ s.t. \sum_{i=1}^n p_i = 1

由拉格朗日乘数法(Lagrange Multiplier),最大变最小时去掉负号
L(p,λ)=i=1npilogpi+λ(1i=1npi)偏导:Lpi=logpi+pi.1piλ令偏导为0得:pi=exp(λ1)\mathcal L(p,\lambda) = \sum_{i=1}^n p_i \log p_i + \lambda(1-\sum_{i=1}^n p_i) \\偏导:\frac{\partial\mathcal L}{\partial p_i} = \log p_i + p_i.\frac{1}{p_i} - \lambda \\ 令偏导为0得:p_i^*=exp(\lambda-1)

因为λ\lambda是一个超参数(常数),所以pip_i^*是一个常数,所以 p1=p2=...=pn=1np_1^*=p_2^*=...=p_n^*=\frac{1}{n}

所以概率分布为一个均匀分布,则熵最大,由此性质我们来证明熵的取值范围:设 p 是一个均匀分布p=1np = \frac{1}{n}
H(p)=i=1n1nlog1n=i=1n1nlogn1=i=1n1nlogn=lognH(p) = -\sum_{i=1}^n \frac{1}{n} \log \frac{1}{n} \\= -\sum_{i=1}^n \frac{1}{n} \log n^{-1} \\= \sum_{i=1}^n \frac{1}{n} \log n \\= \log n
所以:0H(p)logn0 \leq H(p) \leq \log n


已知连续随机变量的均值为μ\mu,方差为σ2\sigma^2,求熵最大对应的概率分布:
arg maxp(x)p(x)logp(x)dxs.t.p(x)dx=1xp(x)dx=μ(xμ)2p(x)dx=σ2\argmax_{p(x)} -\int p(x)\log p(x)dx \\ s.t. \int p(x)dx =1 \\ \int xp(x)dx = \mu \\ \int (x-\mu)^2p(x)dx=\sigma^2
拉格朗日函数
L(p(x),λ1,λ2,λ3)=p(x)logp(x)dx+λ1(p(x)dx1)+λ2(xp(x)dxμ)+λ3((xμ)2p(x)dxσ2)L(p(x),\lambda_1,\lambda_2,\lambda_3) = -\int p(x)\log p(x)dx +\lambda_1(\int p(x)dx - 1)+\lambda_2(\int xp(x)dx - \mu) +\lambda_3(\int (x-\mu)^2p(x)dx - \sigma^2)
F(p)=(logp(x)+λ1+λ2x+λ3(xμ)2)p(x)F(p)=(-\log p(x) + \lambda_{1} +\lambda_{2}x+ \lambda_{3}(x-\mu)^{2})p(x)
求偏导并令其为0(可以把求积分当做求和,这样求偏导就容易想象了)
Lp(x)=[logp(x)+1]+λ1+λ2x+λ3(xμ)2\frac{\partial L}{\partial p(x)} = -[\log p(x)+1]+\lambda_1+\lambda_2x+\lambda_3(x-\mu)^2

p(x)=exp{λ11+λ2x+λ3(xμ)2}p(x) = \exp\{\lambda_1-1+\lambda_2x+\lambda_3(x-\mu)^2\}
把跟x有关的保留,其它的设为常数,有
p(x)=exp{λ11+λ2x+λ3(xμ)2}=e1+λ1eλ2x+λ3(xμ)2=Ceλ2x+λ3(xμ)2=Ceλ3(x22(μλ22λ3)x+u2)=Ceλ3(xμ+λ22λ3)2=C.exp{λ3(xμ+λ22λ3)2}p(x) = \exp\{\lambda_1-1+\lambda_2x+\lambda_3(x-\mu)^2\}\\ =e^{-1+\lambda_{1}}\cdot e^{\lambda_{2}x+ \lambda_{3}(x-\mu)^{2}}=C e^{\lambda_{2}x+ \lambda_{3}(x-\mu)^{2}} \\ = Ce^{\lambda_{3}(x^{2} -2(\mu-\frac{\lambda_{2}}{2\lambda_{3}})x+ u^{2})} = C e^{\lambda_{3}(x -\mu+ \frac{\lambda_{2}}{2\lambda_{3}})^{2}} \\= C.\exp\{\lambda_3(x-\mu+\frac{\lambda_2}{2\lambda_3})^2\}

根据(xμ+λ22λ3)2(x-\mu+\frac{\lambda_2}{2\lambda_3})^2得到p(x)p(x)关于μλ22λ3\mu - \frac{\lambda_{2}}{2\lambda_{3}}对称(偶函数关于x=0对称p(x)=p(x)p(x) = p(-x)),所以E[p(x)]=μλ22λ3=μE[p(x)] = \mu - \frac{\lambda_{2}}{2\lambda_{3}} = \mu,得λ2=0\lambda_{2} = 0

那么
p(x)=Ceλ3(xμ)2p(x)= C e^{\lambda_{3}(x -\mu)^{2}}
因为 p(x)>0p(x)>0 ,所以 C>0,λ3<0C>0,\lambda_{3}<0
λ=λ3\lambda = -\lambda_3
根据积分为1的约束 以及ex22dx=2π\int e^{-\frac{x^2}{2}}dx = \sqrt{2\pi},得:
p(x)dx=1=Ceλ(xμ)2dx=Cπλ\int p(x)dx = 1 = C\int e^{-\lambda(x -\mu)^{2}} dx = C\sqrt{\frac{\pi}{\lambda}}
得到C=λπC = \sqrt{\frac{\lambda}{\pi}}
根据方差的约束,得:
(xμ)2p(x)dx=σ2=(xμ)2eλ(xμ)2dx=Cπλ.12λ=12λ\int (x-\mu)^2p(x)dx=\sigma^2 = \int (x-\mu)^2e^{-\lambda(x -\mu)^{2}}dx = C\sqrt{\frac{\pi}{\lambda}}.\frac{1}{2\lambda} = \frac{1}{2\lambda}
得到λ3=12σ2\lambda_3 = -\frac{1}{2\sigma^2}以及C=λπ=12πσ2C = \sqrt{\frac{\lambda}{\pi}} = \sqrt{\frac{1}{2\pi\sigma^2}}
所以
p(x)=12πσ2e(xμ)22σ2p(x) = \sqrt{\frac{1}{2\pi\sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}


XXYY联合熵Joint Entropy)为:
H(X,Y)=xXyYP(x,y)log2[P(x,y)]=EX,Y[logp(x,y)]=x,yp(x,y)logp(x,y){\displaystyle \mathrm {H} (X,Y)=-\sum _{x\in {\mathcal {X}}}\sum _{y\in {\mathcal {Y}}}P(x,y)\log _{2}[P(x,y)]} \\=\mathbb {E} _{X,Y}[-\log p(x,y)]=-\sum _{x,y}p(x,y)\log p(x ,y)\,

积分形式(连续):
h(X,Y)=X,Yf(x,y)logf(x,y)dxdy{\displaystyle h(X,Y)=-\int _{{\mathcal {X}},{\mathcal {Y}}}f(x,y)\log f(x,y)\,dxdy}

多个随机变量:
H(X1,...,Xn)=x1X1...xnXnP(x1,...,xn)log2[P(x1,...,xn)]{\displaystyle \mathrm {H} (X_{1},...,X_{n})=-\sum _{x_{1}\in {\mathcal {X}}_{1}}...\sum _{x_{n}\in {\mathcal {X}}_{n}}P(x_{1},...,x_{n})\log _{2}[P(x_{1},...,x_{n})]}
多个随机变量的积分形式(连续):
h(X1,,Xn)=f(x1,,xn)logf(x1,,xn)dx1dxn{\displaystyle h(X_{1},\ldots ,X_{n})=-\int f(x_{1},\ldots ,x_{n})\log f(x_{1},\ldots ,x_{n})\,dx_{1}\ldots dx_{n}}

XXYY条件熵Conditional Entropy)为:

H(YX) =xX,yYp(x,y)logp(x,y)p(x){\displaystyle \mathrm {H} (Y|X)\ =-\sum _{x\in {\mathcal {X}},y\in {\mathcal {Y}}}p(x,y)\log {\frac {p(x,y)}{p(x)}}}

证明:
H(YX) xXp(x)H(YX=x)=xXp(x)yYp(yx)logp(yx)=xXyYp(x,y)logp(yx)=xX,yYp(x,y)logp(yx)=xX,yYp(x,y)logp(x,y)p(x).=xX,yYp(x,y)logp(x)p(x,y).{\displaystyle {\begin{aligned}\mathrm {H} (Y|X)\ &\equiv \sum _{x\in {\mathcal {X}}}\,p(x)\,\mathrm {H} (Y|X=x)\\&=-\sum _{x\in {\mathcal {X}}}p(x)\sum _{y\in {\mathcal {Y}}}\,p(y|x)\,\log \,p(y|x)\\&=-\sum _{x\in {\mathcal {X}}}\sum _{y\in {\mathcal {Y}}}\,p(x,y)\,\log \,p(y|x)\\&=-\sum _{x\in {\mathcal {X}},y\in {\mathcal {Y}}}p(x,y)\log \,p(y|x)\\&=-\sum _{x\in {\mathcal {X}},y\in {\mathcal {Y}}}p(x,y)\log {\frac {p(x,y)}{p(x)}}.\\&=\sum _{x\in {\mathcal {X}},y\in {\mathcal {Y}}}p(x,y)\log {\frac {p(x)}{p(x,y)}}.\\\end{aligned}}}

积分形式(连续):
h(XY)=X,Yf(x,y)logf(xy)dxdy{\displaystyle h(X|Y)=-\int _{{\mathcal {X}},{\mathcal {Y}}}f(x,y)\log f(x|y)\,dxdy}

根据定义写作:
H(YX)=H(X,Y)H(X){\displaystyle \mathrm {H} (Y|X)\,=\,\mathrm {H} (X,Y)-\mathrm {H} (X)}
一般形式:
H(X1,X2,,Xn)=i=1nH(XiX1,,Xi1){\displaystyle \mathrm {H} (X_{1},X_{2},\ldots ,X_{n})=\sum _{i=1}^{n}\mathrm {H} (X_{i}|X_{1},\ldots ,X_{i-1})}

证明:
H(YX)=xX,yYp(x,y)log(p(x)p(x,y))=xX,yYp(x,y)(log(p(x))log(p(x,y)))=xX,yYp(x,y)log(p(x,y))+xX,yYp(x,y)log(p(x))=H(X,Y)+xXp(x)log(p(x))=H(X,Y)H(X).{\displaystyle {\begin{aligned}\mathrm {H} (Y|X)&=\sum _{x\in {\mathcal {X}},y\in {\mathcal {Y}}}p(x,y)\log \left({\frac {p(x)}{p(x,y)}}\right)\\[4pt]&=\sum _{x\in {\mathcal {X}},y\in {\mathcal {Y}}}p(x,y)(\log(p(x))-\log(p(x,y)))\\[4pt]&=-\sum _{x\in {\mathcal {X}},y\in {\mathcal {Y}}}p(x,y)\log(p(x,y))+\sum _{x\in {\mathcal {X}},y\in {\mathcal {Y}}}{p(x,y)\log(p(x))}\\[4pt]&=\mathrm {H} (X,Y)+\sum _{x\in {\mathcal {X}}}p(x)\log(p(x))\\[4pt]&=\mathrm {H} (X,Y)-\mathrm {H} (X).\end{aligned}}}

互信息(Mutual information

如果变量 𝑋 和 𝑌 互相独立,它们的互信息为零.

I(X;Y)=EX,Y[SI(x,y)]=x,yp(x,y)logp(x,y)p(x)p(y)I(X;Y)=\mathbb {E} _{X,Y}[SI(x,y)]=\sum _{x,y}p(x,y)\log {\frac {p(x,y)}{p(x)\,p(y)}}
其中 SI(Specific mutual Information)是pointwise mutual information

基本性质:
I(X;Y)=H(X)H(XY)=H(Y)H(YX).I(X;Y)=H(X)-H(X|Y) =H(Y)- H(Y|X).\,
对称性:
I(X;Y)=I(Y;X)=H(X)+H(Y)H(X,Y).I(X;Y)=I(Y;X)=H(X)+H(Y)-H(X,Y).\,

互信息可以表示为给定 Y 值的 X 的后验概率分布 与 X 的先验分布之间的平均 Kullback-Leibler 散度:
I(X;Y)=Ep(y)[DKL(p(XY=y)p(X))].I(X;Y)=\mathbb {E} _{p(y)}[D_{\mathrm {KL} }(p(X|Y=y)\|p(X))].
or
I(X;Y)=DKL(p(X,Y)p(X)p(Y)).I(X;Y)=D_{\mathrm {KL} }(p(X,Y)\|p(X)p(Y)).

统计学习方法中讲到 信息增益等价互信息(74 页),而维基百科中信息增益是Kullback-Leibler 散度

交叉熵(Cross Entropy

在给定 分布 𝑝 的情况下,如果 分布 𝑞 和 分布 𝑝 越接近,交叉熵越小;如果 分布 𝑞 和 分布 𝑝 越远,交叉
熵就越大.
H(p,q)=Ep[logq]=xXp(x)logq(x)=H(p)+DKL(pq){\displaystyle H(p,q)=-\operatorname {E} _{p}[\log q]} =-\sum _{x\in {\mathcal {X}}}p(x)\,\log q(x) = H(p)+D_{\mathrm {KL} }(p\|q)

注意与联合熵H(X,Y){H} (X,Y)的区别,联合熵描述一对随机变量平均所需的信息量,交叉熵描述两个分布之间的差异。
也有说交叉熵H(p,q)H(p,q)是不标准的写法,应该是Hq(p)H_q(p) (交叉熵不是对称的,而联合熵是对称的),参见 Difference of notation between cross entropy and joint entropy 以及Relation between cross entropy and joint entropy
应用:一般在多分类问题中使用交叉熵作为损失函数,如:神经网络,逻辑回归

KL 散度(Kullback–Leibler divergence

KL 散度(Kullback-Leibler Divergence),也叫 KL 距离或相对熵(Relative Entropy),是用概率分布 𝑞 来近似 𝑝 时所造成的信息损失量,KL 散度总是大于等于 0 的。可以衡量两个概率分布之间的距离.KL 散度只有当 𝑝 = 𝑞 时,KL(𝑝, 𝑞) = 0.如果两个分布越接近,KL 散度越小;如果两个分布越远,KL 散度就越大.但 KL 散度并不是一个真正的度量或距离,一是 KL 散度不满足距离的对称性,二是 KL 散度不满足距离的三角不等式性质.

DKL(p(X)q(X))=xXp(x)logq(x)xXp(x)logp(x)=xXp(x)logp(x)q(x)=xXp(x)logq(x)p(x).D_{\mathrm {KL} }(p(X)\|q(X))=\sum _{x\in X}-p(x)\log {q(x)}\,-\,\sum _{x\in X}-p(x)\log {p(x)} \\=\sum _{x\in X}p(x)\log {\frac {p(x)}{q(x)}} = -\sum _{x\in X}p(x)\log {\frac {q(x)}{p(x)}}.
也有写作:
KL(p,q),KL(pq),KL(pq),DKL(p,q)KL(p,q) , KL(p|q) , KL(p\|q) , D_{KL}(p,q)

应用:如:变分推断

JS 散度(Jensen-Shannon divergence

JS 散度(Jensen-Shannon Divergence)是一种对称的衡量两个分布相似度的度量方式,是 KL 散度一种改进.但两种散度都存在一个问题,即如果两个分布 𝑝, 𝑞 没有重叠或者重叠非常少时,KL 散度和 JS 散度都很难衡量两个分布的距离.

DJS(PQ)=12DKL(PM)+12DKL(QM){{\rm {D}_{JS}}}(P\parallel Q)={\frac {1}{2}}D_{KL}(P\parallel M)+{\frac {1}{2}}D_{KL}(Q\parallel M)
其中M=12(P+Q)M={\frac {1}{2}}(P+Q), JS 散度也有写作JSD(PQ),JS(PQ),JS(P,Q)JSD(P\|Q), JS(P\|Q) ,JS(P,Q)等。

属于一种统计距离(Statistical distance

统计距离还有Wasserstein 距离Wasserstein distance,也用于衡量两个分布之间的距
离.对于两个分布μ,νpthWasserstein\mu ,\nu,p^{th}-Wasserstein距离定义为

Wp(μ,ν):=(infγΓ(μ,ν)M×Md(x,y)pdγ(x,y))1/p,{\displaystyle W_{p}(\mu ,\nu ):=\left(\inf _{\gamma \in \Gamma (\mu ,\nu )}\int _{M\times M}d(x,y)^{p}\,\mathrm {d} \gamma (x,y)\right)^{1/p},}

Wasserstein 距离相比 KL 散度和 JS 散度的优势在于:即使两个分布没有重叠或者重叠非常少,Wasserstein 距离仍然能反映两个分布的远近.参见pdf443 页附录

参考文献

[5-1] Olshen R A, Quinlan J R. Induction of decision trees. Machine Learning,1986,1(1): 81–106

[5-2] Olshen R A, Quinlan J R. C4. 5: Programs for Machine Learning. Morgan Kaufmann,1992

[5-3] Olshen R A, Breiman L,Friedman J,Stone C. Classification and Regression Trees. Wadsworth,1984

[5-4] Ripley B. Pattern Recognition and Neural Networks. Cambridge UniversityPress,1996

[5-5] Liu B. Web Data Mining: Exploring Hyperlinks,Contents and Usage Data. Springer-Verlag,2006

[5-6] Hyafil L,Rivest R L. Constructing Optimal Binary Decision Trees is NP-complete.Information Processing Letters,1976,5(1): 15–17

[5-7] Hastie T,Tibshirani R,Friedman JH. The Elements of Statistical Learning: DataMining,Inference,and Prediction. New York: Springer-Verlag,2001

[5-8] Yamanishi K. A learning criterion for stochastic rules. Machine Learning,1992

[5-9] Li H,Yamanishi K. Text classification using ESC-based stochastic decision lists.Information Processing & Management,2002,38(3): 343–361

第 6 章 逻辑斯谛回归与最大熵模型

逻辑斯谛回归logistic regression)(也有称 对数几率回归)是统计学习中的经典分类方法。最大熵是概率模型学习的一个准则,将其推广到分类问题得到最大熵模型maximum entropy model)。逻辑斯谛回归模型与最大熵模型都属于对数线性模型(也有称最大熵分类或对数线性分类,所以这里的模型都是分类模型)。

逻辑斯谛回归

一个事件的几率(odds)是指该事件发生的概率与该事件不发生的概率的比值。如果事件发生的概率是 p,那么该事件的几率是p1p\frac{p}{1-p},该事件的对数几率(log odds)或 logit 函数是:
logit(p)=logp1p(6-1)logit(p) = \log\frac{p}{1-p} \label{6-1}\tag{6-1}

sklearn 中代价函数y{1,+1}y \in \{-1,+1\}
P(Y=+1x)=exp(w.x)1+exp(w.x)=σ(w.x)P(Y=1x)=11+exp(w.x)=1σ(w.x)=σ(w.x)P(Y=+1|x) = \frac{\exp{(w.x)}}{1+\exp{(w.x)}} = \sigma{(w.x)} \\ P(Y=-1|x) = \frac{1}{1+\exp{(w.x)}} = 1 - \sigma{(w.x)} = \sigma{(-w.x)}
两个式子合起来就是:σ(yi.w.xi)\sigma{(y_i.w.x_i)}
negative log likelihood:

logi=1Nσ(yi.w.xi)=i=1Nlogσ(yi.w.xi)=i=1Nlog1σ(yi.w.xi)=i=1Nlog1exp(yi.w.xi)1+exp(yi.w.xi)=i=1Nlog(1+1exp(yi.w.xi))=i=1Nlog(1+exp(yi.w.xi))-\log\prod_{i=1}^N \sigma{(y_i.w.x_i)} = \sum_{i=1}^N-\log\sigma{(y_i.w.x_i)}= \sum_{i=1}^N \log\frac{1}{\sigma{(y_i.w.x_i)}} \\= \sum_{i=1}^N \log\frac{1}{\frac{\exp{(y_i.w.x_i)}}{1+\exp{(y_i.w.x_i)}}}\\= \sum_{i=1}^N \log(1+\frac{1}{\exp{(y_i.w.x_i)}})\\= \sum_{i=1}^N \log(1+\exp{(-y_i.w.x_i)})

当然 sklearn 中加入的正则项。

Softmax 回归是 Logistic 回归的多分类情况。
LogisticRegression 就是一个被 logistic 方程归一化后的线性回归。将预测的输出映射到 0,1 之间。

逻辑斯蒂回归模型的思想跟线性回归模型思想不一样,线性回归模型思想是最小化真实值与模型预测值的误差,而逻辑斯蒂回归模型思想就比较狠了,预测值预测对了损失函数就是 0,错了损失就是无穷大,我个人的理解(一般采用的是-log(h(x)) 这是一个凸函数,刚好满足要求)

最大熵模型

熵的概念在统计学习与机器学习中真是很重要,最大熵模型(maximum entropy model)是概率模型学习中一个准则,其思想为:在学习概率模型时,所有可能的模型中熵最大的模型是最好的模型;若概率模型需要满足一些约束,则最大熵原理(Principle of maximum entropy)就是在满足已知约束的条件集合中选择熵最大模型。
最大熵原理指出,对一个随机事件的概率分布进行预测时,预测应当满足全部已知的约束,而对未知的情况不要做任何主观假设。在这种情况下,概率分布最均匀,预测的风险最小,因此得到的概率分布的熵是最大。

均值和方差也被称为一阶矩和二阶矩(二至四阶中心矩被定义为方差(variance)、偏度(skewness)和峰度(kurtosis))
对于连续分布:给定均值和方差,高斯分布的熵最大(也可以说已知均值和方差时,高斯分布随机性最大 证明
对于连续分布:已知区间,连续均匀分布的熵最大
对于连续分布:已知均值(一阶矩),指数分布的熵最大
对于离散分布:离散均匀分布的熵最大(这里在将熵时有证明过)

参考资料

逻辑回归(非常详细)
机器学习实现与分析之四(广义线性模型)

逻辑回归——Logistic 的起源
Logistic 回归的起源(上)
Logistic 回归的起源(中)
Logistic 回归的起源(下)

6.2 Logistic Regression and the Cross Entropy Cost - Logistic regression - y 属于 0 或 1

6.3 Logistic Regression and the Softmax Cost-Logistic regression sklearn 中的代价函数,这里的 y 属于-1 或 1

附加知识

Generalized Linear Models 广义线性模型

Generalized Linear Models (GLM)

Generalized Linear Models

Generalized Linear Models Explained with Examples

Generalized Linear Model Theory - 推荐

Generalized Linear Models

在线性回归模型中的假设中,有两点需要提出:

  1. 假设因变量服从高斯分布:Y=θTx+ξY={{\theta }^{T}}x+\xi,其中误差项ξN(0,σ2)\xi \sim N(0,{{\sigma }^{2}}),那么因变量YN(θTx,σ2)Y\sim N({{\theta }^{T}}x,{{\sigma }^{2}})
  2. 模型预测的输出为E[Y]E[Y],根据Y=θTx+ξY={{\theta }^{T}}x+\xiE[Y]=E[θTx+ξ]=θTxE[Y]=E[{{\theta }^{T}}x+\xi ]={{\theta }^{T}}x,记η=θTx\eta ={{\theta }^{T}}x,则η=E[Y]\eta =E[Y]

广义线性模型可以认为在以上两点假设做了扩展:

  1. 因变量分布不一定是高斯分布,服从一个指数分布族(Exponential family)即可。
  2. 模型预测输出仍然可以认为是E[Y]E[Y](实际上是E[T(Y)]E[T(Y)],许多情况下T(Y)=YT(Y)=Y),但是YY的分布不一定是高斯分布,E[Y]E[Y]η=θTx\eta ={{\theta }^{T}}x也不一定是简单的相等关系,它们的关系用η=g(E[Y])\eta =g(E[Y])描述,称为连接函数(link function),其中η\eta称为自然参数。

由于以上两点的扩展,广义线性模型的应用比基本线性模型广泛许多。对于广义线性这个术语,可以理解为广义体现在因变量的分布形式比较广,只要是一指数分布族即可,而线性则体现在自然参数η=θTx\eta ={{\theta }^{T}}xθ\theta的线性函数。

这一家族中的模型形式基本上都差不多,不同的就是因变量(Y)不同,如果是连续的,就是多重线性回归,如果是伯努利分布,就是 logistic 回归,如果是 poisson 分布,就是 poisson 回归,如果是负二项分布,就是负二项回归,等等。只要注意区分它们的因变量就可以了。logistic 回归的因变量可以是二分类的(二项逻辑回归),也可以是多分类的(多项逻辑回归或者 softmax 回归),但是二分类的更为常用,也更加容易解释。所以实际中最为常用的就是二分类的 logistic 回归。

根据sklearn 中的广义线性回归 Generalized Linear Regression的第二种方式exponential dispersion model (EDM)

其实就是要让真实 y 与预测 y 之间的差异越小越好:
minw12nsamplesid(yi,y^i)+α2w2\min_{w} \frac{1}{2 n_{\text{samples}}} \sum_i d(y_i, \hat{y}_i) + \frac{\alpha}{2} \|w\|_2

假设 y 分别符合下列分布,求真实 y 与预测 y 之间的差异(Deviance)(log 相减不就是两个概率之间的比吗?不就是对数几率(log odds)吗?对数几率为 0 时不就是概率比为 1 吗?不就是差异最小么!):

上述计算都有一个 2 倍,不知道什么意思所以没有写出来。
还有用 link function 解释的,目前不是很明白-参考广义线性模型(GLM)

S 型函数(Logistic & Sigmoid 函数)

Logistic 函数Logistic function)的公式定义:

f(x)=L1+ek(xx0){\displaystyle f(x)={\frac {L}{1+e^{-k(x-x_{0})}}}}
其中LL是最大值,x0x_0是中心点(位置参数),KK是曲线的倾斜度(形状参数,|K|>0 越大,曲线在中心点附近增长越快)。

逻辑斯谛分布Logistic distribution

μ\mu为均值,s 是一个与标准差(standard deviation)成比例的参数

Sigmoid 函数Sigmoid function)的公式定义:
S(x)=11+ex=exex+1=1S(x).{\displaystyle S(x)={\frac {1}{1+e^{-x}}}={\frac {e^{x}}{e^{x}+1}}=1-S(-x).}

Sigmoid 函数是一个有界、可微的实函数,可以将输入压缩到(0,1)区间,经常用作激活函数和概率输出。
Sigmoid 函数对于小于 0 的值是凸的,对于大于 0 的值是凹的。
Sigmoid 函数是一个标准 Logistic 函数(standard logistic function(一般用σ(x)\sigma(x)):K=1,x0=0,L=1K=1,x_0=0,L=1

导数:
S(x)=S(x)(1S(x))S'(x)= S(x)(1-S(x))
推导:
S(x)=11+ex=(1+ex)1S(x)=(1)(1+ex)2ex(1)=(1+ex)2ex=ex(1+ex)2=1+ex1(1+ex)2=1+ex(1+ex)21(1+ex)2=11+ex1(1+ex)2=11+ex(111+ex)S(x) = {\frac {1}{1+e^{-x}}} = (1+e^{-x})^{-1} \\ S'(x)=(-1)*(1+e^{-x})^{-2}*e^{-x}*(-1) \\= (1+e^{-x})^{-2}*e^{-x} \\= \frac{e^{-x}}{(1+e^{-x})^{2}} \\= \frac{1+e^{-x}-1}{(1+e^{-x})^{2}} \\= \frac{1+e^{-x}}{(1+e^{-x})^{2}} - \frac{1}{(1+e^{-x})^{2}} \\=\frac{1}{1+e^{-x}} - \frac{1}{(1+e^{-x})^{2}} \\= \frac{1}{1+e^{-x}}(1-\frac{1}{1+e^{-x}})

双曲正切函数hyperbolic tangent function):
f(x)=12+12tanh(x2),tanh(x)=2f(2x)1.{\displaystyle f(x)={\frac {1}{2}}+{\frac {1}{2}}\tanh \left({\frac {x}{2}}\right),} \\ {\displaystyle \tanh(x)=2f(2x)-1.}

tanh(x)=exexex+ex=ex(1e2x)ex(1+e2x)=f(2x)e2x1+e2x=f(2x)e2x+111+e2x=2f(2x)1.{\displaystyle {\begin{aligned}\tanh(x)&={\frac {e^{x}-e^{-x}}{e^{x}+e^{-x}}}={\frac {e^{x}\cdot \left(1-e^{-2x}\right)}{e^{x}\cdot \left(1+e^{-2x}\right)}}\\&=f(2x)-{\frac {e^{-2x}}{1+e^{-2x}}}=f(2x)-{\frac {e^{-2x}+1-1}{1+e^{-2x}}}=2f(2x)-1.\end{aligned}}}

归一化指数函数(Softmax 函数)

Softmax 函数Softmax function,也称为归一化指数函数 normalized exponential function)的公式定义:
standard (unit) softmax functionσ:RK[0,1]K{\displaystyle \sigma :\mathbb {R} ^{K}\to [0,1]^{K}}

σ(z)i=ezij=1Kezj     for i=1,,K and z=(z1,,zK)RK.{\displaystyle \sigma (\mathbf {z} )_{i}={\frac {e^{z_{i}}}{\sum _{j=1}^{K}e^{z_{j}}}}\ \ \ \ {\text{ for }}i=1,\dotsc ,K{\text{ and }}\mathbf {z} =(z_{1},\dotsc ,z_{K})\in \mathbb {R} ^{K}.}

一般形式:
σ(z)i=eβzij=1Keβzj or σ(z)i=eβzij=1Keβzj for i=1,,K.{\displaystyle \sigma (\mathbf {z} )_{i}={\frac {e^{\beta z_{i}}}{\sum _{j=1}^{K}e^{\beta z_{j}}}}{\text{ or }}\sigma (\mathbf {z} )_{i}={\frac {e^{-\beta z_{i}}}{\sum _{j=1}^{K}e^{-\beta z_{j}}}}{\text{ for }}i=1,\dotsc ,K.}

性质:
σ(z+c)j=ezj+ck=1Kezk+c=ezjeck=1Kezkec=σ(z)j.{\displaystyle \sigma (\mathbf {z} +\mathbf {c} )_{j}={\frac {e^{z_{j}+c}}{\sum _{k=1}^{K}e^{z_{k}+c}}}={\frac {e^{z_{j}}\cdot e^{c}}{\sum _{k=1}^{K}e^{z_{k}}\cdot e^{c}}}=\sigma (\mathbf {z} )_{j}.}

等式左边的c=(c,,c){\displaystyle \mathbf {c} =(c,\dots ,c)}

如果ziz_i都等于一个参数 C 时会发生什么?从理论上输出为(1K,...,1K)RK(\frac{1}{K},...,\frac{1}{K}) \in \mathbb{R}^K,但是从数值计算上说,当 C 很大时eCe^C会发生上溢,当 C 很小时k=1KeC\sum _{k=1}^{K}e^C会发生下溢,这时我们就可以利用上述性质,将z\mathbf {z}减去maxizi\max_i {z_i},那么最大值就是 0,排除了上溢的可能,同样的分母至少有一个为 1 的项,排除了因为分母下溢而导致被 0 除的可能性。
σ(z)i=ezij=1Kezj=e(zizmax)j=1Ke(zjzmax)=e(zizmax)1+j=1,jmaxKe(zjzmax)\sigma (\mathbf {z} )_{i}={\frac {e^{z_{i}}}{\sum _{j=1}^{K}e^{z_{j}}}} = {\frac {e^{(z_{i}-z_{max})}}{\sum _{j=1}^{K}e^{(z_{j}-z_{max})}}} = {\frac {e^{(z_{i}-z_{max})}}{1+\sum _{j=1,j\neq max}^{K}e^{(z_{j}-z_{max})}}}

log softmax函数在深度学习中也经常遇见,其实就是求完 softmax,再对其求 log,如果直接计算可能会出现问题(当 softmax 很小时,log 会得到-\infty),这时我就要推导出 log softmax 的表达式:
log(σ(z))i=loge(zizmax)j=1Ke(zjzmax)=loge(zizmax)logj=1Ke(zjzmax)=(zizmax)logj=1Ke(zjzmax)\log(\sigma (\mathbf {z} ))_{i}=\log{\frac {e^{(z_{i}-z_{max})}}{\sum _{j=1}^{K}e^{(z_{j}-z_{max})}}} \\= \log e^{(z_{i}-z_{max})} - \log {\sum _{j=1}^{K}e^{(z_{j}-z_{max})}} \\= (z_{i}-z_{max})- \log {\sum _{j=1}^{K}e^{(z_{j}-z_{max})}}

j=1Ke(zjzmax){\sum _{j=1}^{K}e^{(z_{j}-z_{max})}}是大于等于 1 的,并且不会大的离谱,所以不会出问题。

negative log-likelihood(NLL),likelihood 是一个概率(softmax 也是概率),所以 log-likelihood 小于 0,negative log-likelihood 则大于 0,这样就可以最小化 negative log-likelihood 了

参考文献

[6-1] Berger A,Della Pietra SD,Pietra VD. A maximum entropy approach to naturallanguage processing. Computational Linguistics,1996,22(1),39–71

[6-2] Berger A. The Improved Iterative Scaling Algorithm: A Gentle Introduction.http://www.cs.cmu.edu/ afs/cs/user/aberger/www/ps/scaling.ps

[6-3] Hastie T,Tibshirani R,Friedman J. The Elements of Statistical Learning: DataMining,Inference,and Prediction. Springer-Verlag. 2001(中译本:统计学习基础——数据挖掘、推理与预测。范明,柴玉梅,昝红英等译。北京:电子工业出版社,2004)

[6-4] Mitchell TM. Machine Learning. McGraw-Hill Companies,Inc. 1997(中译本:机器学习。北京:机械工业出版社,2003)

[6-5] Collins M,Schapire RE,Singer Y. Logistic Regression,AdaBoost and BregmanDistances. Machine Learning Journal,2004

[6-6] Canu S,Smola AJ. Kernel method and exponential family. Neurocomputing,2005,69:714–720

第 7 章 支持向量机

支持向量机support vector machine,SVM)是一种二类分类模型。它的基本模型是定义在特征空间上的间隔最大的线性分类器,间隔最大使它有别于感知机;支持向量机还包括核技巧,这使它成为实质上的非线性分类器。支持向量机的学习策略就是间隔最大化,可形式化为一个求解凸二次规划(convex quadratic programming)的问题,也等价于正则化的合页损失函数的最小化问题。支持向量机的学习算法是求解凸二次规划的最优化算法。

支持向量机学习方法包含构建由简至繁的模型:线性可分支持向量机(linear supportvector machine in linearly separable case)、线性支持向量机(linear support vectormachine)及非线性支持向量机(non-linear support vector machine)。简单模型是复杂模型的基础,也是复杂模型的特殊情况。当训练数据线性可分时,通过硬间隔最大化(hardmargin maximization),学习一个线性的分类器,即线性可分支持向量机,又称为硬间隔支持向量机;当训练数据近似线性可分时,通过软间隔最大化(soft marginmaximization),也学习一个线性的分类器,即线性支持向量机,又称为软间隔支持向量机;当训练数据线性不可分时,通过使用核技巧(kernel trick)及软间隔最大化,学习非线性支持向量机

当输入空间为欧氏空间或离散集合、特征空间为希尔伯特空间时,核函数(kernelfunction)表示将输入从输入空间映射到特征空间得到的特征向量之间的内积。通过使用核函数可以学习非线性支持向量机,等价于隐式地在高维的特征空间中学习线性支持向量机。这样的方法称为核技巧。核方法(kernel method)是比支持向量机更为一般的机器学习方法。

SVM 有三宝:间隔、对偶、核技巧

函数间隔:γ^=y(w.x+b),y{1,+1}\hat\gamma = y(w.x + b), y\in \{-1,+1\}
几何间隔:γ=1ww.x+b=γ^w\gamma = \frac{1}{\|w\|}|w.x + b| = \frac{\hat\gamma}{\|w\|}

w.x+b=0w.x + b=0

  1. 最大间隔(硬间隔原始问题-凸二次规划)
    maxw,bγ^ws.t.yi(w.xi+b)γ^\begin{aligned} &\max_{w,b} &\frac{\hat\gamma}{\|w\|} \\ &\text{s.t.} &y_i(w.x_i+b) \geq \hat\gamma\end{aligned}
    函数间隔γ^\hat\gamma的大小是可以变的,我们让其等于 1,那么将上述问题改写下:
    minw,b12w2s.t.yi(w.xi+b)1\begin{aligned} &\min_{w,b} &\frac{1}{2}\|w\|^2 \\ &\text{s.t.} &y_i(w.x_i+b) \geq 1\end{aligned}
    这不就是一个标准的凸二次规划问题么!(1w\frac{1}{\|w\|}w=0\|w\| = 0处不可微,arg max1warg min12w2\argmax \frac{1}{\|w\|} 与 \argmin \frac{1}{2}\|w\|^2等价)
    如果数据集线性可分,那么最大间隔分离超平面存在且唯一,具体证明就不证了(见统计学习方法 117 页)。
    在线性可分的情况下,训练集中的样本点与分离超平面距离最近的点称为支持向量,也就是满足yi(w.xi+b)=1y_i(w.x_i +b ) =1的点。
    而超平面w.xi+b=+1,w.xi+b=1w.x_i +b = +1,w.x_i +b = -1称为间隔边界,两个间隔边界之间的距离称为间隔(margin),间隔大小为2w\frac{2}{\|w\|}
    在决定分离超平面时只有支持向量起作用,而其它样本点并不起作用,所以该模型叫做支持向量机。

    • 拉格朗日函数(求最小)
      L(w,b,α)=12w2+i=1Nαi(1yi(w.xi+b))L(w,b,\alpha) = \frac{1}{2}\|w\|^2 + \sum_{i=1}^N \alpha_i(1-y_i(w.x_i+b))
      minw,bmaxαL(w,b,α)\min_{w,b} \max_{\alpha} L(w,b,\alpha)

    • 拉格朗日对偶函数(求最大)
      maxαg(α)=maxαinfw,bL(w,b,α)\max_{\alpha} g(\alpha) = \max_{\alpha} \inf_{w,b} L(w,b,\alpha)

  2. 带正则项的合页损失函数(软间隔原始问题-凸二次规划)
    minw,bi=1Nmax(0,1yi(w.xi+b))hinge loss function+λw2正则化项\min_{w,b} \underbrace{\sum_{i=1}^N\max(0,1-y_i(w.x_i+b))}_{\text{hinge loss function}} + \underbrace{\lambda\|w\|^2}_{\text{正则化项}}
    等价软间隔最大化的优化问题:
    minw,b,ξi=1nξi+λw2subject to yi(wTxi+b)1ξiξi0,for all i.\begin{aligned} &\min_{w,b,\xi} & \sum _{i=1}^{n}\xi _{i}+\lambda \|\mathbf {w} \|^{2} \\ &\displaystyle {\text{subject to }} & y_{i}(\mathbf {w} ^{T}\mathbf {x} _{i}+b)\geq 1-\xi _{i}\\ &&\xi _{i}\geq 0,\,{\text{for all }}i.\end{aligned}
    其中ξi=max(0,1yi(wTxi+b)){\displaystyle \xi _{i}=\max \left(0,1-y_{i}(\mathbf {w} ^{T}\mathbf {x} _{i}+b)\right)}是松弛变量(松弛方法Relaxation有很多,下面只是一种而已)(跟下面松弛变量定义不同,只有统计学习方法中有说到,其它地方没有找到)。
    yi(wTxi+b)1ξiy_{i}(\mathbf {w} ^{T}\mathbf {x} _{i}+b)\geq 1-\xi _{i}相当于分类点可以处于间隔中(不太准确,处于间隔边界的一侧),对于软间隔支持向量机中的支持向量包含了间隔中的向量。

    • 拉格朗日函数(求最小)
      L(w,b,ξ,α,μ)=λw2+i=1Nξi+i=1Nαi(1ξiyi(w.xi+b))+i=1N(μi(ξi))L(w,b,\xi,\alpha,\mu) = \lambda\|w\|^2 + \sum_{i=1}^N \xi_{i} + \sum_{i=1}^N \alpha_i(1-\xi_{i}-y_i(w.x_i+b)) + \sum_{i=1}^N(\mu_i*(-\xi_{i}))
      minw,b,ξmaxα,μL(w,b,ξ,α,μ)\min_{w,b,\xi} \max_{\alpha,\mu} L(w,b,\xi,\alpha,\mu)

    • 拉格朗日对偶函数(求最大)
      maxα,μg(α,μ)=maxα,μinfw,b,ξL(w,b,ξ,α,μ)\max_{\alpha,\mu} g(\alpha,\mu) = \max_{\alpha,\mu} \inf_{w,b,\xi} L(w,b,\xi,\alpha,\mu)

核技巧
首先写出了原形式的对偶形式,然后把xϕ(x)x用\phi (x)代替,最终发现根本不需要知道ϕ(x)\phi (x),只需要核函数就行了,具体证明就不证了,很简单,上面也有介绍了核技巧知识。

书中还介绍了原形式的对偶形式(区别于拉格朗日对偶),也就是等价形式(感知机中 2.3.3 节 44 页 也是等价的意思),这两个地方的等价都是经过基本推导,求出 w 参数,然后对原问题进行了替换。

附加知识

拉格朗日乘数法

最优化:建模、算法与理论/最优化计算方法
Stephen Boyd 的 Convex Optimization - 凸优化
Nonlinear Programming by Dimitri P. Bertsekas - 非线性规划

凸优化Convex optimization):
凸优化问题是目标函数为凸函数,可行集为凸集的优化问题。
标准形式:
minimizexf(x)subject togi(x)0,i=1,,mhi(x)=0,i=1,,p,{\displaystyle {\begin{aligned}&{\underset {\mathbf {x} }{\operatorname {minimize} }}&&f(\mathbf {x} )\\&\operatorname {subject\ to} &&g_{i}(\mathbf {x} )\leq 0,\quad i=1,\dots ,m\\&&&h_{i}(\mathbf {x} )=0,\quad i=1,\dots ,p,\end{aligned}}}
其中xRn\mathbf {x} \in \mathbb {R} ^{n}为优化变量,目标函数(objective function)f:DRnR{\displaystyle f:{\mathcal {D}}\subseteq \mathbb {R} ^{n}\to \mathbb {R} }是凸的,不等式约束gi:RnR{\displaystyle g_{i}:\mathbb {R} ^{n}\to \mathbb {R} }也是凸的,等式约束hi:RnR{\displaystyle h_{i}:\mathbb {R} ^{n}\to \mathbb {R} }仿射affine)的

二次约束二次规划Quadratically constrained quadratic program):
minimize12xTP0x+q0Txsubject to12xTPix+qiTx+ri0for i=1,,m,Ax=b,{\begin{aligned}&{\text{minimize}}&&{\tfrac 12}x^{{\mathrm {T}}}P_{0}x+q_{0}^{{\mathrm {T}}}x\\&{\text{subject to}}&&{\tfrac 12}x^{{\mathrm {T}}}P_{i}x+q_{i}^{{\mathrm {T}}}x+r_{i}\leq 0\quad {\text{for }}i=1,\dots ,m,\\&&&Ax=b,\end{aligned}}
其中P0以及P1,..,PmRn×nP_0以及P_1,..,P_m \in \mathbb{R}^{n \times n},xRn\mathbf {x} \in \mathbb {R} ^{n}为优化变量
如果P0以及P1,..,PmRn×nP_0以及P_1,..,P_m \in \mathbb{R}^{n \times n}是半正定矩阵,那么问题是凸的,如果P1,..,PmP_1,..,P_m为 0,那么约束是线性的,就是二次规划Quadratic programming),即目标函数是二次的,不等式以及等式约束也是线性的;二次规划的前提下,如果P0P_0是半正定矩阵那么就是凸二次规划;如果P0P_0为 0,就是不标准的线性规划Linear programming),即目标函数是线性的,不等式以及等式约束也是线性的。
标准的线性规划:即目标函数是线性的,非负约束(优化变量是非负的),等式约束也是线性的。

线性规划解法有单纯形法等。其它规划的优化算法看这里:内点法,单纯形法等;有线性规划自然也有动态规划Dynamic programming

最小二乘不就是凸二次规划么yf(x)2,f(x)=Ax+c\|y-f(x)\|^2,f(x) = Ax+c):
minimizef(x)=12Axb2=12(Axb)T(Axb)=12(xTATAxxTATbbTAx+btb)=12(xTATAx2xTATb+btb)\text{minimize} f(x) =\frac{1}{2} \|Ax-b\|^2 = \frac{1}{2}(Ax-b)^T(Ax-b) \\= \frac{1}{2}(x^TA^TAx - x^TA^Tb -b^TAx + b^tb) \\= \frac{1}{2}(x^TA^TAx - 2x^TA^Tb + b^tb)
其中ARm×nxRn×1bRm×1A \in \mathbb{R}^{m \times n},x \in \mathbb{R}^{n \times 1},b \in \mathbb{R}^{m \times 1},所以xTATbbTAxx^TA^Tb 和 b^TAx都是相等的实数,btbb^tb也是实数ATAA^TA是半正定矩阵

L1 和 L2 回归也能转化称相应的优化问题,参考:L1L2 正则化和凸优化

共轭函数conjugate function):
设函数f:RnRf:\mathbb{R}^n \to \mathbb{R},定义 f 的共轭函数f:RnRf^*:\mathbb{R}^n \to \mathbb{R}为:
f(y)=supxdomf(yTxf(x))f^*(y) = \sup_{x \in \mathrm{dom} f} (y^Tx - f(x))
共轭函数一定是凸的,supremum 为上界, infimum 为下界。

拉格朗日乘数法Lagrange multiplier):
根据上面标准形式的优化问题,我们来构造一个拉格朗日函数:
L(x,λ,ν)=f(x)+i=1mλigi(x)+i=1pνihi(x)L(x,\lambda,\nu) = f(x) +\sum _{i=1}^{m}\lambda _{i}g_{i}(x)+\sum _{ i=1}^{p}\nu _{i}h_{i}(x)
其中λi,νi\lambda _{i},\nu _{i}分别是不等式和等式对应的 Lagrange 乘子,当然如果用向量λ,ν\lambda,\nu表示,称为原问题的 Lagrange 乘子向量对偶变量

对偶函数Dual function):
g(λ,ν)=infxL(x,λ,ν)g(\lambda,\nu) = \inf_{x} L(x,\lambda,\nu)
对偶函数一定是凹的,又称Lagrange 对偶函数。对偶函数g(λ,ν)pg(\lambda,\nu) \leq p^{\star}, pp^{\star}是原问题的最优值。Lagrange 对偶问题的最优值用dd^{\star}表示,则dpd^{\star}\leq p^{\star},这个性质称为弱对偶性Weak Duality),pdp^{\star}- d^{\star}称为原问题的最优对偶间隙Duality gap),当d=pd^{\star}= p^{\star}时称为强对偶性Strong Duality)。

所以当强对偶性成立时,拉格朗日函数的最小等价于对偶函数的最大值minxmaxλ,νL(x,λ,ν)    maxλ,νg(λ,ν)\min_{x} \max_{\lambda,\nu}L(x,\lambda,\nu) \iff \max_{\lambda,\nu} g(\lambda,\nu)我们一般使用拉格朗日函数求最小化时的参数xx,然后带入其中,利用对偶函数求其最大时的 Lagrange 乘子,即λ,ν\lambda,\nu
这里的maxλ,νL(x,λ,ν)\max_{\lambda,\nu}L(x,\lambda,\nu)是为了让约束起作用(描述不是很准确,其实就是目标函数和约束的交点,参见附录 C.3.2)。

当强对偶性成立时,那么x,λ,νx^{\star},\lambda^{\star},\nu^{\star}分别是原问题和对偶问题的最优解的充分必要条件是满足下面的KKT 条件Karush–Kuhn–Tucker conditions):

xL(x,λ,ν)=0λigi(x)=0,i=1,2,...,m(Complementary slackness)λi0,i=1,2,...,m(Dual feasibility)gi(x)0,i=1,2,...,m(Primal feasibility)hi(x)=0,i=1,2,...,p(Primal feasibility)\nabla_x L(x^{\star},\lambda^{\star},\nu^{\star}) = 0 \\ \lambda_i^{\star}g_i(x^{\star}) = 0,i=1,2,...,m \quad (\text{Complementary slackness})\\ \lambda_i^{\star} \geq 0,i=1,2,...,m \quad (\text{Dual feasibility})\\ g_i(x^{\star}) \leq 0,i=1,2,...,m \quad (\text{Primal feasibility})\\ h_i(x^{\star}) = 0,i=1,2,...,p \quad (\text{Primal feasibility})

可以看到我们只对 x 参数求导,没有对 Lagrange 乘子求导;

其中互补松弛(Complementary slackness)条件用i=0mλigi(x)=0\sum_{i=0}^m\lambda_i^{\star}g_i(x^{\star})=0可能很合理,如果最优解xx^{\star}出现在不等式约束的边界上gi(x)=0g_i(x) = 0,则λi>0\lambda_i^{\star} > 0;如果最优解xx^{\star}出现在不等式约束的内部gi(x)<0g_i(x) < 0,则λi=0\lambda_i^{\star} = 0互补松弛条件说明当最优解出现在不等式约束的内部,则约束失效,所以λi0,i=1,2,...,m\lambda_i^{\star} \geq 0,i=1,2,...,m表示对偶可行性(Dual feasibility)。

如何将不标准的优化问题转换称标准的优化问题(线性规划):参考线性规划问题
A. 如何将不等式约束变成等式约束:

  1. aTxba^Tx \leq b
    只需要加上松弛变量(Slack variable),松弛变量是添加到不等式约束以将其转换为等式的变量,松弛变量特别用于线性规划。松弛变量不能取负值,因为单纯形算法要求它们为正值或零。
    aTx+s=bs0a^Tx +s = b \\ s \geq 0

  2. aTxba^Tx \geq b
    只需要减去剩余变量(surplus variable),剩余变量不能取负值。
    aTxe=be0a^Tx - e = b \\ e \geq 0

B. 无约束变量变成有非负约束变量:
z1=z1+z1z1=z1++z1z1+,z1=0z1+,z10{\begin{aligned}&z_{1}=z_{1}^{+}-z_{1}^{-} \\& |z_{1}| = z_{1}^{+}+\,z_{1}^{-} \\&\braket{z_{1}^{+},z_{1}^{-}} = 0 \\&z_{1}^{+},\,z_{1}^{-} \geq 0\end{aligned}}

如:
a. 利用上述第一条性质
5=505=055 = 5-0 \\ -5 = 0-5
or

(1234)=(1200)(0034)\begin{pmatrix} 1 \\ 2 \\ -3 \\ -4 \end{pmatrix} = \begin{pmatrix} 1 \\ 2 \\ 0 \\ 0 \end{pmatrix} - \begin{pmatrix} 0 \\ 0 \\ 3 \\ 4 \end{pmatrix}

b. 目标函数有带绝对值的(第二条性质)
maxx等价:maxx++xs.t.x+0x0\begin{aligned} \max |x| \\ 等价:\max x^+ + x^- \\ s.t. \quad x^+ \geq 0 \\ x^- \geq 0\end{aligned}

参考文献

[7-1] Cortes C,Vapnik V. Support-vector networks. Machine Learning,1995,20

[7-2] Boser BE,Guyon IM,Vapnik VN. A training algorithm for optimal margin classifiers.In: Haussler D,ed. Proc of the 5th Annual ACM Workshop on COLT. Pittsburgh,PA,1992,144–152

[7-3] Drucker H,Burges CJC,Kaufman L,Smola A,Vapnik V. Support vector regressionmachines. In: Advances in Neural Information Processing Systems 9,NIPS 1996. MITPress,155–161

[7-4] Vapnik Vladimir N. The Nature of Statistical Learning Theory. Berlin: Springer-Verlag,1995(中译本:张学工,译。统计学习理论的本质。北京:清华大学出版社,2000)

[7-5] Platt JC. Fast training of support vector machines using sequential minimaloptimization. Microsoft Research

[7-6] Weston JAE,Watkins C. Support vector machines for multi-class pattern recognition. In: Proceedings of the 7th European Symposium on Articial Neural Networks. 1999

[7-7] Crammer K,Singer Y. On the algorithmic implementation of multiclass kernel-basedmachines. Journal of Machine Learning Research,2001,2(Dec): 265–292

[7-8] Tsochantaridis I,Joachims T,Hofmann T,Altun Y. Large margin methods forstructured and interdependent output variables. JMLR,2005,6: 1453–1484

[7-9] Burges JC. A tutorial on support vector machines for pattern recognition. BellLaboratories,Lucent Technologies. 1997

[7-10] Cristianini N,Shawe-Taylor J. An Introduction to Support Vector Machines andOthre KernerBased Learning Methods. Cambridge University Press,2000(中译本:李国正,等译。支持向量机导论。北京:电子工业出版社,2004)

[7-11] 邓乃扬,田英杰。数据挖掘中的新方法——支持向量机。北京:科学出版社,2004

[7-12] 邓乃扬,田英杰。支持向量机——理论,算法与拓展。北京:科学出版社,2009

[7-13] Scholkpf B,Smola AJ. Learning with Kernels: Support VectorMachines,Regularization,Optimization,and Beyond. MIT Press,2002

[7-14] Herbrich R. Learning Kernel Classifiers,Theory and Algorithms. The MITPress,2002

[7-15] Hofmann T,Scholkopf B,Smola AJ. Kernel methods in machine learning. The Annalsof Statistics,2008,36(3): 1171–1220

第 8 章 提升方法

集成学习Ensemble Learning)也叫集成方法(Ensemble methods)是一种将多种学习算法组合在一起以取得更好表现的一种方法。

利用成员模型的多样性来纠正某些成员模型错误的能力,常用的多样性技术有:

弱分类器:比随机猜测略好,如二分类中,准确率大于 0.5 就可以了(如 0.51)。

参见的集成学习类型:

模型的组合(blending)方法:

aggregation type blending learning
uniform voting/averaging Bagging
non-uniform linear AdaBoost
conditional stacking Decision Tree

参考

AdaBoostAdaBoost是 Adaptive Boosting 的缩写

训练集T={(x1,y1),(x2,y2),...,(xN,yN)},xiRn,yi{1,+1}T = \{(x_1,y_1),(x_2,y_2),...,(x_N,y_N)\}, x_i \in \mathbb{R}^n,y_i\in \{-1,+1\}

Gradient Tree Boosting:梯度提升(Gradient boosting)是一种用于回归、分类和其他任务的机器学习技术,它以弱预测模型(通常是决策树)的集合形式生成预测模型。当决策树是弱学习器时,产生的算法称为梯度提升树(Gradient Tree Boosting or Gradient boosted trees or Gradient Boosted Decision Trees(GBDT)),通常优于随机森林(Random forest它是 Bagging 算法的进化版,也就是说,它的思想仍然是 bagging,但是进行了独有的改进。)

随机森林与一般的 bagging 相比:(参考:Bagging 与随机森林算法原理小结机器学习算法系列(五):bagging 与随机森林对比及随机森林模型参数介绍

因为 Bagging 使用的有放回采样,所以 BaggingClassifier or RandomForestClassifier 都具有 oobscore属性:约有 37%(limn(11/n)n=1/e\lim_{n \to -\infty}(1-1/n)^n=1/e)的样本没有用来训练,这一部分称为 out-of-bag(oob),因为模型没有见过这部分样本,所以可以拿来当验证集合,而不需要再划分验证集或者交叉验证了。 比如我们计算 accuracyscore 时,也可以看下 oob_score的情况

参考文献

[8-1] Freund Y,Schapire RE. A short introduction to boosting. Journal of Japanese Societyfor Artificial Intelligence,1999,14(5): 771–780

[8-2] Hastie T,Tibshirani R,Friedman J. The Elements of Statistical Learning: DataMining,Inference,and Prediction. Springer-Verlag,2001(中译本:统计学习基础——数据挖掘、推理与预测。范明,柴玉梅,昝红英,等译。北京:电子工业出版社,2004)

[8-3] Valiant LG. A theory of the learnable. Communications of the ACM,1984,27(11):1134–1142

[8-4] Schapire R. The strength of weak learnability. Machine Learning,1990,5(2): 197–227

[8-5] Freund Y,Schapire RE. A decision-theoretic generalization of on-line learning and anapplication to boosting. Computational Learning Theory. Lecture Notes in ComputerScience,Vol. 904,1995,23–37 (55, 119-139 (1997)

[8-6] Friedman J,Hastie T,Tibshirani R. Additive logistic regression: a statistical view ofboosting(with discussions). Annals of Statistics,2000,28: 337–407

[8-7] Friedman J. Greedy function approximation: a gradient boosting machine. Annals ofStatistics,2001,29(5)

[8-8] Schapire RE,Singer Y. Improved boosting algorithms using confidence-ratedpredictions. Machine Learning,1999,37(3): 297–336

[8-9] Collins M,Schapire R E,Singer Y. Logistic regression,AdaBoost and Bregmandistances. Machine Learning Journal,2004

第 9 章 EM 算法及其推广

EM 算法(Expectation–maximization algorithm)是一种迭代算法,用于含有隐变量(hidden(unseen or unmeasurable) variable or Latent variable)的概率模型参数的极大似然估计,或极大后验概率估计。EM 算法的每次迭代由两步组成:E 步,求期望(expectation);M 步,求极大(maximization)。所以这一算法称为期望极大算法(expectation maximization algorithm),简称 EM 算法。

概率模型有时既含有观测变量(observable variable),又含有隐变量或潜在变量(latent variable)。如果概率模型的变量都是观测变量,那么给定数据,可以直接用极大似然估计法,或贝叶斯估计法估计模型参数。但是,当模型含有隐变量时,就不能简单地使用这些估计方法。EM 算法就是含有隐变量的概率模型参数的极大似然估计法,或极大后验概率估计法。

本章首先叙述 EM 算法,然后讨论 EM 算法的收敛性;作为 EM 算法的应用,介绍高斯混合模型的学习;最后叙述 EM 算法的推广——GEM 算法(generalized expectation maximization (GEM) algorithm,广义 EM 算法)。

EM 算法是一个优化算法,不是一个统计学习模型。
EM 算法优点:不需要调参数,没有超参数;编程简单,只需要迭代;理论优美,收敛性。
EM 算法的推广:F 函数(F function) 的极大-极大算法 F-MM or MM(maximization maximization algorithm),MCEM,VBEM or VEM,GEM

观测数据X={xi}i=1NX=\{x_i\}_{i=1}^N ,对应的隐含(隐藏)数据Z={zi}i=1NZ=\{z_i\}_{i=1}^N,模型参数θ\theta,完全数据T={(x1,z1),...,(xN,zN)}T=\{(x_1,z_1),...,(x_N,z_N)\}

注意下列公式中|也有用;;表示的,是一个意思,因为θ\theta是一个未知的参数,最开始时会初始化θ(0)\theta^{(0)}
公式中的大 P 和小 p 以及大 X 和小 x 没有统一,不是很严谨

如果不考虑隐藏数据,我们就可以直接使用极大似然估计的方法估计出参数 θ\theta :
θMLE=argmaxθlogp(xθ)=argmaxθiNlogp(xiθ)\theta_{MLE} = \arg\max\limits_\theta\log p(x|\theta) = \arg\max\limits_\theta\sum_i^N\log p(x_i|\theta)

但由于隐藏数据的存在, 我们有 x 的边际似然函数(Marginal Likelihood),在贝叶斯统计的上下文中,边际似然也称为证据(Evidence)
p(xθ)=zp(x,zθ)=Zp(x,zθ)dzp(x|\theta) = \sum_{z} p(x, z|\theta) = \int_Z p(x, z|\theta)dz
将 x 的边际似然函数带入极大似然估计中,在 log 里面会出现积分(求和)符号,导致对似然函数的求导变得困难,无法求解。对于这种无法直接求解的问题,我们通常会采用迭代求解的策略,一步一步逼近最终的结果,在 EM 算法中就是 E 步和 M 步的交替进行,直至收敛。

下面来介绍EM 算法

  1. 随机化参数θ(0)\theta^{(0)}的初始值;

  2. 假设在第 tt 次迭代后,参数的估计值为 θ(t)\theta^{(t)} ,对于第 t+1t+1 次迭代,具体分为两步:
    a. E-step:求期望
    Q 函数的定义:
    Q(θ,θ(t))=ZP(ZX,θ(t))logP(X,Zθ)dZ=EZX,θ(t)[logP(X,Zθ)]Q(\theta,\theta^{(t)}) = \int_Z P(Z|X,\theta^{(t)}) \log P(X,Z|\theta) dZ \\ =\mathbb{E}_{Z|X,\theta^{(t)}}[\log P(X,Z|\theta)]
    全数据的对数似然函数logp(x,zθ)\log p(x, z|\theta)关于在给定观测数据XX和当前参数θ(t)\theta^{(t)}下对未观测数据ZZ的条件概率分布p(zx,θ(t))p(z|x,\theta^{(t)})的期望称为 Q 函数

    注意 Q 函数的第一个参数是自变量,第二个参数是已知的(上一步求得的)

    b. M-step: 最大化Q(θ,θ(t))Q(\theta,\theta^{(t)}) 并求解θ(t+1)\theta^{(t+1)}
    θ(t+1)=argmaxθQ(θ,θ(t))\theta^{(t+1)} = \arg\max\limits_\theta Q(\theta, \theta^{(t)})

  3. 重复 2,直到收敛。
    一般判断收敛有两种方法:1. 判断参数是否收敛θ(t+1)θ(t)ε\theta^{(t+1)}-\theta^{(t)} \leq \varepsilon;2. 判断函数值是否收敛Q(θ(t+1),θ(t))Q(θ(t),θ(t))εQ(\theta^{(t+1)},\theta^{(t)})-Q(\theta^{(t)},\theta^{(t)}) \leq \varepsilon

EM 算法的导出
为什么 EM 算法能近似实现对观测数据的极大似然估计?
书中的推导我这里就不重复了(书中用的 Jensen 不等式推导),这里介绍变分法/ELBO+KL
P(Xθ)=P(X,Zθ)P(ZX,θ)    logP(Xθ)=logP(X,Zθ)logP(ZX,θ)P(X|\theta) = \frac{P(X,Z|\theta)}{P(Z|X, \theta)} \implies \log P(X|\theta) = \log P(X,Z|\theta) - \log P(Z|X,\theta)
根据变分推断的思想:寻找一个简单分布q(z)q(z)来近似条件概率密度p(zx)p(z|x)
logP(Xθ)=logP(X,Zθ)q(Z)logP(ZX,θ)q(Z)\log P(X|\theta) = \log \frac{P(X,Z|\theta)}{q(Z)} - \log \frac{P(Z|X,\theta)}{q(Z)}
然后两边同时求关于变量 ZZ 的期望
EZ[logP(Xθ)]=E_Z[logP(X,Zθ)q(Z)]E_Z[logP(ZX,θ)q(Z)]\mathbb{E}_Z[\log P(X|\theta)] = \mathbb{E}\_Z[\log \frac{P(X,Z|\theta)}{q(Z)}] - \mathbb{E}\_Z[\log \frac{P(Z|X,\theta)}{q(Z)}]
将期望写成积分的形式
Zq(Z)logP(Xθ)dZ=Zq(Z)logP(X,Zθ)q(Z)dZZq(Z)logP(ZX,θ)q(Z)dZ\int_Z q(Z)\log P(X|\theta)dZ = \int_Zq(Z)\log \frac{P(X,Z|\theta)}{q(Z)}dZ - \int_Zq(Z)\log \frac{P(Z|X,\theta)}{q(Z)}dZ
等式左边和ZZ无关(Zq(Z)dZ=1\int_Zq(Z)dZ = 1),所以
logP(Xθ)=Zq(Z)logP(X,Zθ)q(Z)dZZq(Z)logP(ZX,θ)q(Z)dZ=Zq(Z)logP(X,Zθ)q(Z)dZ+DKL(q(Z)P(ZX,θ))=ELBO+DKL(q(Z)P(ZX,θ))\log P(X|\theta) = \int_Zq(Z)\log \frac{P(X,Z|\theta)}{q(Z)}dZ - \int_Zq(Z)\log \frac{P(Z|X,\theta)}{q(Z)}dZ \\ = \int_Zq(Z)\log \frac{P(X,Z|\theta)}{q(Z)}dZ + D_{KL}(q(Z)||P(Z|X,\theta)) \\= ELBO + D*{KL}(q(Z)||P(Z|X,\theta))
我们直接令DKL=0,q(Z)=P(ZX,θ(t))D*{KL} = 0, 即 q(Z)=P(Z|X,\theta^{(t)}),然后最大化ELBO
θ^=arg maxθZq(Z)logP(X,Zθ)q(Z)dZ=arg maxθZP(ZX,θ(t))logP(X,Zθ)P(ZX,θ(t))dZ\hat{\theta} = \argmax_{\theta}\int_Zq(Z)\log \frac{P(X,Z|\theta)}{q(Z)}dZ \\ = \argmax_{\theta}\int_Z P(Z|X,\theta^{(t)})\log \frac{P(X,Z|\theta)}{P(Z|X,\theta^{(t)})}dZ
θ(t)\theta^{(t)}是上一步求出的,可以看作已知的参数
θ^=arg maxθZP(ZX,θ(t))logP(X,Zθ)dZZP(ZX,θ(t))logP(ZX,θ(t))dZ=arg maxθZP(ZX,θ(t))logP(X,Zθ)dZC=arg maxθZP(ZX,θ(t))logP(X,Zθ)dZ=arg maxθQ(θ,θ(t))\hat{\theta} = \argmax_{\theta} \int_Z P(Z|X,\theta^{(t)})\log P(X,Z|\theta)dZ - \int_Z P(Z|X,\theta^{(t)})\log {P(Z|X,\theta^{(t)})}dZ \\ = \argmax_{\theta} \int_Z P(Z|X,\theta^{(t)})\log P(X,Z|\theta)dZ - C \\ = \argmax_{\theta} \int_Z P(Z|X,\theta^{(t)})\log P(X,Z|\theta)dZ \\= \argmax_{\theta} Q(\theta,\theta^{(t)})

参考EM 算法 11.2.2.1 节
参考深入理解 EM 算法(ELBO+KL 形式)
参考深入理解 EM 算法-Jensen 不等式

EM 算法的收敛性
证明p(Xθ)p(\mathbf {X} \mid {\boldsymbol {\theta }})是单调递增的(概率大于等于 0,小于等于 1,单调增一定能收敛),就是证明logp(Xθ)\log p(\mathbf {X} \mid {\boldsymbol {\theta }})是单调递增的,即:
logp(Xθ(t+1))logp(Xθ(t))\log p(\mathbf {X} \mid {\boldsymbol {\theta^{(t+1)} }}) \geq \log p(\mathbf {X} \mid {\boldsymbol {\theta^{(t)} }})

参考深入理解 EM 算法-收敛性证明

广义 EM(generalized expectation maximization (GEM) algorithm)
我们以 EM 算法的 ELBO+KL 形式为例(书中的形式自己了解,这里不做介绍)
logP(Xθ)=Zq(Z)logP(X,Zθ)q(Z)dZ+DKL(q(Z)P(ZX,θ))\log P(X|\theta) = \int_Zq(Z)\log \frac{P(X,Z|\theta)}{q(Z)}dZ + D_{KL}(q(Z)||P(Z|X,\theta))
为了简单起见,我们令
L(q,θ)=Zq(Z)logP(X,Zθ)q(Z)dZ\mathcal{L}(q, \theta) = \int_Zq(Z)\log \frac{P(X,Z|\theta)}{q(Z)}dZ
其实就是变分函数(泛函的极值问题,输入是函数 q)
在前面,我们一直是让DKL=0D_{KL} = 0,然后最大化 L(q,θ)\mathcal{L}(q,\theta) 以增大 logP(Xθ)\log P(X|\theta) 。但是我们现在无法直接让q(Z)=P(ZX,θ)q(Z) = P(Z|X,\theta) ,因为后验本身也比较难求,所以我们希望能够求出某个 q(Z)q(Z) ,能够使得在固定 θ\thetaDKLD_{KL} 尽可能小,尽可能等于 0:
q=argminqDKL(q(Z)P(ZX,θ))=argmaxqL(q,θ)q = \arg\min\limits_q D_{KL}(q(Z)||P(Z|X,\theta))= \arg\max\limits_q \mathcal{L}(q, \theta)

参考深入理解 EM 算法-广义 EM

高斯混合模型Gaussian mixture model):

有样本(观测数据)Data={x1,...,xN}Data = \{x_1,...,x_N\},生成模型,对完全数据T={(x1,z1),...,(xN,zN)}T = \{(x_1,z_1),...,(x_N,z_N)\}建模,给定输入xix_i预测:arg maxk{1,,K}P(Zi=kxi;Θ^)\text{arg max}_{k \in \{1, \dots, K \}} P(Z_i = k \mid \boldsymbol{x}_i ; \hat{\Theta}),想要预测模型,就需要求出模型的参数Θ^\hat{\Theta},其中
P(Zi=kxi;Θ^)=p(xiZi=k;Θ^)P(Zi=k;Θ^)h=1Kp(xiZi=h;Θ^)P(Zi=h;Θ^)=ϕ(xiμ^k,Σ^k)α^kh=1Kϕ(xiμh^,Σ^h)α^h\begin{align*} P(Z_i = k \mid \boldsymbol{x}_i ; \hat{\Theta}) &= \frac{p(\boldsymbol{x}_i \mid Z_i = k ; \hat{\Theta})P(Z_i = k; \hat{\Theta})}{\sum_{h=1}^K p(\boldsymbol{x}_i \mid Z_i = h ; \hat{\Theta})P(Z_i = h; \hat{\Theta})} \\ &= \frac{\phi(\boldsymbol{x}_i \mid \hat{\boldsymbol{\mu}}_k, \hat{\boldsymbol{\Sigma}}_k) \hat{\alpha}_k}{\sum_{h=1}^K \phi(\boldsymbol{x}_i \mid \hat{\boldsymbol{\mu}_h}, \hat{\boldsymbol{\Sigma}}_h) \hat{\alpha}_h} \end{align*}

具体算法求解参考:Gaussian mixture models 以及机器学习-白板推导系列(十一)-高斯混合模型 GMM(Gaussian Mixture Model)
极大似然估计、EM 算法及高斯混合模型
EM 算法与 GMM(高斯混合聚类)Jensen 不等式和变分法两种推导
Expectation-maximization algorithm EM 算法
Mixture model
Gaussian Mixture Model

附加知识

变分推断

变分推断Variational Inference)也称为变分贝叶斯(Variational Bayesian),而变分法主要是研究变分问题,即泛函的极值问题(函数的输入也是函数),根据贝叶斯公式,后验概率
P(ZX)=P(XZ)P(Z)P(X)=P(XZ)P(Z)ZP(X,Z)dZ{\displaystyle P(\mathbf {Z} \mid \mathbf {X} )={\frac {P(\mathbf {X} \mid \mathbf {Z} )P(\mathbf {Z} )}{P(\mathbf {X} )}}={\frac {P(\mathbf {X} \mid \mathbf {Z} )P(\mathbf {Z} )}{\int _{\mathbf {Z} }P(\mathbf {X} ,\mathbf {Z} ')\,d\mathbf {Z} '}}}

上面公式中的积分对于很多情况下是不可行的(所以有些模型忽略了 P(x)),要么积分没有闭式解,要么是指数级别的计算复杂度,所以很难求出后验概率,这时我们需要寻找一个简单分布q(Z)P(ZX){\displaystyle q(\mathbf {Z} )\approx P(\mathbf {Z} \mid \mathbf {X} )},这样推断问题转化成一个泛函优化问题:
q(Z)^=arg minq(Z)候选的概率分布族QKL(q(Z)P(ZX))\hat{q(Z)} = \argmin_{q(Z) \in 候选的概率分布族Q} KL(q(Z)|P(Z|X))

我们有上面EM 算法的导出得到
logP(X)=DKL(qP)+L(q){\displaystyle \log P(\mathbf {X} )=D_{\mathrm {KL} }(q\parallel P)+{\mathcal {L}}(q)}
KL-divergence 大于等于 0,所以logP(X)L(q)\log P(\mathbf {X} ) \geq {\mathcal {L}}(q),所以L(q){\mathcal {L}}(q)称为证据下界(Evidence Lower BOund,ELBO),也就是所谓的变分函数。

KL-divergence 中有后验概率 P(Z|X),本身就是难以计算,才想找个简单分布 q(z)来近似,因此我们不能直接优化 KL-divergence,进而转化成优化 ELBO
q(Z)^=arg maxq(z)QELBO(q,x)\hat{q(Z)} = \argmax_{q(z) \in Q} ELBO(q,x)

分布族 Q 一般选择是平均场(Mean field)分布族,即可以将 Z 拆分为多组相互独立的变量,那么
q(Z)=m=1Mqm(Zm)q({\mathbf {Z}})=\prod _{{m=1}}^{M}q_{m}({\mathbf {Z}}_{m})

那么
ELBO(q,x)=Zq(Z)logP(X,Zθ)q(Z)dZ=Z(m=1Mqm(Zm))logP(X,Zθ)dZZ(m=1Mqm(Zm))logm=1Mqm(Zm)dZELBO(q,x) = \int_Zq(Z)\log \frac{P(X,Z|\theta)}{q(Z)}dZ \\ = \int_Z (\prod _{{m=1}}^{M}q_{m}({\mathbf {Z}}_{m}) ) \log P(X,Z|\theta)dZ - \int_Z (\prod _{{m=1}}^{M}q_{m}({\mathbf {Z}}_{m}) ) \log \prod _{{m=1}}^{M}q_{m}({\mathbf {Z}}_{m})dZ

假设我们只关心其中一个子集(分量)ZjZ_j的近似分布qj(Zj)q_j(Z_j)(先求一个子集,其它子集也就求出来了)
先看减号后面的项(进行展开):
Z(m=1Mqm(Zm))m=1Mlogqm(Zm)dZ=Zq1(Z1).q2(Z2)...qM(ZM).[logq1(Z1)+...+logqM(ZM)]dZ\int_Z (\prod _{{m=1}}^{M}q_{m}({\mathbf {Z}}_{m}) ) \sum_{m=1}^M \log q_{m}({\mathbf {Z}}_{m})dZ \\ = \int_Z q_1(Z_1).q_2(Z_2)...q_M(Z_M).[\log q_1(Z_1)+...+\log q_M(Z_M)]dZ

只看其中一项
Z1Z2...ZMq1(Z1).q2(Z2)...qM(ZM).logq1(Z1)dZ1dZ2...dZM=Z1q1(Z1).logq1(Z1)dZ1.Z2q2(Z2)dZ2....ZMqM(ZM)dZM\int_{Z_1Z_2...Z_M} q_1(Z_1).q_2(Z_2)...q_M(Z_M).\log q_1(Z_1){dZ_1dZ_2...dZ_M} = \int_{Z_1} q_1(Z_1).\log q_1(Z_1)dZ_1.\int_{Z_2}q_2(Z_2)dZ_2....\int_{Z_M}q_M(Z_M)dZ_M

那么最终我们得到减号后面的项
m=1MZmqm(Zm).logqm(Zm)dZm\sum_{m=1}^M \int_{Z_m} q_m(Z_m).\log q_m(Z_m)dZ_m

前面说了,我们只关注其中一个子集ZjZ_j,其它mjm \neq j的对其来说可以看作常数项,得到
Zjqj(Zj).logqj(Zj)dZj+C\int_{Z_j} q_j(Z_j).\log q_j(Z_j)dZ_j + C

再看减号前面的项
Z1Z2...ZM(m=1Mqm(Zm))logP(X,Zθ)dZ1dZ2...dZM=Zjqj(Zj)(ZiijMqi(Zi)logP(X,Zθ)dZi)dZj=Zjqj(Zj)(EijMqi(Zi)[logP(X,Zθ)])dZj=Zjqj(Zj)(logp~(X,Zj))dZj\int_{Z_1Z_2...Z_M} (\prod _{{m=1}}^{M}q_{m}({\mathbf {Z}}_{m}) ) \log P(X,Z|\theta){dZ_1dZ_2...dZ_M} \\= \int_{Z_j}q_j(Z_j) \bigg(\int_{Z_i} \prod_{i \neq j}^M q_i(Z_i) \log P(X,Z|\theta) dZ_i \bigg)dZ_j \\= \int_{Z_j}q_j(Z_j) \bigg(E_{\prod_{i \neq j}^M q_i(Z_i)}[\log P(X,Z|\theta)] \bigg)dZ_j \\= \int_{Z_j}q_j(Z_j) \bigg(\log \tilde{p}(X,Z_j) \bigg)dZ_j

最终:
ELBO(q,x)=Zjqj(Zj)(logp~(X,Zj))dZjZjqj(Zj).logqj(Zj)dZjC=KL(qj(Zj)p~(X,Zj))CELBO(q,x) =\int_{Z_j}q_j(Z_j) \bigg(\log \tilde{p}(X,Z_j) \bigg)dZ_j - \int_{Z_j} q_j(Z_j).\log q_j(Z_j)dZ_j - C \\= -KL(q_j(Z_j) \| \tilde{p}(X,Z_j)) -C
所以:
arg maxELBO(q,x)=arg minKL(qj(Zj)p~(X,Zj))\argmax ELBO(q,x) = \argmin KL(q_j(Z_j) \| \tilde{p}(X,Z_j))
qj(Zj)=p~(X,Zj)q_j(Z_j) = \tilde{p}(X,Z_j)取 KL 最小值

具体求最优算法这里不做介绍,可以参考变分推断(二)—— 进阶 以及 【一文学会】变分推断及其求解方法

参考11.4 节变分推断

参考文献

[9-1] Dempster AP,Laird NM,Rubin DB. Maximum-likelihood from incomplete data via theEM algorithm. J. Royal Statist. Soc. Ser. B.,1977,39

[9-2] Hastie T,Tibshirani R,Friedman J. The Elements of Statistical Learning: DataMining,Inference,and Prediction. Springer-Verlag,2001(中译本:统计学习基础——数据挖掘、推理与预测。范明,柴玉梅,昝红英等译。北京:电子工业出版社,2004)

[9-3] McLachlan G,Krishnan T. The EM Algorithm and Extensions. New York: John Wiley& Sons,1996

[9-4] 茆诗松,王静龙,濮晓龙。高等数理统计。北京:高等教育出版社;海登堡:斯普林格出版社,1998

[9-5] Wu CFJ. On the convergence properties of the EM algorithm. The Annals ofStatistics,1983,11: 95–103

[9-6] Radford N,Geoffrey H,Jordan MI. A view of the EM algorithm that justifiesincremental,sparse,and other variants. In: Learning in Graphical Models. Cambridge,MA: MITPress,1999,355–368

第 10 章 隐马尔可夫模型

隐马尔可夫模型Hidden Markov Model,HMM)是可用于标注问题的统计学习模型,描述由隐藏的马尔可夫链随机生成观测序列的过程,属于生成模型。

马尔可夫模型是关于时序的概率模型,描述由一个隐藏的马尔可夫链随机生成不可观测的状态随机序列,再由各个状态生成一个观测而产生观测随机序列的过程。
隐藏的马尔可夫链随机生成的状态的序列,称为状态序列(state sequence);每个状态生成一个观测,而由此产生的观测的随机序列,称为观测序列(observation sequence)。序列的每一个位置又可以看作是一个时刻。

隐马尔可夫模型由初始概率分布、状态转移概率分布以及观测概率分布确定。
隐马尔可夫模型的形式定义如下:
设 Q 是所有可能的状态的集合,N 是可能的状态数;V 是所有可能的观测的集合,M 是可能的观测数。
Q={q1,q2,...,qN},V={v1,v2,...,vM}Q = \{q_1,q_2,...,q_N\} , V= \{v_1,v_2,...,v_M\}
长度为 T 的状态序列I=(i1,i2,...,iT)I = (i_1,i_2,...,i_T)以及与状态序列对应的长度为 T 的观测序列O=(o1,o2,...,oT)O = (o_1,o_2,...,o_T)

状态转移矩阵(状态转移概率分布):(就是初始化参数transmat_prior,也可以用 params 和求出的属性 transmat_)
A=[aij]N×NA=[a_{ij}]_{N\times N}
其中aij=P(i_t+1=qjit=qi),下标i,j=1,...,Na_{ij} = P(i\_{t+1} = q_j | i_t = q_i) ,下标 i,j = 1,...,N表示在时刻tt处于状态qiq_i的条件下 在时刻t+1t+1转移到状态qjq_j的概率

观测矩阵(观测概率分布):(对于 MultinomialHMM 用 params 和求出的属性 emissionprob_,叫发生概率矩阵;对于 GMMHMM 有 n_mix 、means_prior、covars_prior ;对于 GaussianHMM 有 means_prior、covars_prior )
B=[bj(k)]N×MB = [b_j(k)]_{N \times M}
其中bj(k)=P(ot=vkit=qj),k=1,...,M,j=1,...,Nb_j(k) = P(o_t = v_k | i_t = q_j) ,k = 1,...,M,j = 1,...,N表示在时刻tt处于状态qjq_j的条件下生成观测vkv_k的概率

初始状态概率向量(初始概率分布):(就是初始化参数startprob_prior和求出的属性 startprob_ )
π=(πi)\pi = (\pi_i)
其中πi=P(i1=qi),下标i=1,...,N\pi_i = P(i_1 =q_i) ,下标i = 1,...,N表示时刻t=1t=1时 处于状态qiq_i的概率

隐马尔可夫模型由初始状态概率向量π\pi、状态转移概率矩阵AA和观测概率矩阵BB决定。
π\piAA决定状态序列,BB决定观测序列。
因此,一个隐马尔可夫模型可以用三元符号表示,即
λ=(A,B,π)\lambda = (A,B,\pi)
称为隐马尔可夫模型的三要素。

状态转移概率矩阵AA与初始状态概率向量π\pi确定了隐藏的马尔可夫链,生成不可观测的状态序列。观测概率矩阵BB确定了如何从状态生成观测,与状态序列综合确定了如何产生观测序列。

从定义可知,隐马尔可夫模型作了两个基本假设

  1. 齐次马尔可夫性假设,即假设隐藏的马尔可夫链在任意时刻 t 的状态只依赖于其前一时刻的状态,与其他时刻的状态及观测无关,也与时刻 t 无关。
    P(itit1,ot1,...,i1,o1)=P(itit1),t=1,2,...,TP(i_{t}|i_{t-1},o_{t-1},...,i_{1},o_{1}) = P(i_{t}|i_{t-1}), t=1,2,...,T
  2. 观测独立性假设,即假设任意时刻的观测只依赖于该时刻的马尔可夫链的状态,与其他观测及状态无关。
    P(otiT,oT,iT1,oT1,...,it+1,ot+1,it,ot,it1,ot1,...,i1,o1)=P(otit)P(o_{t}|i_{T},o_{T},i_{T-1},o_{T-1},...,i_{t+1},o_{t+1},i_{t},o_{t},i_{t-1},o_{t-1},...,i_{1},o_{1}) = P(o_{t}|i_{t})

隐马尔可夫模型的三个基本问题

  1. 概率计算问题。给定模型λ=(A,B,π)\lambda = (A,B,\pi)和观测序列O=(o1,o2,...,oT)O = (o_1,o_2,...,o_T),计算在模型λ\lambda下观测序列OO出现的概率P(Oλ)P(O|\lambda)

  2. 学习问题。已知观测序列O=(o1,o2,...,oT)O = (o_1,o_2,...,o_T),估计模型λ=(A,B,π)\lambda = (A,B,\pi)参数,使得在该模型下观测序列概率P(Oλ)P(O|\lambda)最大。即用极大似然估计的方法估计参数。(λMLE=arg maxλP(Oλ)\lambda_{MLE}=\argmax_{\lambda}P(O|\lambda),使用 EM 算法求解。)

  3. 预测问题,也称为解码(decoding)问题。已知模型λ=(A,B,π)\lambda = (A,B,\pi)和观测序列O=(o1,o2,...,oT)O = (o_1,o_2,...,o_T),求对给定观测序列条件概率P(IO)P(I|O)最大的状态序列I=(i1,i2,...,iT)I = (i_1,i_2,...,i_T)。即给定观测序列,求最有可能的对应的状态序列。(Viterbi 算法求I^=arg maxIP(IO,λ)\hat{I}=\argmax_{I}P(I|O,\lambda)

概率计算问题
引入隐变量,对完全数据建模(这里还是一样P(Oλ),P(O;λ)P(O|\lambda),P(O;\lambda)是一样的,λ\lambda是参数)
P(Oλ)=IP(O,Iλ)=IP(OI,λ)P(Iλ)P(O|\lambda) = \sum_{I}P(O,I|\lambda)= \sum_{I}P(O|I,\lambda)P(I|\lambda)
根据乘法规则(概率论基础教程 51 页,注意P(i1λ)=P(i1)P(i_1|\lambda) = P(i_1))以及马尔可夫假设有:
P(Iλ)=P(i1,i2,...,iTλ)=P(i1).P(i2i1,λ).P(i3i1,i2,λ)...P(iTi1,i2,...,iT1,λ)=P(i1)t=2TP(iti1,i2,...,it1,λ)=P(i1)t=2TP(itit1,λ)=πi1t=2Tait1itP(I|\lambda) = P(i_1,i_2,...,i_T|\lambda)=P(i_1).P(i_2|i_1,\lambda).P(i_3|i_1,i_2,\lambda)...P(i_T|i_1,i_2,...,i_{T-1},\lambda) \\= P(i_1)\prod_{t=2}^T P(i_t|i_1,i_2,...,i_{t-1},\lambda) \\= P(i_1)\prod_{t=2}^T P(i_t|i_{t-1},\lambda) \\= \pi_{i_1}\prod_{t=2}^T a_{i_{t-1}i_{t}}
根据乘法规则以及观测独立性假设有:
P(OI,λ)=P(o1,o2,...,oTi1,i2,...,iT,λ)=P(o1i1,i2,...,iT,λ).P(o2o1,i1,i2,...,iT,λ).P(o3o1,o2,i1,i2,...,iT,λ)...P(oTo1,o2,...,oT1,i1,i2,...,iT,λ)=P(o1i1,λ).P(o2i2,λ)...P(oTiT,λ)=t=1Tbit(ot)P(O|I,\lambda) = P(o_1,o_2,...,o_T|i_1,i_2,...,i_{T},\lambda) \\= P(o_1|i_1,i_2,...,i_{T},\lambda).P(o_2|o_1,i_1,i_2,...,i_{T},\lambda).P(o_3|o_1,o_2,i_1,i_2,...,i_{T},\lambda)...P(o_T|o_1,o_2,...,o_{T-1},i_1,i_2,...,i_{T},\lambda) \\ = P(o_1|i_1,\lambda).P(o_2|i_2,\lambda)...P(o_T|i_T,\lambda) \\= \prod_{t=1}^Tb_{i_t}(o_t)
那么

P(O,Iλ)=P(OI,λ)P(Iλ)=πi1t=2Tait1itt=1Tbit(ot)=πi1bi1(o1).ai1i2bi2(o2)...aiT1iTbiT(oT)=πi1bi1(o1)t=2Tait1itbit(ot)P(O,I|\lambda) = P(O|I,\lambda)P(I|\lambda) = \pi_{i_1}\prod_{t=2}^T a_{i_{t-1}i_{t}}\prod_{t=1}^Tb_{i_t}(o_t) \\= \pi_{i_1}b_{i_1}(o_1) .a_{i_1i_2}b_{i_2}(o_2)...a_{i_{T-1}i_T}b_{i_T}(o_T) = \pi_{i_1}b_{i_1}(o_1)\prod_{t=2}^T a_{i_{t-1}i_{t}}b_{i_t}(o_t)

概率计算问题- 直接由上面计算概率可得
P(Oλ)=IP(O,Iλ)=IP(OI,λ)P(Iλ)=i1,i2,...,iTπi1bi1(o1)t=2Tait1itbit(ot)=i1N...iTNπi1bi1(o1)t=2Tait1itbit(ot)P(O|\lambda) = \sum_{I}P(O,I|\lambda)= \sum_{I}P(O|I,\lambda)P(I|\lambda) \\= \sum_{i_1,i_2,...,i_T} \pi_{i_1}b_{i_1}(o_1)\prod_{t=2}^T a_{i_{t-1}i_{t}}b_{i_t}(o_t) \\= \sum_{i_1 \in N}...\sum_{i_T\in N} \pi_{i_1}b_{i_1}(o_1)\prod_{t=2}^T a_{i_{t-1}i_{t}}b_{i_t}(o_t)
时间复杂度O(TNT)O(TN^{T}),所以不可行。

上面说过直接求不好求,有以下方法可求得:
概率计算问题- 前向计算
首先我们定义前向概率αt(i)=P(o1,o2,...,ot,it=qiλ)\alpha_t(i) = P(o_1,o_2,...,o_t,i_t=q_i | \lambda),表示时刻tt部分观测序列为o1,o2,...,oto_1,o_2,...,o_t且状态为qiq_i的概率,那么
P(Oλ)=i=1NP(O,iT=qiλ)=i=1NP(o1,...,oT,iT=qiλ)=i=1NαT(i)P(O|\lambda) = \sum_{i=1}^N P(O,i_T=q_i|\lambda) = \sum_{i=1}^N P(o_1,...,o_T,i_T=q_i|\lambda) = \sum_{i=1}^N \alpha_T(i)

其实P(Oλ)=j=1NP(O,i1=qjλ)=...=j=1NP(O,it=qjλ)=i=1Nj=1NP(O,i1=qi,i2=qjλ)P(O|\lambda) = \sum_{j=1}^N P(O,i_1=q_j|\lambda) =...= \sum_{j=1}^N P(O,i_t=q_j|\lambda) = \sum_{i=1}^N\sum_{j=1}^N P(O,i_1=q_i,i_2=q_j|\lambda),注意这里是小 t,只不过我们定义了前向概率,并且O=(o1,...,oT)O=(o_1,...,o_T)

所以我们只要求出αT(i)\alpha_T(i),如何求?依次α1(i)...αt+1(i)...αT(i)\alpha_1(i) ... \alpha_{t+1}(i) ... \alpha_T(i)

α1(i)=P(o1,i1=qiλ)=P(i1=qiλ)P(o1i1=qi,λ)=πibi(o1)αt+1(i)=P(o1,o2,...,ot,ot+1,it+1=qiλ)=j=1NP(o1,o2,...,ot,ot+1,it+1=qi,it=qjλ)=j=1NP(ot+1o1,..,ot,it+1=qi,it=qj,λ)P(o1,o2,...,ot,it+1=qi,it=qjλ)=j=1NP(ot+1it+1=qi)P(o1,o2,...,ot,it+1=qi,it=qjλ)=j=1NP(ot+1it+1=qi)P(it+1=qio1,o2,...,ot,it=qj,λ)P(o1,o2,...,ot,it=qjλ)=j=1NP(ot+1it+1=qi)P(it+1=qiit=qj,λ)P(o1,o2,...,ot,it=qjλ)=j=1NP(ot+1it+1=qi)P(it+1=qiit=qj,λ)αt(j)=P(ot+1it+1=qi)j=1NP(it+1=qiit=qj,λ)αt(j)=[j=1Nαt(j)aji]bi(ot+1)\alpha_1(i) = P(o_1,i_1=q_i | \lambda) =P(i_1=q_i | \lambda)P(o_1|i_1=q_i , \lambda) = \pi_ib_i(o_1) \\ \vdots\\ \alpha_{t+1}(i) = P(o_1,o_2,...,o_t,o_{t+1},i_{t+1}=q_i | \lambda) \\=\sum_{j=1}^N P(o_1,o_2,...,o_t,o_{t+1},i_{t+1}=q_i,i_{t}=q_j | \lambda) \\ =\sum_{j=1}^NP(o_{t+1}|o_1,..,o_t,i_{t+1}=q_i,i_{t}=q_j,\lambda)P(o_1,o_2,...,o_t,i_{t+1}=q_i,i_{t}=q_j | \lambda) \\=\sum_{j=1}^NP(o_{t+1}|i_{t+1}=q_i)P(o_1,o_2,...,o_t,i_{t+1}=q_i,i_{t}=q_j | \lambda) \\= \sum_{j=1}^NP(o_{t+1}|i_{t+1}=q_i)P(i_{t+1}=q_i | o_1,o_2,...,o_t,i_{t}=q_j,\lambda)P(o_1,o_2,...,o_t,i_{t}=q_j | \lambda) \\=\sum_{j=1}^NP(o_{t+1}|i_{t+1}=q_i)P(i_{t+1}=q_i | i_{t}=q_j,\lambda)P(o_1,o_2,...,o_t,i_{t}=q_j | \lambda) \\=\sum_{j=1}^NP(o_{t+1}|i_{t+1}=q_i)P(i_{t+1}=q_i | i_{t}=q_j,\lambda)\alpha_t(j) \\= P(o_{t+1}|i_{t+1}=q_i)\sum_{j=1}^NP(i_{t+1}=q_i | i_{t}=q_j,\lambda)\alpha_t(j) \\= \bigg[\sum_{j=1}^N\alpha_t(j)a_{ji} \bigg] b_i(o_{t+1})

概率计算问题- 后向计算
首先我们定义后向概率βt(i)=P(ot+1,ot+2,...,oTit=qi,λ)\beta_t(i) = P(o_{t+1},o_{t+2},...,o_T|i_t=q_i , \lambda),表示时刻状态为qiq_i的条件下,从t+1t+1TT的部分观测序列为ot+1,ot+2,...,oTo_{t+1},o_{t+2},...,o_T概率,那么
P(Oλ)=i=1NP(O,i1=qiλ)=i=1NP(o1,...,oT,i1=qiλ)=i=1NP(o1,...,oTi1=qi,λ)P(i1=qiλ)=i=1NP(o1o2,...,oT,i1=qi,λ)P(o2,...,oTi1=qi,λ)P(i1=qiλ)=i=1NP(o1i1=qi,λ)P(o2,...,oTi1=qi,λ)P(i1=qiλ)=i=1Nbi(o1)β1(i)πiP(O|\lambda) = \sum_{i=1}^N P(O,i_1=q_i|\lambda) = \sum_{i=1}^N P(o_1,...,o_T,i_1=q_i|\lambda) \\= \sum_{i=1}^N P(o_1,...,o_T|i_1=q_i,\lambda)P(i_1=q_i|\lambda) \\= \sum_{i=1}^N P(o_1|o_2,...,o_T,i_1=q_i,\lambda)P(o_2,...,o_T|i_1=q_i,\lambda)P(i_1=q_i|\lambda) \\ = \sum_{i=1}^N P(o_1|i_1=q_i,\lambda)P(o_2,...,o_T|i_1=q_i,\lambda)P(i_1=q_i|\lambda) \\= \sum_{i=1}^N b_i(o_1)\beta_1(i)\pi_i
所以我们只要求出β1(i)\beta_1(i),如何求?依次βT(i)...β1t1(i)...β1(i)\beta_T(i) ... \beta_1{t-1}(i) ... \beta_1(i)

βT(i)=P(iT=qi,λ)=1βt(i)=P(ot+1,ot+2,...,oTit=qi,λ)=j=1NP(ot+1,ot+2,...,oT,it+1=qjit=qi,λ)=j=1NP(ot+1,ot+2,...,oTit+1=qj,it=qi,λ)P(it+1=qjit=qi,λ)条件前面没有ot(根据概率图也能得出给定it+1时,itot+1,...,oT无关)=j=1NP(ot+1,ot+2,...,oTit+1=qj,λ)P(it+1=qjit=qi,λ)=j=1NP(ot+1ot+2,...,oT,it+1=qj,λ)P(ot+2,...,oTit+1=qj,λ)P(it+1=qjit=qi,λ)=j=1NP(ot+1it+1=qj,λ)P(ot+2,...,oTit+1=qj,λ)P(it+1=qjit=qi,λ)=j=1Nbj(ot+1)βt+1(j)aij\beta_T(i) = P(i_T = q_i,\lambda) = 1 \\ \vdots \\ \beta_t(i) = P(o_{t+1},o_{t+2},...,o_T|i_t=q_i , \lambda) \\= \sum_{j=1}^N P(o_{t+1},o_{t+2},...,o_T,i_{t+1}=q_j|i_t=q_i , \lambda) \\= \sum_{j=1}^N P(o_{t+1},o_{t+2},...,o_T|i_{t+1}=q_j,i_t=q_i , \lambda) P(i_{t+1}=q_j|i_t=q_i , \lambda) \\ 条件前面没有o_t(根据概率图也能得出给定i_{t+1}时,i_t与o_{t+1},...,o_T无关) \\= \sum_{j=1}^N P(o_{t+1},o_{t+2},...,o_T|i_{t+1}=q_j, \lambda) P(i_{t+1}=q_j|i_t=q_i , \lambda) \\= \sum_{j=1}^N P(o_{t+1}|o_{t+2},...,o_T,i_{t+1}=q_j, \lambda)P(o_{t+2},...,o_T|i_{t+1}=q_j, \lambda) P(i_{t+1}=q_j|i_t=q_i , \lambda)\\= \sum_{j=1}^N P(o_{t+1}|i_{t+1}=q_j, \lambda)P(o_{t+2},...,o_T|i_{t+1}=q_j, \lambda) P(i_{t+1}=q_j|i_t=q_i , \lambda) \\ =\sum_{j=1}^N b_j(o_{t+1}) \beta_{t+1}(j) a_{ij}

预测问题,也称为解码(decoding)问题
维特比算法(Viterbi algorithm)实际是用动态规划解隐马尔可夫模型预测问题,即用动态规划(dynamic programming)求概率最大路径(最优路径),这里的最优路径就是最优状态序列II

请参考书籍和机器学习-白板推导系列(十四)-隐马尔可夫模型 HMM(Hidden Markov Model)

这一类模型需要求解的问题的大体框架为:
其中XX代表观测序列,ZZ代表隐变量序列,λ\lambda代表参数。

{RepresentationProbabilistic graphical modelLearningλMLE=argmaxλP(Xλ)Baum Welch Algorithm(EM)Inference{DecodingZ=argmaxZP(ZX,λ)orP(z1,z2,,ztx1,x2,,xt,λ)Viterbi AlgorithmProb of evidenceP(Xλ)Forward Algorithm,Backward AlgorithmFilteringP(ztx1,x2,,xt,λ)(online)Forward AlgorithmSmothingP(ztx1,x2,,xT,λ)(offline)Forward-Backward AlgorithmPrediction{P(zt+1,zt+2,...x1,x2,,xt,λ)P(xt+1,xt+2,...x1,x2,,xt,λ)}Forward Algorithm\begin{cases} Representation & \text{Probabilistic graphical model} \\ Learning & \lambda_{MLE}=arg \underset{\lambda}{\max} P(X|\lambda) \boxed{\text{Baum Welch Algorithm(EM)}}\\ Inference & \begin{cases} Decoding & Z=arg\underset{Z}{\max}P(Z|X,\lambda) or P(z_1,z_2,\cdots,z_t|x_1,x_2,\cdots,x_t,\lambda) \boxed{\text{Viterbi Algorithm}}\\ \text{Prob of evidence} & P(X|\lambda) \boxed{\text{Forward Algorithm,Backward Algorithm}} \\ Filtering & P(z_t|x_1,x_2,\cdots,x_t,\lambda) \boxed{\text{(online)Forward Algorithm}}\\ Smothing & P(z_t|x_1,x_2,\cdots,x_T,\lambda) \boxed{\text{(offline)Forward-Backward Algorithm}}\\Prediction & \begin{Bmatrix} P(z_{t+1},z_{t+2},...|x_1,x_2,\cdots,x_t,\lambda) \\ P(x_{t+1},x_{t+2},...|x_1,x_2,\cdots,x_t,\lambda) \end{Bmatrix} \boxed{\text{Forward Algorithm}} \end{cases}\\ \end{cases}

Filtering problem (stochastic processes)
Smoothing problem (stochastic processes)

附加知识

随机过程

随机过程Stochastic process):

(Ω,F,P)(\Omega ,{\mathcal {F}},P)为一个概率空间Probability space),Ω\Omega样本空间sample space), F\mathcal {F} 是(Sigma-algebra),PP 是(Probability measure);
(S,Σ)(S,\Sigma )为可测量的空间(measurable space),SS为随机变量的集合
{X(t):tT}or{X(t,ω):tT}{\displaystyle \{X(t):t\in T\}} or {\displaystyle \{X(t,\omega ):t\in T\}}

其中X(t)X(t)是一个随机变量,(在自然科学的许多问题中tt表示时间,那么X(t)X(t)表示在时刻tt观察到的值);ωΩ\omega \in \OmegaTT是指标集 or 参数集(index set or parameter set),一般表示时间或空间,如:离散T={0,1,2,...}T=\{0,1,2,...\}一般称为随机序列或时间序列,连续T=[a,b],a可以取0或者,b可以取+T=[a,b] ,a可以取0或者 -\infty,b可以取+\infty

映射X(t,ω):T×ΩRX(t,\omega):T \times \Omega \to R,即X(.,.)X(.,.)是定义在T×ΩT \times \Omega上的二元值函数;
tT\forall t \in T(固定tTt \in T),X(t,.)X(t,.)是定义在样本空间Ω\Omega上的函数,称为随机变量;
ωΩ\forall \omega \in \Omega,映射X(.,ω):TSX(.,\omega):T \to S(其实就是固定ωΩ\omega \in \Omega,变成关于T的函数),被称为样本函数(sample function),特别是当TT表示时间时,称为随机过程{X(t,ω):tT}{\displaystyle \{X(t,\omega ):t\in T\}}样本路径(sample path)。

傻傻分不清楚的马尔可夫

参见书中第19章

马尔可夫性质Markov property):
如果随机过程Stochastic process)的未来状态的条件概率分布(以过去和现在值为条件)仅取决于当前状态,则随机过程具有马尔可夫性质;与此属性的过程被认为是马氏马尔可夫过程Markov process)。最著名的马尔可夫过程是马尔可夫链Markov chain)。布朗运动Brownian motion)是另一个著名的马尔可夫过程。马尔可夫过程是不具备记忆特质的(Memorylessness

一阶 离散
P(Xn=xnXn1=xn1,,X0=x0)=P(Xn=xnXn1=xn1).{\displaystyle P(X_{n}=x_{n}\mid X_{n-1}=x_{n-1},\dots ,X_{0}=x_{0})=P(X_{n}=x_{n}\mid X_{n-1}=x_{n-1}).}

m 阶 离散
Pr(Xn=xnXn1=xn1,Xn2=xn2,,X1=x1)=Pr(Xn=xnXn1=xn1,Xn2=xn2,,Xnm=xnm) for n>m{\displaystyle {\begin{aligned}&\Pr(X_{n}=x_{n}\mid X_{n-1}=x_{n-1},X_{n-2}=x_{n-2},\dots ,X_{1}=x_{1})\\={}&\Pr(X_{n}=x_{n}\mid X_{n-1}=x_{n-1},X_{n-2}=x_{n-2},\dots ,X_{n-m}=x_{n-m}){\text{ for }}n>m\end{aligned}}}

时间齐次(Time-homogeneous)
Pr(Xt+s=xXt+s1=y)=Pr(Xt=xXt1=y){\displaystyle \Pr(X_{t+s}=x\mid X_{t+s-1}=y)=\Pr(X_{t}=x\mid X_{t-1}=y)}

马尔可夫模型Markov model):

. System state is fully observable System state is partially observable
System is autonomous Markov chain Hidden Markov model
System is controlled Markov decision process Partially observable Markov decision process

马尔可夫模型是具有马尔可夫性假设的随机过程,最简单的马尔可夫模型是马尔可夫链Markov chain

. Countable state space(对应随机过程的离散Ω\Omega Continuous or general state space(对应随机过程的连续Ω\Omega
Discrete-time(对应随机过程的离散TT (discrete-time) Markov chain on a countable or finite state space Markov chain on a measurable state space (for example, Harris chain)
Continuous-time(对应随机过程的连续TT Continuous-time Markov process or Markov jump process Any continuous stochastic process with the Markov property (for example, the Wiener process)

{Xnn=1,2,}\{X_n|n=1,2,\cdots\}是有限个或可数个可能值的随机过程。除非特别提醒,这个随机过程的可能值的集合都将记为非负整数集{0,1,2,}\{0,1,2,\cdots\}
如果 Xn=iX_n=i,那么称该过程在时刻 tt 在状态 ii ,并假设Pi,jP_{i,j} 称为(单步(one-step))转移概率(transition probability),表示处在ii状态的随机变量下一时刻处在jj 状态的概率,如果对所有的状态 i0,i1,,in1,i,ji_0,i_1,\cdots,i_{n-1},i,j及任意n0n\ge 0P{Xn+1=jXn=i,Xn1=in1,,X1=i1,X0=i0}=Pi,jP\{X_{n+1}=j|X_n=i,X_{n-1}=i_{n-1},\cdots,X_1=i_1,X_0=i_0\}=P_{i,j},这样的过程称为马尔可夫链(Markov chain)。

一个随机过程{X(t):t0}\{X(t):t \geq 0\}如果tR+t \in \mathbb{R}_+则称为连续时间的马尔科夫链,如果tN0t \in \mathbb{N}_0则称为离散时间的马尔科夫链

根据 Pi,jP_{i,j}的定义显然有Pi,j0,  i,j0;    j=0Pi,j=1,  i=0,1,P_{i,j}\ge0,\;i,j\ge0;\;\;\sum_{j=0}^\infty P_{i,j}=1,\;i=0,1,\cdots,
Pi,jP_{i,j} 记录从 iijj 的(单步)转移(概率)矩阵transition matrix)也称为随机矩阵、概率矩阵、转移矩阵、替代矩阵或马尔可夫矩阵
Pi,j=(Pin,in+1)=[P0,0P0,1P0,2P1,0P1,1P1,2Pi,0Pi,1Pi,2]\mathbf{P}_{i,j}=(P_{i_{n},i_{n+1}}) =\begin{bmatrix}P_{0,0}&P_{0,1}&P_{0,2}&\cdots\\P_{1,0}&P_{1,1}&P_{1,2}&\cdots\\\vdots&\vdots&\vdots\\P_{i,0}&P_{i,1}&P_{i,2}&\cdots\\\vdots&\vdots&\vdots\end{bmatrix}
现在定义 n 步(n-step)转移概率Pi,jnP_{i,j}^nPi,jn=P{Xn+k=jXk=i},  n0,i,j0P_{i,j}^n=P\{X_{n+k=j}|X_k=i\},\;n\ge 0,i,j\ge 0

右随机矩阵是一个非负实数的方阵,每个行总和为 1。
左随机矩阵是一个非负实数的方阵,每个列的总和为 1。
双随机矩阵是一个非负实数的方阵,每个行和每个列的总和为 1。

假设AA是马尔可夫矩阵,其性质有:

  1. 矩阵AA的 k 次幂AkA^k也是马尔可夫矩阵。
  2. 至少有一个特征值为 1,其特征值在[-1,1]区间,特征值为 1 对应的特征向量π\pi称为平稳概率向量(stationary probability vector)。
  3. 对于任意概率向量Probability vector)或者随机向量π0\pi_0limkAkπ0=π\lim_{k \to \infty} A^k \pi_0 = \pi (这里是在没有-1 特征值的情况下)。
  4. 对于任意概率向量μ0\mu_0μ1=Aμ0\mu_1 = A \mu_0也是概率向量。

特征值的求解:det(AλI)=0\det(A-\lambda I)=0

由于 AA 的每一列相加等于 1,所以 AIA−I 的每一列相加等于 0,这也就是说 AIA−I 的行是相关的,其行列式det(AI)=0\det(A-I)=0为零,所以 AIA−I奇异矩阵,所以 1 是 AA 的一个特征值。

对角化 A=PΛP1A = P \Lambda P^{-1} (参见线性代数及其应用279页,特征值相同特征向量不一定相同),其中Λ\Lambda是由AA的特征值组成的对角矩阵
μk=Akμ0=(PΛP1)kμ0=PΛkP1μ0\mu_k = A^k \mu_0 = (P \Lambda P^{-1})^k \mu_0 = P \Lambda^k P^{-1} \mu_0
不妨设 AA的特征向量和相应的特征值分别为 x1,...,xn{x_1},...,{x_n}λ1,...,λn\lambda_1,...,\lambda_n,可以用特征向量来做一组基,可以把空间中任何向量写成它的线性组合:μ0=c1x1++cnxn\mu_0 = c_1{x_1} + \cdots + c_n{x_n}
那么:
Akμ0=Akc1x1++Akcnxn=c1Akx1++cnAkxn=c1Ak1Ax1++cnAk1Axn=c1Ak1λ1x1++cnAk1λnxn=c1λ1kx1++cnλnkxn=i=1nciλikxiA^k\mathbf{\mu_0} = A^kc_1\mathbf{x_1} + \cdots + A^kc_n\mathbf{x_n}\\ = c_1A^k\mathbf{x_1} + \cdots + c_nA^k\mathbf{x_n} \\= c_1A^{k-1}A\mathbf{x_1} + \cdots + c_nA^{k-1}A\mathbf{x_n} \\= c_1A^{k-1}\lambda_1\mathbf{x_1} + \cdots + c_nA^{k-1}\lambda_n\mathbf{x_n}\\=c_1\lambda_1^k\mathbf{x_1} + \cdots + c_n\lambda_n^k\mathbf{x_n}\\=\sum_{i=1}^n{c_i\lambda_i^k\bm{x_i}}

不妨令λ1=1\lambda_1=1, 有λi1|\lambda_i|\leq 1,那么:
u=limkAku0=limki=1kciλikxi=c1x1\bm{u_\infty}=\lim_{k\to\infty}{A^k\bm{u_0}}=\lim_{k\to\infty}{\sum_{i=1}^k{c_i\lambda_i^k\bm{x_i}}}=c_1\bm{x_1}

因为uu_\infty是概率向量,而特征值为 1 对应的特征向量也是概率向量,所以c1=1c_1=1,得到u=x1\bm{u_\infty}=\bm{x_1}

值得注意的是,除包含 λ=1\lambda=1 的情形还包含 λ=1\lambda=-1 的情形。
上式如果除λ1=1\lambda_1=1还有λ2=1\lambda_2=-1,那么就有:
u=limki=1kciλikxi=c1x1+(1)kc2x2\bm{u_\infty}=\lim_{k\to\infty}{\sum_{i=1}^k{c_i\lambda_i^k\bm{x_i}}}=c_1\bm{x_1}+(-1)^k c_2\bm{x_2}

u=x1+(1)kx2\bm{u_\infty}=\bm{x_1}+(-1)^k\bm{x_2},此时kk为奇数和偶数结果是不同的,造成的结果就是在两种结果之间反复横跳,永远达不到稳态。

如:
A=[0110]A=\begin{bmatrix}0&1\\1&0\\\end{bmatrix}
其特征值为λ1=1λ2=1\lambda_1=1,\lambda_2=-1

也可以参考第 21 章 PageRank 算法

规划论

规划论又称数学规划,运筹学(Operations research)的一个分支。 规划论是指在既定条件(约束条件)下,按照某一衡量指标(目标函数)在多种 方案中寻求最优方案(取最大或最小值)。规划论包括线性规划、非线性规划和动态规划等,是一种优化算法或方法(Optimization algorithms and methods

数学优化(Mathematical optimization

优化技术

线性规划

当目标函数与约束条件都是线形的,则称为线性规划(Linear programming‎)。

求解方法:图解法(graphical method)、单纯形法(simplex algorithm)、对偶单纯形法等

非线性规划

除去线性规划,则为非线性规划(Nonlinear programming)。其中,凸规划(前面的章节有讲到凸优化)、二次规划(Quadratic programming)、几何规划都是一种特殊的非线性规划。

求解方法:拉格朗日乘子法、可行方向法、制约函数法(constrained function method )等。

内点法(Interior point methods)是一种求解线性规划或非线性凸优化问题的算法。

无约束优化问题

去除带约束的规划问题,则为无约束优化问题(Unconstrained convex optimization,对应的有约束优化(Constrained optimization))。

求解方法: 1、 最速下降法(也叫梯度下降) 2、 共轭梯度下降 3、 牛顿法 4、 拟牛顿法

动态规划

若规划问题与时间有关,则称为动态规划(Dynamic programming‎);

把多阶段过程转化为一系列单阶段问题,逐个求解,解决这类问题的方法称为动态规划,它是一种方法、考察问题的一种途径,但不是一种特殊的算法。 没有统一的标准模型,也没有构造模型的通用方法,甚至还没有判断一个问题能否构造动态规划模型的准则。这样就只能对每类问题进行具体分析,构造具体的模型。对于较复杂的问题在选择状态、决策、确定状态转移规律等方面需要丰富的想象力和灵活的技巧性,这就带来了应用上的局限性。

动态规划一般可分为线性动规,区域动规,树形动规,背包动规(Knapsack problem)四类。
线性动规:拦截导弹,合唱队形,挖地雷,建学校,剑客决斗等;
区域动规:石子合并, 加分二叉树,统计单词个数,炮兵布阵等;
树形动规:贪吃的九头龙,二分查找树,聚会的欢乐,数字三角形等;
背包问题:背包问题,完全背包问题,分组背包问题,二维背包,装箱问题,挤牛奶

随机规划

若规划问题与随机变量有关,则称为随机规划(Stochastic programming)。

随机动态规划

Stochastic dynamic programming

组合规划

若规划问题与有限个事物的排列组合有关,则称为组合规划(combinatorial optimization)

参考文献

[10-1] Rabiner L,Juang B. An introduction to hidden markov Models. IEEE ASSPMagazine,January 1986

[10-2] Rabiner L. A tutorial on hidden Markov models and selected applications in speechrecognition. Proceedings of IEEE,1989

[10-3] Baum L,et al. A maximization technique occuring in the statistical analysis of probabilistic functions of Markov chains. Annals of Mathematical Statistics,1970,41: 164–171

[10-4] Bilmes JA. A gentle tutorial of the EM algorithm and its application to parameter estimation for Gaussian mixture and hidden Markov models.

[10-5] Lari K,Young SJ. Applications of stochastic context-free grammars using the Inside-Outside algorithm,Computer Speech & Language,1991,5(3): 237–257

[10-6] Ghahramani Z. Learning Dynamic Bayesian Networks. Lecture Notes in ComputerScience,Vol. 1387,1997,168–197

以下来自隐马尔可夫模型

[10-7] J. Li, A. Najmi, R. M. Gray, Image classification by a two dimensional hidden Markov model,IEEE Transactions on Signal Processing , 48(2):517-33, February 2000. [2-D HMM] (download)

[10-8] J. Li, R. M. Gray, R. A. Olshen, Multiresolution image classification by hierarchical modeling with two dimensional hidden Markov models, IEEE Transactions on Information Theory , 46(5):1826-41, August 2000. [2-D MHMM] (download)

[10-9] J. Li, W. Miller, Significance of inter-species matches when evolutionary rate varies, Journal of Computational Biology , 10(3-4):537-554, 2003. [HMMMO] (download)

[10-10] J. Li, J. Z. Wang, Studying digital imagery of ancient paintings by mixtures of stochastic models, IEEE Transactions on Image Processing, 12(3):340-353, 2004. [Mixture of 2-D MHMMs] (download)

第 11 章 条件随机场

条件随机场Conditional random field, CRF)条件随机场(CRFs)是一类常用的统计建模方法(statistical modeling methods),常用于模式识别(pattern recognition)和机器学习,并用于结构预测(structured prediction)。

相关的机器学习库有PyStructpython-crfsuite

这里推荐学习:机器学习-白板推导系列(十七)-条件随机场CRF(Conditional Random Field) 以及论文Conditional Random Fields: Probabilistic Models for Segmenting and Labeling Sequence Data

条件随机场是在无向图上的判别模型。

条件随机场是给定一组输入随机变量条件下另一组输出随机变量的条件概率分布模型,其特点是假设输出随机变量构成马尔可夫随机场。
条件随机场可以用于不同的预测问题,本书仅论及它在标注问题的应用。因此主要讲述线性链(linear chain)条件随机场,这时,问题变成了由输入序列对输出序列预测的判别模型,形式为对数线性模型,其学习方法通常是极大似然估计或正则化的极大似然估计。

条件随机场(conditional random field)是给定随机变量X条件下,随机变量Y的马尔可夫随机场。
XXYY是随机变量,P(YX)P(Y|X)是在给定X的条件下YY的条件概率分布。若随机变量YY构成一个由无向图G(V,E)G=(V,E)表示的马尔可夫随机场,即
p(YvX,Yw,wv)=p(YvX,Yw,wv)p(\boldsymbol{Y}_v |\boldsymbol{X}, \boldsymbol{Y}_w, w \neq v) = p(\boldsymbol{Y}_v |\boldsymbol{X}, \boldsymbol{Y}_w, w \sim v)
对任意结点vv成立,则称条件概率分布P(YX)P(Y|X)为条件随机场。式中wvw \sim v表示在图G(V,E)G=(V,E)中与结点vv有边连接的所有结点wwwvw \neq v表示结点vv以外的所有结点,YvY_vYuY_uYwY_w为结点vvuuww对应的随机变量。

线性链条件随机场(linear chain conditional random field)假设X和Y有相同的图结构。

条件随机场在定义中并没有要求X和Y具有相同的结构。现实中,一般假设X和Y有相同的图结构。

X(X1,X2,...,Xn)Y(Y1Y2,...,Yn)X=(X_1,X_2,...,X_n),Y=(Y_1,Y_2,...,Y_n)均为线性链表示的随机变量序列,若在给定随机变量序列XX的条件下,随机变量序列YY的条件概率分布P(YX)P(Y|X)构成条件随机场,即满足马尔可夫性
P(YiX,Y1,...,Yi1,Yi+1,...,Yn)=P(YiX,Yi1,Yi+1)i=1,2,...,n(i=1n时只考虑单边)P(Y_i|X,Y_1,...,Y_{i-1},Y_{i+1},...,Y_n) = P(Y_i|X,Y_{i-1},Y_{i+1})\\ i=1,2,...,n(当i=1和n时只考虑单边)
则称P(YX)P(Y|X)为线性链条件随机场。

graph LR Y1(("Y₁")) Y2(("Y₂")) Yi(("Yᵢ")) Yn(("Yₙ")) Xg(("X₁:ₙ")) Y1---Y2-.-Yi-.-Yn Y1---Xg Y2---Xg Yi---Xg Xg---Yn style Y1 fill:#fff style Y2 fill:#fff style Yi fill:#fff style Yn fill:#fff style Xg fill:#f96

线性链条件随机场可以用于标注等问题。
在标注问题中,XX表示输入观测序列,YY表示对应的输出标记序列或状态序列。

这时,在条件概率模型P(Y|X)中,Y是输出变量,表示标记序列,X是输入变量,表示需要标注的观测序列。也把标记序列称为状态序列。
学习时,利用训练数据集通过极大似然估计或正则化的极大似然估计得到条件概率模型P^(YX)\hat{P}(Y|X)
预测时,对于给定的输入序列x,求出条件概率P^(YX)\hat{P}(Y|X)最大的输出序列y^\hat{y}

根据无向图的因子分解,得:
P(YX)=1Zexpi=1KFi(xci)P(Y|X) = \frac{1}{Z} exp\sum_{i=1}^K F_i(x_{ci})
根据CRF的概率无向图表示,我们可以将其实际的节点带入进去(最大团yt1,yt,x1:Ty_{t-1},y_t,x_{1:T}),(为了表达的方便,凑一个y0)有:
P(YX)=1Zexpt=1TF(yt1,yt,x1:T)P(Y|X) = \frac{1}{Z} exp\sum_{t=1}^T F(y_{t-1},y_t,x_{1:T})
F(yt1,yt,x1:T)F(y_{t-1},y_t,x_{1:T})分解为2个部分,即:x1:Tx_{1:T}(已知的)对yty_t的影响以及yt1,yty_{t-1},y_t间的影响。数学化表示为:
F(yt1,yt,x1:T)=yt,x1:T+yt1,yt,x1:TF(y_{t-1},y_t,x_{1:T})=\triangle y_{t},x_{1:T} + \triangle y_{t-1},y_{t},x_{1:T}
其实还有个yt1,x1:T\triangle y_{t-1},x_{1:T},因为这是上一个的状态,对于t时刻是已知的,这里忽略了。
其中,yt,x1:T\triangle y_t,x_{1:T}状态函数,即表示为在tt位置上的节点yty_t状态;
yt1,yt,x1:T\triangle y_{t-1},y_t,x_{1:T}转移函数,即表示当前节点yty_t与上一个节点yt1y_{t-1}的相关性。
定义在𝑌上下文的局部特征函数tkt_k,这类特征函数只和当前节点和上一个节点有关,即为上面的
yt1,yt,x1:T=k=1Kλktk(yi1,yi,X,i),k=1,2,..,K\triangle y_{t-1},y_t,x_{1:T} =\sum_{k=1}^K \lambda_k t_k(y_{i-1},y_i,X,i),k=1,2,..,K
其中KK是定义在该节点的局部特征函数的总个数,ii是当前节点在序列的位置。λk\lambda_k为特征函数的信任度。
定义在𝑌节点上的节点特征函数,这类特征函数只和当前节点有关,即为上面的
yt,x1:T=l=1Lμlsl(yi,X,i),l=1,2,,L\triangle y_t,x_{1:T} =\sum_{l=1}^L \mu_l s_l(y_i,X,i),l=1,2,…,L
其中LL是定义在该节点的节点特征函数的总个数,ii是当前节点在序列的位置。μl\mu_l为特征函数的信任度。
为了使特征函数make sense(有道理,合乎情理; 可以理解;讲得通),一般是指示函数,即取值非0即1。无论是节点特征函数还是局部特征函数,它们的取值只能是0或者1。即满足特征条件或者不满足特征条件。
如:
tk{yt1=名词,yt=动词,x1:T}=1tk{yt1=名词,yt=助词,x1:T}=0t_k\{y_{t-1}=名词, y_t=动词, x_{1:T}\} = 1 \\ t_k\{y_{t-1}=名词, y_t=助词, x_{1:T}\} = 0
所以linear-chain-CRF的参数化形式为:
P(YX)=1Z(X)expi=1T(k=1Kλktk(yi1,yi,X,i)+l=1Lμlsl(yi,X,i))P(Y|X)=\frac{1}{Z(X)}exp \sum_{i=1}^ T \bigg (\sum_{k=1}^K \lambda_k t_k (y_{i-1},y_i,X,i) +\sum_{l=1}^L \mu_l s_l (y_i,X,i)\bigg )
YY表示的是标注序列,是一个列向量,长度为TTX=x1:TX = x_{1:T}表示的词语序列,也是一个列向量,长度也为TT
其中,Z(X)Z(X)为规范化因子:
Z(X)=Yexpi=1T(k=1Kλktk(yi1,yi,X,i)+lLμlsl(yi,X,i))Z(X)=\sum_Y exp \sum_{i=1}^T \bigg(\sum_{k=1}^K\lambda_k t_k (y_{i-1},y_i,X,i) +\sum_{l}^L\mu_l s_l (y_i,X,i)\bigg)
模型的简化表示-数值表示
假设,共有K=K1+K2K=K_1+K_2个特征函数,其中,K1K_1个局部特征函数tkt_kK2K_2个节点特征函数sls_l。我们用1个特征函数fk(yi1,yi,X,i)f_k(y_{i-1},y_i,X,i)来统一表示:
fk(yi1,yi,X,i)={tk(yi1,yi,X,i)k=1,2,..,K1sl(yi,X,i)k=K1+l,l=1,2,,K2\begin{aligned}f_k(y_{i-1},y_i,X,i)=\left\{\begin{aligned} & t_k (y_{i-1},y_i,X,i) \qquad k = 1,2,..,K_1 \\ & s_l (y_i,X,i) \qquad k = K_1+l,l=1,2,…,K_2 \end{aligned}\right.\end{aligned}
fk(yi1,yi,X,i)f_k(y_{i-1},y_i,X,i)在各个序列位置求和得到:
fk(Y,X)=i=1Tfk(yi1,yi,X,i)\begin{aligned}f_k(Y,X)=\sum_{i=1}^T f_k(y_{i-1},y_i,X,i)\end{aligned}
同时也统一fk(yi1,yi,x,i)f_k(y_{i-1},y_i,x,i)对应的权重系数wkw_k
wk={λkk=1,2,..,K1μlk=K1+l,l=1,2,,K2\begin{aligned}w_k=\left\{ \begin{aligned} & \lambda_k \qquad k = 1,2,..,K_1 \\ & \mu_l \qquad k = K_1+l,l=1,2,…,K_2 \end{aligned}\right.\end{aligned}
这样,Linear-chain-CRF的简化工作就到这里结束啦:
P(YX)=1Z(X)expk=1Kwkfk(Y,X)\begin{aligned}P(Y|X)=\frac{1}{Z(X)}exp\sum_{k=1}^K w_kf_k(Y,X)\end{aligned}
其中,规范化因子:
Z(X)=Yexpk=1Kwkfk(Y,X)\begin{aligned}Z(X)=\sum_Y exp\sum_{k=1}^Kw_kf_k(Y,X)\end{aligned}
模型的简化表示-向量表示
如果对fk(Y,X)f_k(Y,X)wkw_k进行向量化表示,F(Y,X)F(Y,X)WW都是K×1K \times 1的列向量:
W=[w1w2wK]\begin{aligned}W =\left [ \begin{aligned} w_1\\ w_2\\ …\\ w_K \end{aligned}\right]\end{aligned}
F(Y,X)=[f1(Y,X)f2(Y,X)fK(Y,X)]\begin{aligned}F(Y,X) =\left[ \begin{aligned} f_1(Y,X)\\ f_2(Y,X)\\ ………\\ f_K(Y,X) \end{aligned}\right]\end{aligned}
那么Linear-chain-CRF的向量内积形式可以表示为:
PW(YX)=exp(WF(Y,X))Z(X,W)=exp(WF(Y,X))Yexp(WF(Y,X))\begin{aligned}P_W(Y|X) = \frac{exp(W \bullet F(Y,X))}{Z(X,W)}\\ = \frac{exp(W \bullet F(Y,X))}{\sum_Y exp(W \bullet F(Y,X))}\end{aligned}
向量化的意义:
就是为了干掉连加的形式,为后面的训练提供更加合理的计算支持。
要解决的三个问题

  1. Inference(概率计算问题):计算条件概率分布,即给定X序列,算出序列中每个位置所对应标注的概率,即:P(ytX)P(y_t|X)
  2. Learning:把参数学习出来(parameter estimation),也就是给定NN个训练数据,求上面向量表示的WW的参数值,即:W^=argmaxi=1NP(Y(i)X(i))\hat{W}=argmax\prod_{i=1}^N P(Y^{(i)}|X^{(i)})
  3. Decoding:给定X序列,找到一个最有可能的标注序列,即:Y^=argmaxP(YX)\hat{Y}=argmax P(Y|X),其中,Y=y1y2..yTY=y_1y_2..y_T

Inference:条件概率(前向-后向)
Learning(参数估计)
Decoding(Vitebi)

参考【NLP】从隐马尔科夫到条件随机场 以及视频机器学习-白板推导系列(十七)-条件随机场CRF(Conditional Random Field)

附加知识

随机场

随机场Random field, RF)是由若干个位置组成的整体,当给每一个位置中按照某种分布(或者是某种概率)随机赋予一个值之后,其全体就叫做随机场。

以词性标注为例:

假如我们有10个词形成的句子需要做词性标注。这10个词每个词的词性可以在我们已知的词性集合(名词,动词…)中去选择。当我们为每个词选择完词性后,这就形成了一个随机场。

马尔科夫随机场Markov random field, MRF)是随机场的特例,它假设随机场中某一个位置的赋值仅仅与和它相邻的位置的赋值有关,和与其不相邻的位置的赋值无关。
换一种表示方式,把马尔科夫随机场映射到无向图中。此无向图中的节点都与某个随机变量相关,连接着节点的边代表与这两个节点有关的随机变量之间的关系。
继续词性标注为例:(还是10个词的句子)
如果我们假设所有词的词性仅与和它相邻的词的词性有关时,这个随机场就特化成一个马尔科夫随机场。
比如第3个词的词性除了与自己本身的位置有关外,只与第2个词和第4个词的词性有关。

条件随机场(CRF)是马尔科夫随机场的特例,它假设马尔科夫随机场中只有𝑋和𝑌两种变量,𝑋一般是给定的,而𝑌一般是在给定𝑋的条件下我们的输出。这样马尔科夫随机场就特化成了条件随机场。

在我们10个词的句子词性标注的例子中,𝑋是词,𝑌是词性。因此,如果我们假设它是一个马尔科夫随机场,那么它也就是一个CRF。
对于CRF,我们给出准确的数学语言描述:
设𝑋与𝑌是随机变量,P(𝑌|𝑋)是给定𝑋时𝑌的条件概率分布,若随机变量𝑌构成的是一个马尔科夫随机场,则称条件概率分布P(𝑌|𝑋)是条件随机场。

线性链条件随机场(Linear-CRF)
注意在CRF的定义中,我们并没有要求𝑋和𝑌有相同的结构。当𝑋和𝑌有相同结构,即:
X=(x1,x2,,xT),Y=(y1,y2,,yT)X=(x_1,x_2,…,x_T),Y=(y_1,y_2,…,y_T)
这个时候,𝑋和𝑌有相同的结构的CRF就构成了线性链条件随机场。

MEMM(Maximum Entropy Markov Model)

判别模型
Maximum Entropy Markov Models for Information Extraction and Segmentation
Maximum Entropy Markov Models
Hidden Markov Model and Naive Bayes relationship
Maximum Entropy Markov Models and Logistic Regression

Maximum-Entropy Markov Model


MEMM与HMM

概率图模型

介绍概率图模型(Probabilistic Graphical Model)之前,先简单了解下结构学习Structured Learning),相比于回归,输出一个标量或者预测,输出一个向量,结构化学习的输出更加复杂,可以是图像,可以是语句,可以是树结构,等。
那么与概率图模型有什么关系呢?
概率图形模型形成了大量的结构化预测模型。特别是,贝叶斯网络和随机场很受欢迎。参见

什么是结构化学习?What is structured learning?
结构化预测是监督学习、分类和回归标准范式的概括。所有这些都可以被认为是找到一个函数来最小化训练集上的一些损失。区别在于使用的函数类型和损失。
在分类中,目标域是离散的类标签,损失通常是0-1的损失,即对误分类进行计数。在回归中,目标域是实数,损失通常是均方误差。在结构化预测中,目标域和损失或多或少都是任意的。这意味着目标不是预测标签或数字,而是可能更复杂的对象,如序列或图形。

概率图模型Probabilistic Graphical Model,PGM),简称图模型(Graphical Model,GM),是指一种用图结构来描述多元随机变量之间条件独立关系的概率模型,从而给研究高维空间中的概率模型带来了很大的便捷性。
很多机器学习模型都可以归结为概率模型,即建模输入和输出之间的条件概率分布.因此,图模型提供了一种新的角度来解释机器学习模型,并且这种角度有很多优点,比如了解不同机器学习模型之间的联系,方便设计新模型(Developing Bayesian networks)等.在机器学习中,图模型越来越多地用来设计和分析各种学习算法.

图模型有三个基本问题

  1. 表示(Representation)问题:对于一个概率模型,如何通过图结构来描述变量之间的依
    赖关系.
  2. 学习(Learning)问题:图模型的学习包括图结构的学习和参数的学习.在本章中,
    我们只关注在给定图结构时的参数学习,即参数估计问题.
  3. 推断(Inference)问题:在已知部分变量时,计算其他变量的条件概率分布

{Representation(表示){有向图 Bayesian Network无向图 Markov NetworkLearning(学习){参数学习{完备数据隐变量EM结构学习Inference(推断){精确推断近似推断{确定性近似变分推断随机近似MCMC\begin{cases} Representation(表示) & \begin{cases} \text{有向图 Bayesian Network} \\ \text{无向图 Markov Network} \end{cases} \\ Learning(学习) & \begin{cases} \text{参数学习} & \begin{cases} \text{完备数据} \\ \text{隐变量} \to EM \end{cases} \\ \text{结构学习} \end{cases}\\ Inference(推断) & \begin{cases} \text{精确推断} \\ \text{近似推断} & \begin{cases} \text{确定性近似} \to 变分推断 \\ \text{随机近似} \to MCMC \end{cases} \end{cases} \\ \end{cases}

图的表示
图可以用G=(V,E)G=(V,E)表示,VV是顶点vertices(nodes or points)集合,
E{(x,y)(x,y)V2  and  xy}{\displaystyle E\subseteq \{(x,y)\mid (x,y)\in V^{2}\;{\textrm {and}}\;x\neq y\}}是边的集合edges;对于有向图而言,边是有向的(directed edges, directed links, directed lines, arrows or arcs)它们是有序的顶点对,代表着方向;对于无向图而言,边是无向的。

也有些地方有向边一般用尖括号表示<>;而无向边一般用弧形括号表示();如:
有向图:
G1=(V,E)V(G1)={v1,v2,v3}E(G1)={v1,v2,v1,v3,v2,v3}G1=(V,E) \\ V(G1)=\{v1,v2,v3\}\\ E(G1)=\{\braket{v1,v2},\braket{v1,v3},\braket{v2,v3}\}

graph LR v1(("v1"))-->v2(("v2")) v1-->v3(("v3")) v2-->v3

无向图:
G2=(V,E)V(G2)={v1,v2,v3,v4}E(G2)={(vl,v2),(v1,v3),(v1,v4),(v2,v3),(v2,v4),(v3,v4)}G2=(V,E) \\ V(G2)=\{v1,v2,v3,v4\} \\ E(G2)=\{(vl,v2),(v1,v3),(v1,v4),(v2,v3),(v2,v4),(v3,v4)\}

graph LR v1(("v1"))---v2(("v2")) v1---v3(("v3")) v1---v4(("v4")) v2---v3 v2---v4 v3---v4
(概率)有向图模型

有向图模型(Directed Graphical Model)又称贝叶斯网络(Bayesian Network)或信念网络(Belief Network,BN)或(causal networks)是一类用有向图(Directed Graphical)来描述随机向量概率分布的模型.

这里是 有向无环图(DAG)

定义和概率 Definitions and concepts:

parent 父节点
descendants 后代
non-descendants 非后代(不包括父代,也就是all-parent-descendants)

graph LR x1(("x₁"))-->x2(("x₂"))-->x4(("x₄")) x1-->x3(("x₃")) x2-->x3 x3-->x5(("x₅"))

X=x1,x2,x3,x4,x5X=x_1,x_2,x_3,x_4,x_5
V={x1,x2,x3,x4,x5}V=\{x_1,x_2,x_3,x_4,x_5\}
E={x1,x2,x1,x3,x2,x3,x2,x4},x3,x5E=\{\braket{x_1,x_2},\braket{x_1,x_3},\braket{x_2,x_3},\braket{x_2,x_4}\},\braket{x_3,x_5}
G=(V,E)G=(V,E)
有向图对应的概率分布可以分解为
p(X)=p(x1,x2,x3,x4,x5)=p(x1)p(x2x1)p(x3x1,x2)p(x4x2)p(x5x3)p(X) = p(x_1,x_2,x_3,x_4,x_5) = p(x_1)p(x_2|x_1)p(x_3|x_1,x_2)p(x_4|x_2)p(x_5|x_3)

Pattern Model 条件独立性
Chain(间接因果关系/tail to head) XYZX\rightarrow Y\rightarrow Z 已知Y时,X和Z为条件独立,即 X ⁣ ⁣ ⁣ZYX \perp \!\!\!\perp Z\mid Y
Fork(共因关系/tail to tail) XYZX\leftarrow Y\rightarrow Z 已知Y时,X和Z为条件独立,即 X ⁣ ⁣ ⁣ZYX \perp \!\!\!\perp Z \mid Y (Y未知时,X和Z为不独立)
Collider(共果关系/head to head) XYZX\rightarrow Y\leftarrow Z 已知Y时,X和Z为不独立,即 X ⁣ ⁣ ⁣ ⁣ ⁣ ⁣ ⁣ ⁣ ⁣/    ZYX \perp \!\!\!\perp \!\!\!\!\!\!/ \;\; Z \mid Y(Y未知时,X和Z为独立)
graph TD y(("y")) y-->x1(("x₁")) y-->x2(("x₂")) y-->x3(("x₃")) style y fill:#fff style x1 fill:#f96 style x2 fill:#f96 style x3 fill:#f96

由图可得
P(y,x1,x2,x3)=P(y)P(x1y)P(x2y)P(x3y)=P(x1,x2,x3y)P(y)P(x1,x2,x3y)=P(x1y)P(x2y)P(x3y)P(y,x_1,x_2,x_3) = P(y)P(x_1|y)P(x_2|y)P(x_3|y) = P(x_1,x_2,x_3|y)P(y) \\ \Darr\\ P(x_1,x_2,x_3|y)=P(x_1|y)P(x_2|y)P(x_3|y)
这不就是朴素贝叶斯的条件相互独立的假设么?P(Xy)=i=1nP(xiy)P(X|y) = \prod_{i=1}^n P(x_i|y)
而这个独立假设太强了,每个特征之间没有任何关系(独立同分布i.i.d.);
那么我们假设当前只与前一时刻有关,与其它无关,那么我们就有了Markov假设,如隐马尔可夫模型:
其中y为隐变量,x为观测变量

graph LR y1(("y₁")) y1-->x1(("x₁")) y1-->y2(("y₂")) y2-->x2(("x₂")) y2-->y3(("y₃")) y3-->x3(("x₃")) y3-->y4(("y₄")) y4-->x4(("x₄")) style y1 fill:#fff style y2 fill:#fff style y3 fill:#fff style y4 fill:#fff style x1 fill:#f96 style x2 fill:#f96 style x3 fill:#f96 style x4 fill:#f96

我们能从图中直接得到
P(ytyt1,...,y1,xt1,...,x1)=P(ytyt1)P(y_t|y_{t-1},...,y_1,x_{t-1},...,x_1) = P(y_t|y_{t-1}),即Markov假设;
P(xtxT,...,xt+1,xt1,...,x1,Y)=P(xtyt)P(x_t|x_{T},...,x_{t+1},x_{t-1},...,x_1,Y) = P(x_t|y_{t}),即观测独立性假设;

(概率)无向图模型

无向图模型(Undirected Graphical Model)又称马尔可夫随机场(Markov random field, MRF)或马尔可夫网络(Markov network)是一类用无向图(Undirected Graphical)来描述一组具有局部马尔可夫性质的随机向量𝑿的联合概率分布的模型.
以下定义是等价的
Global Markov    Local Markov    Pair MarkovHammesleyClifford因子分解\text{Global Markov} \iff \text{Local Markov}\iff\text{Pair Markov}\xLeftrightarrow{Hammesley−Clifford } 因子分解

graph TD x1(("x₁")) x2(("x₂")) x3(("x₃")) x4(("x₄")) x1-->x4 x2-->x4 x3-->x4 y1(("x₁")) y2(("x₂")) y3(("x₃")) y4(("x₄")) y1---y2 y1---y3 y1---y4 y2---y3 y2---y4 y3---y4

道德化的过程中,原有的一些条件独立性会丢失。

以上内容只是讲到了概率图的表示。

参考文献

[11-1] Bishop M. Pattern Recognition and Machine Learning. Springer-Verlag,2006

[11-2] Koller D,Friedman N. Probabilistic Graphical Models: Principles and Techniques.MIT Press,2009

[11-3] Lafferty J,McCallum A,Pereira F. Conditional random fields: probabilistic models for segmenting and labeling sequence data. In: International Conference on Machine Learning,2001

[11-4] Sha F,Pereira F. Shallow parsing with conditional random fields. In: Proceedings ofthe 2003 Conference of the North American Chapter of Association for ComputationalLinguistics on Human Language Technology,Vol.1,2003

[11-5] McCallum A,Freitag D,Pereira F. Maximum entropy Markov models for informationextraction and segmentation. In: Proc of the International Conference on Machine Learning,2000

[11-6] Taskar B,Guestrin C,Koller D. Max-margin Markov networks. In: Proc of the NIPS2003,2003

[11-7] Tsochantaridis I,Hofmann T,Joachims T. Support vector machine learning forinterdependent and structured output spaces. In: ICML,2004

第 12 章 监督学习方法总结

参考:生成模型和判别模型

分类 方法 适用问题 模型特点 模型类别 学习策略 损失函数 学习算法
监督 感知机 二分类 分离超平面 判别模型 极小化误分点到超平面距离 误分点到超平面距离 随机梯度下降
监督 k 近邻法 多分类、回归 特征空间、样本点 判别模型
监督 朴素贝叶斯 多分类 特征与类别的联合概率分布,条件独立假设 生成模型 极大似然估计、最大后验概率估计 对数似然损失 概率计算公式、EM 算法
监督 决策树 多分类、回归 分类树、回归树 判别模型 正则化的极大似然估计 对数似然损失 特征选择、生成、剪枝
监督 逻辑斯蒂回归 多分类 特征条件下类别的条件概率分布,对数线性模型 判别模型 极大似然估计,正则化的极大似然估计 逻辑斯蒂损失 改进的迭代尺度算法,梯度下降,拟牛顿法
监督 支持向量机 二分类 分离超平面,核技巧 判别模型 极小化正则化合页损失,软间隔最大化 合页损失 序列最小最优化算法(SMO)
监督 提升方法 (Boosting) 二分类 弱分类器的线性组合 判别模型 极小化加法模型的指数损失 指数损失 前向分布加法算法
监督 EM 算法 概率模型参数估计 含隐变量概率模型 极大似然估计,最大后验概率估计 对数似然损失 迭代算法
监督 隐马尔科夫模型(HMM) 标注 观测序列与状态序列的联合概率分布模型 生成模型 极大似然估计,最大后验概率估计 对数似然损失 概率计算公式,EM 算法
监督 最大熵马尔科夫模型(MEMM) 标注 - 判别模型
监督 条件随机场(CRF) 标注 状态序列条件下观测序列的条件概率分布,对数线性模型 判别模型 极大似然估计,正则化极大似然估计 对数似然损失 改进的迭代尺度算法,梯度下降,拟牛顿法
监督 马尔可夫随机场 Markov Random Fields - - 生成模型

模型
分类问题与标注问题的预测模型都可以认为是表示从输入空间到输出空间的映射。它们可以写成条件概率分布P(Y|X)或决策函数Y=f(X)的形式。前者表示给定输入条件下输出的概率模型,后者表示输入到输出的非概率模型。有时,模型更直接地表示为概率模型,或者非概率模型;但有时模型兼有两种解释。

朴素贝叶斯法、隐马尔可夫模型是概率模型。感知机、k近邻法、支持向量机、提升方法是非概率模型。而决策树、逻辑斯谛回归与最大熵模型、条件随机场既可以看作是概率模型,又可以看作是非概率模型。

直接学习条件概率分布P(Y|X)或决策函数Y=f(X)的方法为判别方法,对应的模型是判别模型。感知机、k近邻法、决策树、逻辑斯谛回归与最大熵模型、支持向量机、提升方法、条件随机场是判别方法。首先学习联合概率分布P(X,Y),从而求得条件概率分布P(Y|X)的方法是生成方法,对应的模型是生成模型。朴素贝叶斯法、隐马尔可夫模型是生成方法。

可以用非监督学习的方法学习生成模型。具体地,应用EM算法可以学习朴素贝叶斯模型以及隐马尔可夫模型。

决策树是定义在一般的特征空间上的,可以含有连续变量或离散变量。感知机、支持向量机、k近邻法的特征空间是欧氏空间(更一般地,是希尔伯特空间)。提升方法的模型是弱分类器的线性组合,弱分类器的特征空间就是提升方法模型的特征空间。

感知机模型是线性模型,而逻辑斯谛回归与最大熵模型、条件随机场是对数线性模型。k近邻法、决策树、支持向量机(包含核函数)、提升方法使用的是非线性模型。

学习策略
在二类分类的监督学习中,支持向量机、逻辑斯谛回归与最大熵模型、提升方法各自使用合页损失函数、逻辑斯谛损失函数、指数损失函数。3种损失函数分别写为
[1yf(x)]+log[1exp(yf(x))]exp(yf(x))[1-yf(x)]_+ \\ \log[1-\exp(-yf(x))] \\ \exp(-yf(x))
这3种损失函数都是0-1损失函数的上界,具有相似的形状。所以,可以认为支持向量机、逻辑斯谛回归与最大熵模型、提升方法使用不同的代理损失函数(surrogate loss function)表示分类的损失,定义经验风险或结构风险函数,实现二类分类学习任务。学习的策略是优化以下结构风险函数:
minfH1Ni=1NL(yi,f(xi))+λJ(f)\min_{f \in H} \frac{1}{N}\sum_{i=1}^N L(y_i,f(x_i)) +\lambda J(f)
这里,第1项为经验风险(经验损失),第2项为正则化项,L(Y,f(X))为损失函数,J(f)为模型的复杂度,≥0为系数。

支持向量机用L2范数表示模型的复杂度。原始的逻辑斯谛回归与最大熵模型没有正则化项,可以给它们加上L2范数正则化项。提升方法没有显式的正则化项,通常通过早停止(early stopping)的方法达到正则化的效果。

以上二类分类的学习方法可以扩展到多类分类学习以及标注问题,比如标注问题的条件随机场可以看作是分类问题的最大熵模型的推广。

概率模型的学习可以形式化为极大似然估计或贝叶斯估计的极大后验概率估计。这时,学习的策略是极小化对数似然损失或极小化正则化的对数似然损失。对数似然损失可以写成
logP(yx)-\log P(y|x)
极大后验概率估计时,正则化项是先验概率的负对数。

决策树学习的策略是正则化的极大似然估计,损失函数是对数似然损失,正则化项是决策树的复杂度。

逻辑斯谛回归与最大熵模型、条件随机场的学习策略既可以看成是极大似然估计(或正则化的极大似然估计),又可以看成是极小化逻辑斯谛损失(或正则化的逻辑斯谛损失)。

朴素贝叶斯模型、隐马尔可夫模型的非监督学习也是极大似然估计或极大后验概率估计,但这时模型含有隐变量。

学习算法

统计学习的问题有了具体的形式以后,就变成了最优化问题。有时,最优化问题比较简单,解析解存在,最优解可以由公式简单计算。但在多数情况下,最优化问题没有解析解,需要用数值计算的方法或启发式的方法求解。

朴素贝叶斯法与隐马尔可夫模型的监督学习,最优解即极大似然估计值,可以由概率计算公式直接计算。

感知机、逻辑斯谛回归与最大熵模型、条件随机场的学习利用梯度下降法、拟牛顿法等。这些都是一般的无约束最优化问题的解法。

支持向量机学习,可以解凸二次规划的对偶问题。有序列最小最优化算法等方法。

决策树学习是基于启发式算法的典型例子。可以认为特征选择、生成、剪枝是启发式地进行正则化的极大似然估计。

提升方法利用学习的模型是加法模型、损失函数是指数损失函数的特点,启发式地从前向后逐步学习模型,以达到逼近优化目标函数的目的。

EM算法是一种迭代的求解含隐变量概率模型参数的方法,它的收敛性可以保证,但是不能保证收敛到全局最优。

支持向量机学习、逻辑斯谛回归与最大熵模型学习、条件随机场学习是凸优化问题,全局最优解保证存在。而其他学习问题则不是凸优化问题。

第 13 章 无监督学习概论

基本问题
聚类(clustering)是将样本集合中相似的样本归到相同的类,相似的定义一般用距离度量。
如果一个样本只能属于一个类,则称为硬聚类(hard clustering),如k-means;如果一个样本可以属于多个类,每一个样本以概率属于每一个类i=1Np(zixi)=1\sum_{i=1}^N p(z_i|x_i) =1,则称为软聚类(soft clustering),如GMM。
聚类主要用于数据分析,也可用于监督学习的前处理。聚类可以帮助发现数据中的统计规律。

降维(dimensionality reduction)是将样本集合中的样本(实例)从高维空间转换到低维空间。
高维空间通常是高维的欧氏空间,而低维空间是低维的欧氏空间或流形(manifold)。低维空间是从数据中自动发现的。降维有线性降维和非线性降维,降维方法有主成分分析。
降维的好处有:节省存储空间、加速计算、解决维度灾难(前面章节有讲到)等
降维主要用于数据分析,也可用于监督学习的前处理。降维可以帮助发现高维数据中的统计规律。

概率模型估计(probability model estimation),简称概率估计,假设训练数据由一个概率模型生成,由训练数据学习概率模型的结构和参数。
概率模型的结构类型或者概率模型的集合事先给定,而模型的具体结构与参数从数据中自动学习。假设数据由GMM生成(已知结构),学习的目标是估计这个模型的参数。
概率模型包括混合模型、概率图模型等。概率图模型又包括有向图模型和无向图模型(前面章节有讲到)。

无监督学习方法

同监督学习一样,无监督学习也有三要素:模型、策略、算法
模型 就是函数z=gθ(x)z=g_\theta(x),条件概率分布Pθ(zx)P_\theta(z |x),或Pθ(xz)P_\theta(x|z),在聚类、降维、概率模型估计中拥有不同的形式

策略 在不同的问题中有不同的形式,但都可以表示为目标函数的优化

算法 通常是迭代算法,通过迭代达到目标函数的最优化,比如,梯度下降法。

参考文献

[13-1] Hastie T, Tibshirani R, Friedman J. The elements of statistical learning:data mining, inference, and prediction. Springer. 2001.

[13-2] Bishop M. Pattern Recognition and Machine Learning. Springer, 2006.

[13-3] Koller D, Friedman N. Probabilistic graphical models: principles and techniques. Cambridge, MA: MIT Press, 2009.

[13-4] Goodfellow I,Bengio Y,Courville A. Deep learning. Cambridge, MA: MIT Press, 2016.

[13-5] Michelle T M. Machine Learning. McGraw-Hill Companies, Inc. 1997.(中译本:机器学习。北京:机械工业出版社,2003.)

[13-6] Barber D. Bayesian reasoning and machine learning, Cambridge, UK:Cambridge University Press, 2012.

[13-7] 周志华. 机器学习. 北京:清华大学出版社,2017.

第 14 章 聚类方法

聚类分析(Cluster analysis or clustering)是针对给定的样本,根据它们特征点的相似度(距离),将其归并到若干个类(簇)中的分析问题。聚类分析本身不是一种特定的算法,而是要解决的一般任务。

算法分类

Hard clustering
Soft clustering (also: fuzzy clustering)


Connectivity models:如 hierarchical clustering 基于距离连通性(based on distance connectivity)
Centroid models:如 k-means
Distribution models:如GMM
Density models:如 DBSCAN and OPTICS
Subspace models:如 biclustering
Group models:如
Graph-based models:如
Signed graph models:如
Neural models:如

聚类分析算法(Cluster analysis algorithms):sklearn中的聚类算法和介绍
基于连接的聚类(层次聚类) Connectivity-based clustering(Hierarchical clustering):
层次聚类有聚合Agglomerative(自下而上"bottom-up")和分裂 Divisive(自上而下"top-down")两种方法。

  1. 距离度量或相似度
  2. 合并规则(cluster间的距离规则,cluster间的距离可以是最短、最长、中心距离、平均距离等)
  3. 停止条件(达到k值,也就是cluster的个数达到阈值;cluster的直径达到阈值;)

基于质心的聚类 Centroid-based clustering:

基于分布的聚类 Distribution-based clustering:

基于密度的聚类 Density-based clustering:

基于网格的聚类 Grid-based clustering:
其它聚类

高效聚类

高维数据的聚类Clustering high-dimensional data):
子空间聚类Subspace clustering

投影聚类Projected clustering
基于投影的聚类Projection-based clustering
混合方法Hybrid approaches

算法评估Evaluation

聚类结果的评估与聚类本身一样困难(并不像计算错误数量或监督分类算法的精度和召回率那么简单)。一般分为Internal evaluation和External evaluation,当两种评估效果都不好时就需要human evaluation。


先作一些定义:

数据集D={x1,x2,...,xi,...,xj,...,xn}D=\{x_1,x_2,...,x_i,...,x_j,...,x_n\}

TT为给定的正数,若集合CpC_p中任意两个样本间的距离dist(xi,xj)Tdist(x_i,x_j) \leq T,则称CpC_p为一个类或簇(cluster)

C={C1,C2,...,Ck}C = \{C_1,C_2,...,C_k\}表示(预测的)类或簇(cluster)
C={C1,C2,...,Cs}C^* = \{C_1^*,C_2^*,...,C_s^*\}表示参考模型的类或簇(cluster)
λ\lambda表示簇CC的标记(预测)向量,如:λ=[0,1,...,k],λ=[0,1,...,s]\lambda = [0,1,...,k],\lambda^* = [0,1,...,s],长度为样本数量nn
λi\lambda_i为样本xix_i的预测或标记值
a=TP,TP={(xi,xj)λi=λj,λi=λj,i<j}a = TP, TP=\{(x_i,x_j)\mid\lambda_i = \lambda_j, \lambda_i^* = \lambda_j^* ,i \lt j\} ,表示“样本对”在CC中属于相同的簇且在CC^*中也属于相同的簇的数量(true positive)
b=TN,TN={(xi,xj)λi=λj,λiλj,i<j}b = TN, TN=\{(x_i,x_j)\mid\lambda_i = \lambda_j, \lambda_i^* \neq \lambda_j^* ,i \lt j\} ,表示“样本对”在CC中属于相同的簇且在CC^*中也属于不同的簇的数量(true negative)
c=FP,FP={(xi,xj)λiλj,λi=λj,i<j}c = FP, FP=\{(x_i,x_j)\mid\lambda_i \neq \lambda_j, \lambda_i^* = \lambda_j^* ,i\lt j\} ,表示“样本对”在CC中属于不同的簇且在CC^*中也属于相同的簇的数量(false positive)
d=FN,FN={(xi,xj)λiλj,λiλj,i<j}d = FN, FN=\{(x_i,x_j) \mid \lambda_i \neq \lambda_j, \lambda_i^* \neq \lambda_j^* ,i\lt j\} ,表示“样本对”在CC中属于不同的簇且在CC^*中也属于不同的簇的数量(false negative)

注意:labels_pred = [0, 0, 1, 1] 与 labels_true = [0, 0, 1, 1] 以及 labels_pred = [0, 0, 1, 1] 与 labels_true = [1, 1, 0, 0] 是没有区别的,他们都正确的聚类了。

样本对的数量为Cn2=(n2)=n(n1)/2=a+b+c+dC_n^2 = \binom{n}{2} =n(n-1)/2 = a+b+c+d,这里的CC是排列组合的意思

dij=dist(xi,xj)d_{ij} = dist(x_i,x_j)表示样本xi,xjx_i,x_j之间的距离

np=Cpn_p = |C_p|表示簇CpC_p中的样本数量,
xˉp=1npxiCpxi\bar{x}_p = \frac{1}{n_p}\sum_{x_i \in C_p}x_i分别表示簇CpC_p的质心(中心、均值、中心点、centroid)
diam(Cp)=max{dist(xi,xj)xi,xjCp}diam(C_p) = \max \{dist(x_i,x_j) \mid x_i,x_j \in C_p\}表示簇的直径diam或者簇类样本间的最远距离
avg(Cp)=2np(np1)1i<jnpdist(xi,xj)avg(C_p) = \frac{2}{n_p(n_p-1)} \sum_{1\leq i \lt j\leq n_p }dist(x_i,x_j)表示簇类样本间的平均距离
ACp=xiCp(xixˉp)(xixˉp)TA_{C_p} = \sum_{x_i \in C_p} (x_i-\bar{x}_p)(x_i-\bar{x}_p)^T表示簇类样本散布矩阵(scatter matrix
SCp=1np1ACpS_{C_p} = \frac{1}{n_p-1}A_{C_p}表示簇类样本协方差矩阵(covariance matrix

两个簇之间的距离主要有以下表示方法:
dmin(Cp,Cq)=min{dist(xi,xj)xiCp,xjCq}d_{min}(C_p,C_q) = \min\{dist(x_i,x_j) \mid x_i \in C_p,x_j \in C_q\}表示两个簇间的最短距离
dmax(Cp,Cq)=max{dist(xi,xj)xiCp,xjCq}d_{max}(C_p,C_q) = \max\{dist(x_i,x_j) \mid x_i \in C_p,x_j \in C_q\}表示两个簇间的最长距离
dcen(Cp,Cq)=dist(xˉp,xˉq)d_{cen}(C_p,C_q) = dist(\bar{x}_p,\bar{x}_q)表示两个簇中心间的距离
dmean(Cp,Cq)=1npnqxiGpxjGqdist(xi,xj)d_{mean}(C_p,C_q) = \frac{1}{n_p n_q}\sum_{x_i \in G_p}\sum_{x_j \in G_q} dist(x_i,x_j)表示两个簇 任意两个样本之间距离的平均值 为两个簇的距离


聚类标准(Clustering criteria):簇内相似度(intra-cluster similarity)高,簇间相似度(inter-cluster similarity)低

确定数据集中的簇数Determining the number of clusters in a data set),也就是K值的选取。
对于某类聚类算法(特别是k-means, k-medoids),有一个通常称为k的参数指定要检测的聚类数。其他算法如DBSCAN和OPTICS 算法不需要指定该参数;层次聚类完全避免了这个问题。
簇数是数据集中重要的概括统计量。经验值:kn/2k \thickapprox \sqrt{n/2}

Internal evaluation:
无监督的方法,无需基准数据。类内聚集程度和类间离散程度。簇内相似度(intra-cluster similarity)高,簇间相似度(inter-cluster similarity)低。

metrics.davies_bouldin_score

metrics.silhouette_score

metrics.calinski_harabasz_score

External evaluation(就是需要人为标记每个样本所属的类)
有监督的方法,需要基准数据或者参考模型。用一定的度量评判聚类结果与基准数据的符合程度。(基准是一种理想的聚类,通常由专家构建)

Precision and recall
metrics.f1_score

metrics.jaccard_score

metrics.homogeneity_score, metrics.completeness_score, metrics.v_measure_score
三个指数一起返回metrics.homogeneity_completeness_v_measure

衡量两个数据聚类之间相似性的度量(也可以是标记数据和预测数据之间的相似性)
1.0 是完美匹配分数。未调整 Rand 指数的得分范围为 [0, 1],调整后(adjusted)的 Rand 指数为 [-1, 1]
metrics.rand_score 、metrics.adjusted_rand_score

metrics.fowlkes_mallows_score

from sklearn.metrics.cluster import contingency_matrix
x = ["a", "a", "a", "b", "b", "b"]
y = [0, 0, 1, 1, 2, 2]
contingency_matrix(x, y)
array([[2, 1, 0],
       [0, 1, 2]])

列数代表y中不重复的个数;行代表x中不重复的个数;
第一行分别表示:2个a属于0,1个a属于1,0个a属于2
第二行分别表示:0个b属于0,1个b属于1,2个b属于2

metrics.cluster.contingency_matrix

行Predicted (expectation)
列Actual (observation)
TN表示预测negative正确
TP表示预测positive正确
FN表示预测negative错误
FP表示预测positive错误

Actual\Predicted P N
P(positive) TP FN
N(negative) FP TN

https://en.jinzhao.wiki/wiki/Template:Diagnostic_testing_diagram

注意正常的confusion matrix中的四个元素相加为Cn2=n(n1)/2C_n^2=n(n-1)/2,而pair_confusion_matrix是(n1)(n-1),并且混淆矩阵
C=[C00(FN)C01(FP)C10(TN)C11(TP)]\begin{split}C = \left[\begin{matrix} C_{00}(FN) & C_{01}(FP) \\ C_{10}(TN) & C_{11}(TP) \end{matrix}\right]\end{split}
metrics.cluster.pair_confusion_matrix
metrics.plot_confusion_matrix可以绘制混淆矩阵

参考Adjusted mutual information
metrics.adjusted_mutual_info_score,metrics.normalized_mutual_info_score,metrics.mutual_info_score

Cluster tendency(聚类趋势)

参考Clustering performance evaluation以及Evaluation and assessment

k-means: 样本集合X={x1,x2,...,xn},xiRmX=\{x_1,x_2,...,x_n\},x_i \in \mathbb{R}^m,算法的目标是将n个样本分到不同的cluster中C={C1,...,Ck},k<n,CiCj=,i=1kCi=XC = \{ C_1,...,C_k\},k \lt n,C_i \cap C_j =\empty , \cup_{i=1}^kC_i =X;
F:xil,l{1,...,k}F: x_i \to l,l\in \{1,...,k\}表示划分函数,输入样本,输出所在的cluster

n个样本分到k个cluster的所有分法的种类有S(n,k)S(n,k)种,这个数字是指数级的,所以最优问题是个NP困难问题
S(n,k)=1k!l=1k(1)kl(kl)knS(n,k) = \frac{1}{k!}\sum_{l=1}^k(-1)^{k-l}\dbinom{k}{l}k^n

  1. 随机选择k个中心(m1,m2,...,mk)(m_1,m_2,...,m_k),
  2. 将样本分别划分到与其最近的中心的cluster中
  3. 更新每个cluster的均值(m1,m2,...,mk)(m_1,m_2,...,m_k)作为cluster的新的中心
  4. 重复2,3直到收敛(中心变化很小)

K-means算法有以下不足:

  1. 算法对初始值的选取依赖性极大。初始值不同,往往得到不同的局部极小值。
  2. 由于将均值点作为聚类中心进行新一轮计算,远离数据密集区的孤立点和噪声点会导致聚类中心偏离真正的数据密集区,所以K-均值算法对噪声点和孤立点很敏感。

K-mediods算法优缺点
K-中心点轮换算法是一种使目标函数下降最快的方法,它属于启发式搜索算法,能从n个对象中找出以k个中心点为代表的一个局部优化划分聚类。与K-均值算法比较,K-中心点轮换算法解决了K-均值算法本身的缺陷:

  1. 解决了K-均值算法对初始值选择依赖度大的问题。K-均值算法对于不同的初始值,结果往往得到不同的局部极小值。而K-中心点轮换算法采用轮换替换的方法替换中心点,从而与初始值的选择没有关系。
  2. 解决了K-均值算法对噪声和离群点的敏感性问题。由于该算法不使用平均值来更改中心点而是选用位置最靠近中心的对象作为中心代表点,因此并不容易受极端数据的影响,具有很好的鲁棒性。

K-中心点轮换算法也存有以下缺点:

  1. 由于K-中心点轮换算法是基于划分的一种聚类算法,仍然要求输入要得到的簇的数目k,所以当k的取值不正确时,对聚类的结果影响甚大。
  2. 从以上的时间复杂度也可以看出,当n和k较大时,计算代价很高,所以将该算法应用于大数据集时不是很理想。

参考文献

[14-1] Jain A, Dubes R. Algorithms for clustering data. Prentice-Hall, 1988.

[14-2] Aggarwal C C, Reddy C K. Data clustering: algorithms and applications. CRC Press, 2013.

[14-3] MacQueen J B. Some methods for classification and analysis of multivariate observations. Proceedings of the 5th Berkeley Symposium on Mathematical Statistics and Probability. Volume 1,pp.396-410. 1967.

[14-4] Hastie T,Tibshirani R,Friedman J. The Elements of Statistical Learning: DataMining,Inference,and Prediction. Springer. 2001(中译本:统计学习基础——数据挖掘、推理与预测。范明,柴玉梅,昝红英等译。北京:电子工业出版社,2004)

[14-5] Pelleg D, Moore A W. X-means: Extending K-means with Efficient Estimation of the Number of Clusters. Proceedings of ICML, pp. 727-734, 2000.

[14-6] Ester M, Kriegel H, Sander J, et al. A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise. Proceedings of ACM SIGKDD, pp. 226-231, 1996.

[14-7] Shi J, Malik J. Normalized cuts and image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2000,22(8):888-905.

[14-8] Dhillon I S. Co-clustering documents and words using bipartite spectral graph partitioning. Proceedings of ACM SIGKDD, pp. 269-274. 2001.

第 15 章 奇异值分解

奇异值分解(Singular Value Decomposition, SVD)是在机器学习领域广泛应用的算法,它不光可以用于降维算法中的特征分解,还可以用于推荐系统,以及自然语言处理等领域。是很多机器学习算法的基石。也是矩阵分解(Matrix decomposition)的一种。

先了解下特征值分解(Eigenvalues and eigenvectors以及Eigendecomposition of a matrix)以及对角化(Diagonalizable matrix
特征值(有些方阵是没有特征值的):
Au=λu    (AIλ)u=0{\displaystyle A\mathbf {u} =\lambda \mathbf {u} \implies (A-I\lambda)\mathbf {u} = 0}
特征值分解:
An×nA_{n\times n},是具有n个线性无关的特征向量qiq_i(不一定是不同特征值,可以有相同的特征值),那么A可以分解为:
A=QΛQ1{\displaystyle A=Q\Lambda Q^{-1}}
其中A是方阵,Λ\Lambda是由特征值组成的对角矩阵(diagonal matrix);
qiq_i通常是标准化的,但不是必须的;
因为Q的列是线性无关的,所以 Q 是可逆的;
如果A的特征值都不为0那么A是可逆的(invertibleA1=QΛ1Q1{\displaystyle \mathbf {A} ^{-1}=\mathbf {Q} \mathbf {\Lambda } ^{-1}\mathbf {Q} ^{-1}}

注意只有可对角化的矩阵分解才能用这种方式:如以下矩阵不可被对角化
A=[1101]{\displaystyle A = \left[{\begin{matrix}1&1\\0&1\end{matrix}}\right]}
其特征值为[1,1][1,1],特征向量为[1,0]T,[1,0]T[1,0]^T , [-1,0]^T

如果 A\mathbf {A} 是对称矩阵(symmetric matrix),因为Q\mathbf {Q} 由特征向量构成 A\mathbf {A} 它保证是一个正交矩阵(orthogonal matrix),有Q1=QT{\displaystyle \mathbf {Q} ^{-1}=\mathbf {Q} ^{\mathrm {T} }}

Q其实也是酉矩阵(Unitary Matrix),它是 正交矩阵 在复数上的推广。

不可对角化的矩阵称为有缺陷的(defective)。对于有缺陷的矩阵,特征向量的概念推广到广义特征向量(generalized eigenvectors),特征值的对角矩阵推广到Jordan 范式(Jordan normal form)。在代数闭域上,任何矩阵A都具有Jordan 范式,因此允许广义特征向量的基和分解为广义特征空间(generalized eigenspaces)。

附加知识

矩阵性质

这里介绍一些参见的矩阵,以及其性质。

共轭转置(Conjugate transpose

共轭(Complex conjugate)是复数上的概念
对于一个复数z=a+biz =a + bi,其共轭为zˉ=abi\bar{z} = a-bi,所以有zzˉ=a2+b2z\bar{z} = a^2 + b^2
共轭转置也有其它叫法,如:Hermitian conjugate, bedaggered matrix, adjoint matrix or transjugate。值得注意的是adjoint matrix而不是 这个Adjugate matrix,虽然有时候他们都用AA^*表示。这里为了统一我用AHA^H表示A的共轭转置矩阵。可以参考共轭转置矩阵与伴随矩阵都用A*表示合理吗?

有个神奇的公式:欧拉公式 eπi+1=0e^{\pi i}+1=0

(AH)ij=Aji{\displaystyle \left({\boldsymbol {A}}^{\mathrm {H} }\right)_{ij}={\overline {{\boldsymbol {A}}_{ji}}}}

AH=(A)T=AT{\displaystyle {\boldsymbol {A}}^{\mathrm {H} }=\left({\overline {\boldsymbol {A}}}\right)^{\mathsf {T}}={\overline {{\boldsymbol {A}}^{\mathsf {T}}}}}

例如:
A=[12i51+ii42i],AT=[11+i2ii542i],AH=[11i2+ii54+2i]{\displaystyle {\boldsymbol {A}}={\begin{bmatrix}1&-2-i&5\\1+i&i&4-2i\end{bmatrix}}} , {\displaystyle {\boldsymbol {A}}^{\mathsf {T}}={\begin{bmatrix}1&1+i\\-2-i&i\\5&4-2i\end{bmatrix}}} , {\displaystyle {\boldsymbol {A}}^{\mathrm {H} }={\begin{bmatrix}1&1-i\\-2+i&-i\\5&4+2i\end{bmatrix}}}

性质:

埃尔米特矩阵(Hermitian matrix

Hermitian matrix (or self-adjoint matrix)
A是复方阵
A Hermitian    A=AH{\displaystyle A{\text{ Hermitian}}\quad \iff \quad A=A^{\mathsf {H}}}
例子:
A=[0aibcida+ib1minc+idm+in2]A = {\displaystyle {\begin{bmatrix}0&a-ib&c-id\\a+ib&1&m-in\\c+id&m+in&2\end{bmatrix}}}
埃尔米特矩阵是对称矩阵在复数上的推广

其矩阵有很多性质,具体见维基百科。

Skew-Hermitian matrix:A skew-Hermitian    AH=A{\displaystyle A{\text{ skew-Hermitian}}\quad \iff \quad A^{\mathsf {H}}=-A}

Normal matrix

A是复方阵
A normal    AHA=AAH{\displaystyle A{\text{ normal}}\quad \iff \quad A^{H}A=AA^{H}}

例子:
A=[110011101],AAH=[211121112]=AHA.{\displaystyle A={\begin{bmatrix}1&1&0\\0&1&1\\1&0&1\end{bmatrix}}} , {\displaystyle AA^{H}={\begin{bmatrix}2&1&1\\1&2&1\\1&1&2\end{bmatrix}}=A^{H}A.}

Normal matrix一定是可对角化的A=UΛUHA = U\Lambda U^HUU是酉矩阵,Λ=diag(λ1,...)\Lambda = diag(\lambda_1,...)AA的特征值组成的对角矩阵

对于复矩阵,所有的unitary, Hermitian, and skew-Hermitian 矩阵都是normal矩阵
对应的对于实矩阵,所有的 orthogonal, symmetric, and skew-symmetric 矩阵也都是normal矩阵

酉矩阵(Unitary matrix

U是复方阵
UH=U1U^{H} = U^{-1}
性质:

酉矩阵它是正交矩阵在复数上的推广

矩阵分解(因子分解)

sympy.Matrix除了分解还有diagonalize对角化(也是一种矩阵分解),eig特征值(其实也可以特征值分解),rref行简化阶梯型,det行列式,inv逆矩阵,广义逆矩阵pinv;更多参考
scipy.linalg中也有很多关于线性代数的方法:scipy.linalg.svd,以及各种矩阵分解的方法;更多参考
numpy.linalg中也有很多关于线性代数的方法:np.linalg.svd;更多参考

除了SVD和PCA,还有很多矩阵分解(Matrix decomposition)的方法。不过有很多分解是有要求的,比如必须是方阵(特征值分解就要求必须是方阵)等。

类似的可以定义QL、RQ 和 LQ,L是下三角矩阵(left(lower) triangular matrix

参考文献

[15-1] Leon S J. Linear algebra with applications. Pearson, 2009(中译本:线性代数。张文博,张丽静 译. 北京:机械工业出版社)

[15-2] Strang G. Introduction to linear algebra. Fourth Edition. Wellesley-Cambridge Press, 2009.

[15-3] Cline A K. Dhillon I S. Computation of the singular value decomposition, Handbook of linear algebra. CRC Press, 2006.

[15-4] 徐树方. 矩阵计算的理论与方法。北京:北京大学出版社, 1995.

[15-5] Kolda T G,Bader B W. Tensor decompositions and applications. SIAM Review, 2009, 51(3):455-500.

第 16 章 主成分分析

主成分分析(Principal component analysis, PCA)是一种常用的无监督学习方法,PCA利用正交变换把由线性相关变量表示的观测数据转换为少数几个由线性无关变量表示的数据,线性无关的变量称为主成分。

主成分分析步骤如下:

  1. 对给定数据进行规范化,使得数据每一变量的平均值为0,方差为1(StandardScaler)。

  2. 对数据进行正交变换,原来由线性相关变量表示的数据,通过正交变换变成由若干个线性无关的新变量表示的数据。

新变量是可能的正交变换中变量的方差的和(信息保存)最大的,方差表示在新变量上信息的大小。将新变量依次称为第一主成分、第二主成分等。

我们通常表示一个样本是在实数空间中用正交坐标系表示,规范化的数据分布在原点附近

主成分分析就是对数据进行正交变换,对原坐标系进行旋转变换,并将数据在新坐标系中表示;我们将选择方差最大的方向作为新坐标系的第一坐标轴。方差最大代表着在该方向上的投影(不就是在这个坐标系的坐标轴上的表示么)分散的最开。

根据方差的定义,每个方向上方差就是在该坐标系(变换后的新坐标系)上表示所对应的维度的方差var(a)=1N1i=1N(aiμ)2var(a) = \frac{1}{N-1}\sum_{i=1}^N (a_i - \mu)^2(用第一个方向来说明, N个样本的第一个维度组成向量aa);由于我们已经对数据进行的规范化,所以均值为0;var(a)=1N1i=1N(ai)2var(a) = \frac{1}{N-1}\sum_{i=1}^N (a_i)^2 ;aia_i就是第ii个样本x(i)x^{(i)}与第一个方向的内积。

我们的目的就是为了var(a)var(a)最大,我们要求的就是找到变换后的新坐标系,假设方差最大的方向的单位向量为v1v_1,数据集T={x(1),x(2),...,x(N)},x={x1,...,xm}TT = \{x^{(1)},x^{(2)},...,x^{(N)}\} , x=\{x_1,...,x_m\}^T,m维

max1N1i=1Nx(i),v12=1N1i=1Nx(i)T.v12=1N1i=1N(x(i)T.v1)T(x(i)T.v1)=1N1i=1Nv1Tx(i)x(i)Tv1=1N1v1Ti=1N[x(i)x(i)T]v1\max \frac{1}{N-1}\sum_{i=1}^N \braket{x^{(i)},v_1}^2 = \frac{1}{N-1}\sum_{i=1}^N \|{x^{(i)}}^{T}.v_1\|^2 \\= \frac{1}{N-1}\sum_{i=1}^N ({x^{(i)}}^{T}.v_1)^T({x^{(i)}}^{T}.v_1) \\= \frac{1}{N-1}\sum_{i=1}^N v_1^T{x^{(i)}}{x^{(i)}}^{T}v_1 \\= \frac{1}{N-1} v_1^T \sum_{i=1}^N[{x^{(i)}}{x^{(i)}}^{T}]v_1
设矩阵X=[x(1),x(2),...,x(N)]X = [x^{(1)},x^{(2)},...,x^{(N)}]那么XXT=i=1N[x(i)x(i)T]XX^T =\sum_{i=1}^N[{x^{(i)}}{x^{(i)}}^{T}],得到
max1N1v1TXXTv1s.t.v1Tv1=1\max \frac{1}{N-1} v_1^T XX^T v_1 \\ s.t. \quad v_1^Tv_1 =1
拉格朗日函数(参见矩阵求导
L=1N1v1TXXTv1+λ1(v1Tv11)L = - \frac{1}{N-1} v_1^T XX^T v_1 + \lambda_1(v_1^Tv_1 - 1)
Lv1=21N1v1TXXT+2λ1v1T=0\frac{\partial L}{\partial v_1} = -2\frac{1}{N-1}v_1^T XX^T + 2\lambda_1 v_1^T =0
1N1v1TXXT=λ1v1T    1N1XXTv1=λ1v1\frac{1}{N-1}v_1^T XX^T = \lambda_1 v_1^T \implies \frac{1}{N-1}XX^Tv_1 = \lambda_1 v_1
其实1N1XXT\frac{1}{N-1}XX^T就是Xm×NX_{m \times N}样本的协方差矩阵Σm×m\Sigma_{m \times m}

λ1\lambda_1Σm×m\Sigma_{m \times m}的特征值,v1v_1(列向量)是其对应的特征值向量;

接着求第二个主成分v2v_2,主成分是相互正交的
max1N1v2TXXTv2s.t.v2Tv2=1,v2Tv1=0\max \frac{1}{N-1} v_2^T XX^T v_2 \\ s.t. \quad v_2^Tv_2 =1 ,v_2^Tv_1 =0
注意到
v1TXXTv2=λ1v1Tv2=0=v2TXXTv1=λ1v2Tv1v_1^T XX^T v_2 = \lambda_1 v_1^T v_2 = 0 = v_2^T XX^T v_1 = \lambda_1 v_2^T v_1
依次求得其它成分。

最终有主成分组成的矩阵
Vm×m=[v1,v2,...,vm]V_{m \times m } = [v_1,v_2,...,v_m]
降维到k维就是一次取前k个向量组成的矩阵与X作乘积,那么降维后的数据:

Yk×N=Vm×kTXm×NY_{k \times N} = V_{m \times k }^T X_{m \times N}

前面学习了SVD需要求ATAA^TA的特征值分解;而PCA需要求1N1XXT\frac{1}{N-1}XX^T的特征值分解;
只需要取A=XTN1A = \frac{X^T}{\sqrt{N-1}}就可以将PCA问题可以转化为SVD问题求解
其实,PCA只与SVD的右奇异向量的压缩效果相同。
一般 XX 的维度很高,XXTXX^T 的计算量很大,并且方阵的特征值分解计算效率不高,SVD除了特征值分解这种求解方式外,还有更高效且更准确的迭代求解法,避免了XXTXX^T的计算

稀疏主成分分析Sparse PCA
稀疏 PCA 问题有许多不同的公式,以下是使用Structured Sparse Principal Component Analysis以及Online Dictionary Learning for Sparse Coding

(U,V)=arg min U,V12XUV22+αV1subject to Uk2=1 for all 0k<ncomponents\begin{split}(U^*, V^*) = \underset{U, V}{\operatorname{arg\,min\,}} & \frac{1}{2} ||X-UV||_2^2+\alpha||V||_1 \\ \text{subject to } & ||U_k||_2 = 1 \text{ for all } 0 \leq k < n_{components}\end{split}

意思就是求UV让其近似等于X,然后得到一个稀疏矩阵V
sklearn.decomposition.SparsePCA.components_ 就是其稀疏的矩阵VV
SPCA的含义参考 Matrix decomposition

附加知识

基变换

我们常说的向量(3,2)其实隐式引入了一个定义:以 x 轴和 y 轴上正方向长度为 1 的向量为标准。向量 (3,2) 实际是说在 x 轴投影为 3 而 y 轴的投影为 2。注意投影是一个标量,所以可以为负。
而x 轴和 y 轴的方向的单位向量就是(1,0)和(0,1),所以(1,0)和(0,1)就是坐标系的基

如:另一组基(单位向量)(12,12)(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})(12,12)(-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})
那么(3,2)在该坐标系中如何表示呢?我们知道一个向量aa在一个方向(单位向量x单位向量x)上的投影可以用内积表示a,x=a.xcosθ=acosθ\braket{a,x} = \|a\|.\|x\|\cos \theta = \|a\|\cos \theta,其中θ\theta表示两个向量的夹角。

a=(3,2)a=(3,2)x=(12,12)x=(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})这个方向的投影为a,x=52\braket{a,x} = \frac{5}{\sqrt{2}}
a=(3,2)a=(3,2)y=(12,12)y=(-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})这个方向的投影为a,y=12\braket{a,y} = -\frac{1}{\sqrt{2}}

所以新坐标为(52,12)(\frac{5}{\sqrt{2}},-\frac{1}{\sqrt{2}})

我们也可以用矩阵来表示(x,y是行向量;a用列向量)
[xy][32]=[12,1212,12][32]=[5212]\begin{bmatrix} x \\ y\end{bmatrix}\begin{bmatrix} 3 \\ 2\end{bmatrix} = \begin{bmatrix} \frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}} \\ -\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}\end{bmatrix}\begin{bmatrix} 3 \\ 2\end{bmatrix} = \begin{bmatrix} \frac{5}{\sqrt{2}} \\ -\frac{1}{\sqrt{2}}\end{bmatrix}

协方差

协方差(Covariance)的定义:

cov(X,Y)=E[(XE[X])(YE[Y])]cov(X,Y)=1ni=1n(xiE(X))(yiE(Y)).{\displaystyle \operatorname {cov} (X,Y)=\operatorname {E} {{\big [}(X-\operatorname {E} [X])(Y-\operatorname {E} [Y]){\big ]}}} \\ {\displaystyle \operatorname {cov} (X,Y)={\frac {1}{n}}\sum _{i=1}^{n}(x_{i}-E(X))(y_{i}-E(Y)).}

性质:
cov(X,Y)=E[(XE[X])(YE[Y])]=E[XYXE[Y]E[X]Y+E[X]E[Y]]=E[XY]E[X]E[Y]E[X]E[Y]+E[X]E[Y]=E[XY]E[X]E[Y],{\displaystyle {\begin{aligned}\operatorname {cov} (X,Y)&=\operatorname {E} \left[\left(X-\operatorname {E} \left[X\right]\right)\left(Y-\operatorname {E} \left[Y\right]\right)\right]\\&=\operatorname {E} \left[XY-X\operatorname {E} \left[Y\right]-\operatorname {E} \left[X\right]Y+\operatorname {E} \left[X\right]\operatorname {E} \left[Y\right]\right]\\&=\operatorname {E} \left[XY\right]-\operatorname {E} \left[X\right]\operatorname {E} \left[Y\right]-\operatorname {E} \left[X\right]\operatorname {E} \left[Y\right]+\operatorname {E} \left[X\right]\operatorname {E} \left[Y\right]\\&=\operatorname {E} \left[XY\right]-\operatorname {E} \left[X\right]\operatorname {E} \left[Y\right],\end{aligned}}}

cov(X,a)=0cov(X,X)=var(X)cov(X,Y)=cov(Y,X)cov(aX,bY)=abcov(X,Y)cov(X+a,Y+b)=cov(X,Y)cov(aX+bY,cW+dV)=accov(X,W)+adcov(X,V)+bccov(Y,W)+bdcov(Y,V){\displaystyle {\begin{aligned}\operatorname {cov} (X,a)&=0\\\operatorname {cov} (X,X)&=\operatorname {var} (X)\\\operatorname {cov} (X,Y)&=\operatorname {cov} (Y,X)\\\operatorname {cov} (aX,bY)&=ab\,\operatorname {cov} (X,Y)\\\operatorname {cov} (X+a,Y+b)&=\operatorname {cov} (X,Y)\\\operatorname {cov} (aX+bY,cW+dV)&=ac\,\operatorname {cov} (X,W)+ad\,\operatorname {cov} (X,V)+bc\,\operatorname {cov} (Y,W)+bd\,\operatorname {cov} (Y,V)\end{aligned}}}

协方差矩阵

协方差矩阵(Covariance matrix)的定义:对称的方阵
XX是个随机向量(random vector),每个实体(随机变量)就是一个列向量,就是矩阵用列向量表示;
X=(X1,X2,...,Xn)T{\displaystyle \mathbf {X} =(X_{1},X_{2},...,X_{n})^{\mathrm {T} }}
XX的协方差矩阵用KXX{\displaystyle \operatorname {K} _{\mathbf {X} \mathbf {X} }}表示,矩阵中的每个元素KXiXj=cov[Xi,Xj]=E[(XiE[Xi])(XjE[Xj])]{\displaystyle \operatorname {K} _{X_{i}X_{j}}=\operatorname {cov} [X_{i},X_{j}]=\operatorname {E} [(X_{i}-\operatorname {E} [X_{i}])(X_{j}-\operatorname {E} [X_{j}])]}

KXX=[E[(X1E[X1])(X1E[X1])]E[(X1E[X1])(X2E[X2])]E[(X1E[X1])(XnE[Xn])]E[(X2E[X2])(X1E[X1])]E[(X2E[X2])(X2E[X2])]E[(X2E[X2])(XnE[Xn])]E[(XnE[Xn])(X1E[X1])]E[(XnE[Xn])(X2E[X2])]E[(XnE[Xn])(XnE[Xn])]]{\displaystyle \operatorname {K} _{\mathbf {X} \mathbf {X} }={\begin{bmatrix}\mathrm {E} [(X_{1}-\operatorname {E} [X_{1}])(X_{1}-\operatorname {E} [X_{1}])]&\mathrm {E} [(X_{1}-\operatorname {E} [X_{1}])(X_{2}-\operatorname {E} [X_{2}])]&\cdots &\mathrm {E} [(X_{1}-\operatorname {E} [X_{1}])(X_{n}-\operatorname {E} [X_{n}])]\\\\\mathrm {E} [(X_{2}-\operatorname {E} [X_{2}])(X_{1}-\operatorname {E} [X_{1}])]&\mathrm {E} [(X_{2}-\operatorname {E} [X_{2}])(X_{2}-\operatorname {E} [X_{2}])]&\cdots &\mathrm {E} [(X_{2}-\operatorname {E} [X_{2}])(X_{n}-\operatorname {E} [X_{n}])]\\\\\vdots &\vdots &\ddots &\vdots \\\\\mathrm {E} [(X_{n}-\operatorname {E} [X_{n}])(X_{1}-\operatorname {E} [X_{1}])]&\mathrm {E} [(X_{n}-\operatorname {E} [X_{n}])(X_{2}-\operatorname {E} [X_{2}])]&\cdots &\mathrm {E} [(X_{n}-\operatorname {E} [X_{n}])(X_{n}-\operatorname {E} [X_{n}])]\end{bmatrix}}}

样本的协方差(无偏)
cov(X,Y)=1n1i=1n(XiXˉ)(YiYˉ)cov(X,Y) = \frac{1}{n - 1}\sum_{i=1}^{n}\left ( X_{i} - \bar{X} \right )\left ( Y_{i} - \bar{Y} \right )

样本的协方差矩阵
cov[Xn×p]n×n=cov[X1,...,Xn]=1n1KXXcov[X_{n \times p}]_{n \times n} = cov[X_1,...,X_n] = \frac{1}{n-1}{K} _{\mathbf {X} \mathbf {X} }
这里很多协方差函数都有参数,可以设置到底是按行向量还是列向量计算协方差。
也有些地方是用1n\frac{1}{n},就是无偏与有偏的区别。

np.cov(X3×2,rowvar=False)\text{np.cov}(X_{3 \times 2},rowvar = False)输出2×22 \times 2(rowvar = False表示一列是一个特征);默认是输出3×33 \times 3(默认是行表示一个特征)
np.cov(x,y,z)\text{np.cov}(x,y,z)输出3×33 \times 3

期望

期望(Expectation)的定义:
E[X]=i=1kxipi=x1p1+x2p2++xkpk.p1+p2++pk=1,{\displaystyle \operatorname {E} [X]=\sum _{i=1}^{k}x_{i}\,p_{i}=x_{1}p_{1}+x_{2}p_{2}+\cdots +x_{k}p_{k}.} \\ {\displaystyle p_{1}+p_{2}+\cdots +p_{k}=1,}

性质:
E[X+Y]=E[X]+E[Y],E[aX]=aE[X],{\displaystyle {\begin{aligned}\operatorname {E} [X+Y]&=\operatorname {E} [X]+\operatorname {E} [Y],\\\operatorname {E} [aX]&=a\operatorname {E} [X],\end{aligned}}}

如果X,YX,Y是相互独立的,那么E[XY]=E[X]E[Y]{\displaystyle \operatorname {E} [XY]=\operatorname {E} [X]\operatorname {E} [Y]}

常数的期望等于常数本身E(a)=aE(a) = a

xix_i是随机变量XX的一个实例,XX服从什么分布,xix_i也是服从什么分布的,所以E(xi)=E(X),D(xi)=D(X)E(x_i) = E(X),D(x_i) = D(X)

E(X)E(X)为一阶矩
E(X2)E(X^2)为二阶矩

原点矩(raw moment)和中心矩(central moment
E(Xk)E(X^k)为k阶远点矩,一阶原点矩是数学期望
E(XE(X))kE(X-E(X))^k为k阶中心矩,二阶原点矩是方差(以E(X)E(X)为中心)

方差

方差(Variance)的定义:
Var(X)=E[(Xμ)2]μ=E[X]\operatorname {Var} (X)=\operatorname {E} \left[(X-\mu )^{2}\right] \\ {\displaystyle \mu =\operatorname {E} [X]}

性质:
D(X)D(X)Var(X)Var(X) 都是表示方差;方差大于等于0;参数的方差为0;

Var(X)=Cov(X,X).\operatorname {Var} (X)=\operatorname {Cov} (X,X).
Var(X)=E[(XE[X])2]=E[X22XE[X]+E[X]2]=E[X2]2E[X]E[X]+E[X]2=E[X2]E[X]2{\displaystyle {\begin{aligned}\operatorname {Var} (X)&=\operatorname {E} \left[(X-\operatorname {E} [X])^{2}\right]\\[4pt]&=\operatorname {E} \left[X^{2}-2X\operatorname {E} [X]+\operatorname {E} [X]^{2}\right]\\[4pt]&=\operatorname {E} \left[X^{2}\right]-2\operatorname {E} [X]\operatorname {E} [X]+\operatorname {E} [X]^{2}\\[4pt]&=\operatorname {E} \left[X^{2}\right]-\operatorname {E} [X]^{2}\end{aligned}}}

Var(X+Y)=E[(X+Y)2](E[X+Y])2=E[X2+2XY+Y2](E[X]+E[Y])2=E[X2]+2E[XY]+E[Y2](E[X]2+2E[X]E[Y]+E[Y]2)=E[X2]+E[Y2]E[X]2E[Y]2=Var(X)+Var(Y){\displaystyle {\begin{aligned}\operatorname {Var} (X+Y)&=\operatorname {E} \left[(X+Y)^{2}\right]-(\operatorname {E} [X+Y])^{2}\\[5pt]&=\operatorname {E} \left[X^{2}+2XY+Y^{2}\right]-(\operatorname {E} [X]+\operatorname {E} [Y])^{2} \\&=\operatorname {E} \left[X^{2}\right]+2\operatorname {E} [XY]+\operatorname {E} \left[Y^{2}\right]-\left(\operatorname {E} [X]^{2}+2\operatorname {E} [X]\operatorname {E} [Y]+\operatorname {E} [Y]^{2}\right)\\[5pt]&=\operatorname {E} \left[X^{2}\right]+\operatorname {E} \left[Y^{2}\right]-\operatorname {E} [X]^{2}-\operatorname {E} [Y]^{2}\\[5pt]&=\operatorname {Var} (X)+\operatorname {Var} (Y)\end{aligned}}}

Var(X+a)=Var(X).\operatorname {Var} (X+a)=\operatorname {Var} (X).
Var(aX)=a2Var(X).\operatorname {Var} (aX)=a^{2}\operatorname {Var} (X).
Var(aX+bY)=a2Var(X)+b2Var(Y)+2abCov(X,Y),\operatorname {Var} (aX+bY)=a^{2}\operatorname {Var} (X)+b^{2}\operatorname {Var} (Y)+2ab\,\operatorname {Cov} (X,Y),
Var(aXbY)=a2Var(X)+b2Var(Y)2abCov(X,Y),\operatorname {Var} (aX-bY)=a^{2}\operatorname {Var} (X)+b^{2}\operatorname {Var} (Y)-2ab\,\operatorname {Cov} (X,Y),

Var(XY)=E(X2)E(Y2)[E(X)]2[E(Y)]2.{\displaystyle \operatorname {Var} (XY)=\operatorname {E} \left(X^{2}\right)\operatorname {E} \left(Y^{2}\right)-[\operatorname {E} (X)]^{2}[\operatorname {E} (Y)]^{2}.}

有偏估计和无偏估计

参数估计需要未知参数的估计量和一定置信度
估计方法:用点估计估计一个值;用区间估计估计值的可能区间和是该值的可能性。
估计的偏差的定义:
bias(θ^m)=E(θ^m)θbias(\hat{\theta}_m) = E(\hat{\theta}_m)-\theta
θ\theta是数据分布真实参数(完美)
θ\theta的估计量或统计量θ^m\hat{\theta}_m

对估计值的评价标准:

无偏估计
例如样本均值的估计为μ^=i=1mxi\hat{\mu} = \sum_{i=1}^m x_i,真实的均值为μ\mu,如何知道这个估计是有偏还是无偏?根据定义判断偏差是否为0;bias(μ^)=E(μ^)μbias(\hat{\mu}) = E(\hat{\mu}) - \mu
这里的证明再有效性已经证明过了。

有偏估计
如果样本方差的估计为σ^2=1mi=1m(xiμ^)2\hat{\sigma}^2 = \frac{1}{m}\sum_{i=1}^m(x_i-\hat{\mu})^2,μ^=i=1mxi\hat{\mu} = \sum_{i=1}^m x_i,那么偏差bias(σ^2)=E[σ^2]σ2bias(\hat{\sigma}^2) = E[\hat{\sigma}^2] - \sigma^2不为0,就证明这个估计是有偏的。
我们来求E[σ^2]E[\hat{\sigma}^2]
E[σ^2]=E[1mi=1m(xiμ^)2]=E[1mi=1m(xi22μ^xi+μ^2)]=E[1mi=1mxi21mi=1m2μ^xi+1mi=1mμ^2]=E[1mi=1mxi22μ^1mi=1mxi+1mi=1mμ^2]=E[1mi=1mxi22μ^μ^+μ^2]=E[1mi=1mxi2μ^2]=E[1mi=1mxi2]E[μ^2]=1mE[i=1mxi2]E[μ^2]E[\hat{\sigma}^2] = E[\frac{1}{m}\sum_{i=1}^m(x_i-\hat{\mu})^2] \\= E[\frac{1}{m}\sum_{i=1}^m(x_i^2- 2\hat{\mu}x_i +\hat{\mu}^2)] \\= E[\frac{1}{m}\sum_{i=1}^mx_i^2-\frac{1}{m}\sum_{i=1}^m 2\hat{\mu}x_i +\frac{1}{m}\sum_{i=1}^m\hat{\mu}^2]\\= E[\frac{1}{m}\sum_{i=1}^mx_i^2-2\hat{\mu}\frac{1}{m}\sum_{i=1}^m x_i +\frac{1}{m}\sum_{i=1}^m\hat{\mu}^2]\\= E[\frac{1}{m}\sum_{i=1}^mx_i^2-2\hat{\mu}\hat{\mu} +\hat{\mu}^2]= E[\frac{1}{m}\sum_{i=1}^mx_i^2-\hat{\mu}^2] \\= E[\frac{1}{m}\sum_{i=1}^mx_i^2]-E[\hat{\mu}^2] = \frac{1}{m}E[\sum_{i=1}^mx_i^2]-E[\hat{\mu}^2]

E[μ^2]=D(μ^)+E(μ^)2E[\hat{\mu}^2] = D(\hat{\mu}) + E(\hat{\mu})^2
D(μ^)=1mσ2E(μ^)=μD(\hat{\mu}) = \frac{1}{m} \sigma^2 和 E(\hat{\mu}) = \mu在上面的有效性中已经证明了
所以E[μ^2]=1mσ2+μ2E[\hat{\mu}^2] = \frac{1}{m} \sigma^2 + \mu^2

E[xi2]=D(xi)+E(xi)2E[{x_i}^2] = D(x_i) + E(x_i)^2
D(xi)=σ2E(xi)=μD(x_i) = \sigma^2 和 E(x_i) = \mu在上面的有效性中已经证明了
所以1mE[i=1mxi2]=1mi=1mE[xi2]=σ2+μ2\frac{1}{m}E[\sum_{i=1}^mx_i^2] = \frac{1}{m}\sum_{i=1}^m E[x_i^2]= \sigma^2 + \mu^2

所以E[σ^2]=σ2+μ2(1mσ2+μ2)=m1mσ2E[\hat{\sigma}^2] = \sigma^2 + \mu^2 - (\frac{1}{m} \sigma^2 + \mu^2) =\frac{m-1}{m}\sigma^2
所以估计是有偏估计。

所以方差的无偏估计为σ~2=1m1i=1m(xiμ^)2\tilde{\sigma}^2 = \frac{1}{m-1}\sum_{i=1}^m(x_i-\hat{\mu})^2

limm1mi=1m(xiμ^)2=1m1i=1m(xiμ^)2\lim_{m \to \infty}\frac{1}{m}\sum_{i=1}^m(x_i-\hat{\mu})^2 = \frac{1}{m-1}\sum_{i=1}^m(x_i-\hat{\mu})^2,也就是说有偏估计是一个渐近无偏估计。

无偏估计不一定是最好的估计!

因子分析FA

因子分析(Factor analysis, FA

每个变量都可以表示成公共因子的线性函数与特殊因子之和
Xi=ai1F1+ai2F2+...++aimFm+ϵi,(i=1,2,...,p)X_i = a_{i1}F_1 + a_{i2}F_2 +...++ a_{im}F_m + \epsilon_i ,(i=1,2,...,p)
式中的F1,F2,…,Fm称为公共因子,εi称为Xi的特殊因子。该模型可用矩阵表示为:X = AF+ε
X 表示原始数据,矩阵A称为因子载荷矩阵,F表示公共因子, ε是特殊因子
aij称为因子“载荷”,是第i个变量在第j个因子上的负荷,如果把变量Xi看成m维空间中的一个点,则aij表示它在坐标轴Fj上的投影。
X=[X1,X2...Xp]TX = [X_1,X_2...X_p]^T
A=[a11a12...a1ma21a22...a2m............ap1ap2...apm]A= \begin{bmatrix} a_{11} & a_{12} & ... & a_{1m} \\\\ a_{21} & a_{22} & ... & a_{2m} \\\\ ... & ... & ... & ... \\\\ a_{p1} & a_{p2} & ... & a_{pm} \\\\ \end{bmatrix}
F=[F1,F2...Fm]TF = [F_1,F_2...F_m]^T
ϵ=[ϵ1,ϵ2...ϵp]T\epsilon = [\epsilon_1,\epsilon_2...\epsilon_p]^T

主成分分析,是分析维度属性的主要成分表示。
因子分析,是分析属性们的公共部分的表示。
二者均应用于高斯分布的数据,非高斯分布的数据采用独立成分分析ICA算法

独立成分分析ICA

独立成分分析(Independent component analysis, ICA

X=AS
Y=WX=WAS , A = inv(W)
ICA(Independent Component Correlation Algorithm)是一种函数,X为n维观测信号矢量,S为独立的m(m<=n)维未知源信号矢量,矩阵A被称为混合矩阵。
ICA的目的就是寻找解混矩阵W(A的逆矩阵),然后对X进行线性变换,得到输出向量Y。

ICA是找出构成信号的相互独立部分(不需要正交),对应高阶统计量分析。ICA理论认为用来观测的混合数据阵X是由独立元S经过A线性加权获得。
ICA理论的目标就是通过X求得一个分离矩阵W,使得W作用在X上所获得的信号Y是独立源S的最优逼近,

独立成分分析 (ICA) 应用参考(Origin来做ICA分析)

独立成分分析 - 讲解的原理

Independent Component Analysis (ICA)

参考文献

[16-1] 方开泰. 实用多元统计分析. 上海:华东师范大学出版社, 1989.
[16-2] 夏绍玮,杨家本,杨振斌. 系统工程概论. 北京:清华大学出版社,1995.
[16-3] Jolliffe I. Principal component analysis, Sencond Edition. John Wiley & Sons, 2002.
[16-4] Shlens J. A tutorial on principal component analysis. arXiv preprint arXiv: 14016.1100, 2014.
[16-5] Schölkopf B, Smola A, Müller K-R. Kernel principal component analysis. Artificial Neural Networks-ICANN'97. Springer, 1997:583-588.
[16-6] Hardoon D R, Szedmak S, Shawe-Taylor J. Canonical correlation analysis: an overview with application to learning methods. Neural Computation, 2004, 16(12):2639-2664.
[16-7] Candes E J, Li X D, Ma Y, et al. Robust Principal component analysis? Journal of the ACM(JACM), 2011, 58(3):11.

第 17 章 潜在语义分析

我们先介绍下文本信息处理中的一些问题:

  1. 一词多义(多义现象)polysemy
    分类时:比如bank 这个单词如果和mortgage, loans, rates 这些单词同时出现时,bank 很可能表示金融机构的意思。可是如果bank 这个单词和lures, casting, fish一起出现,那么很可能表示河岸的意思。

  2. 一义多词(同义现象)synonymy
    检索时:比如用户搜索“automobile”,即汽车,传统向量空间模型仅仅会返回包含“automobile”单词的页面,而实际上包含“car”单词的页面也可能是用户所需要的。

LSA能解决同义(语义相似度)问题:发现单词与主题之间的关系,这里主题是汽车;也能解决一定程度的多义问题,同一个单词在不同文档中表示不同话题

潜在语义分析(Latent semantic analysis, LSA)旨在 解决这种方法不能准确表示语义的问题,试图从大量的文本数据中发现潜在的话题,以话题向量表示文本的语义内容,以话题向量空间的度量更准确地表示文本之间的语义相似度。

文本doc集合D={d1,d2,...,dn}D = \{d_1,d_2,...,d_n\}
文本集合中出现的单词word集合W={w1,w2,...,wm}W = \{w_1,w_2,...,w_m\}
单词-文本矩阵(word-document matrix)
X=[x11x12x1nx21x22x2nxm1xm2xmn]X = \begin{bmatrix} x_{11} & x_{12} & \cdots & x_{1n} \\ x_{21} & x_{22} & \cdots & x_{2n} \\ \vdots & \vdots & & \vdots \\ x_{m1} & x_{m2} & \cdots & x_{mn} \\ \end{bmatrix}
每一列表示一个文本;xijx_{ij}表示单词wiw_i在文本djd_j中出现的频数或权值。
每个文本中不可能出现所有单词,所以该矩阵是稀疏矩阵。

权值通常用单词频率-逆文档频率term frequency–inverse document frequency,TFIDF)表示,定义为:
tf-idf(t,d,D)=tf(t,d)idf(t,D){\displaystyle \text {tf-idf} (t,d,D)=\mathrm {tf} (t,d)\cdot \mathrm {idf} (t,D)}
t为某一个单词(term,word);d为某一个文档(document);xij=tf-idf(wi,dj,D)x_{ij} = \text {tf-idf} (w_i,d_j,D);D表示文档集合,N=DN = |D|表示文档总数;
tf(t,d)=ft,dtdft,d=td中出现的频数d中出现的所有单词的频数和{\displaystyle \mathrm {tf} (t,d)={\frac {f_{t,d}}{\sum _{t'\in d}{f_{t',d}}}} = \frac{t在d中出现的频数}{d中出现的所有单词的频数和}}
idf(t,D)=logN{dD:td}=log文档总数含有单词t的文本总数\mathrm{idf}(t, D) = \log \frac{N}{|\{d \in D: t \in d\}| } = \log \frac{文档总数}{含有单词t的文本总数}

直观上理解:
一个单词在一个文本中出现的频数越高,这个单词在这个文本中的重要度(TF)就越高;
一个单词在整个文档集合中出现的文档数越少,这个单词就越能表示其所在文档的特点,重要度(TDF)就越高;
两种重要度的积,表示综合重要度。
意思就是重要的单词在一个文本中出现的越多越重要,在越少的文本中出现越重要;如:的,可能在每个文档中出现都很多(TF大),并且每个文档都有出现过(TDF小),所以反而不重要了。

相似度(余弦相似度)Cosine similarity)可以表示两个文本之间的语义相似度,越大越相似。
我们知道训练点积AB=ABcosθ{\displaystyle \mathbf {A} \cdot \mathbf {B} =\left\|\mathbf {A} \right\|\left\|\mathbf {B} \right\|\cos \theta },而相似度就是向量之间的夹角similarity=cos(θ)=ABAB=i=1nAiBii=1nAi2i=1nBi2,{\displaystyle {\text{similarity}}=\cos(\theta )={\mathbf {A} \cdot \mathbf {B} \over \|\mathbf {A} \|\|\mathbf {B} \|}={\frac {\sum \limits _{i=1}^{n}{A_{i}B_{i}}}{{\sqrt {\sum \limits _{i=1}^{n}{A_{i}^{2}}}}{\sqrt {\sum \limits _{i=1}^{n}{B_{i}^{2}}}}}},}

文档用向量表示:di=x.i=[x1ix2ixmi]d_i = x_{.i} = \begin{bmatrix} x_{1i} \\ x_{2i} \\ \vdots \\ x_{mi} \end{bmatrix}
那么di,djd_i,d_j之间的相似度为
similarity=x.ix.jx.ix.j\text{similarity} = \frac{x_{.i}\cdot x_{.j}}{\|x_{.i}\|\|x_{.j}\|}
这里比较相似度用在了单词向量空间(word vector space model)中,有一个问题就是多义和同义现象,这时我们就可以考虑话题向量空间(topic vector space model):
假设用一个向量表示文档,该向量的每一个分量表示一个话题,其数值为该话题在该文本中的权值,然后比较两个文档相似度(一般话题数远小于单词数)。
潜在语义分析就是构建这样一个话题向量空间的方法。

单词-文本矩阵
X=[x11x12x1nx21x22x2nxm1xm2xmn]X = \begin{bmatrix} x_{11} & x_{12} & \cdots & x_{1n} \\ x_{21} & x_{22} & \cdots & x_{2n} \\ \vdots & \vdots & & \vdots \\ x_{m1} & x_{m2} & \cdots & x_{mn} \\ \end{bmatrix}
假设所有文档共含有k个话题,话题向量空间TT(单词-话题矩阵)
T=[t11t12t1kt21t22t2ktm1tm2tmk]T = \begin{bmatrix} t_{11} & t_{12} & \cdots & t_{1k} \\ t_{21} & t_{22} & \cdots & t_{2k} \\ \vdots & \vdots & & \vdots \\ t_{m1} & t_{m2} & \cdots & t_{mk} \\ \end{bmatrix}

那么文档在话题向量空间的表示(话题-文本矩阵)
Y=[y11y12y1ny21y22y2nyk1yk2ykn]Y = \begin{bmatrix} y_{11} & y_{12} & \cdots & y_{1n} \\ y_{21} & y_{22} & \cdots & y_{2n} \\ \vdots & \vdots & & \vdots \\ y_{k1} & y_{k2} & \cdots & y_{kn} \\ \end{bmatrix}

参考文献

[17-1] Deerwester S C, Dumais S T, Landauer T K, et al. Indexing by latent semantic analysis. Journal of the Association for Information Science and Technology ,1990, 41: 391-407.
[17-2] Landauer T K. Latent semantic analysis. In: Encyclopedia of Cognitive Science, Wiley, 2006.
[17-3] Lee D D, Seung H S. Learning the parts of objects by non-negative matrix factorization. Nature, 1999, 401(6755):788-791.
[17-4] Lee D D, Seung H S. Algorithms for non-negative matrix factorization. Advances in Neural Information Processing Systems, 2001: 556-562.
[17-5] Xu W, Liu X, Gong Y. Document clustering based on non-negative matrix factorization. Proceedings of the 26th Annual International ACM SIGIR Conference in Research and Development in Information Retrieval, 2003.
[17-6] Wang Q, Xu J, Li H, et al. Regularized latent semantic indexing. Proceedings of the 34th International ACM SIGIR Conference in Research and Development in Information Retrieval, 2011.

第 18 章 概率潜在语义分析

生成模型,用隐变量表示话题

概率潜在语义分析(Probabilistic latent semantic analysis, PLSA

概率有向图模型:

阴影圆表示观测变量,空心圆表示隐变量;箭头表示概率关系;方框表示多次重复,方框内的字母表示重复次数;
文档d是一个观测变量;话题变量z是隐变量(话题的个数是超参数);单词变量w是一个观测变量;

  1. 依据概率分布P(d),从文本集合中随机选取一个文本d,共生成N个文本;针对每个文本,执行下一步操作
  2. 在文本d给定条件下,依据条件概率分布P(z|d),从话题集合中随机选取一个话题z,共生成L个话题,这里L是文本长度(每个话题生成一个单词,所以生成的文本d的长度是L)
  3. 在话题z给定条件下,依据概率分布P(w|z),从单词集合中随机选取一个单词w

文本-单词共现数据(矩阵T)的生成概率为所有单词-文本对(w,d)的生成概率乘积,
P(T)=(w,d)P(w,d)n(w,d)T=[n(wi,dj)]i=1,2,...,M;j=1,2,...,NP(T) = \prod_{(w,d)} P(w,d)^{n(w,d)} \\ T = [n(w_i,d_j)] i=1,2,...,M; j=1,2,...,N
n(w,d)n(w,d)表示(w,d)出现的次数;
矩阵T的行表示单词,列表示文本,元素表示(w,d)出现的次数,用n(w,d)n(w,d)表示;
单词-文本对 出现的总次数为N×L

这里假设文本的长度都是等长的,正常情况是第一个文本的长度是L1,...,第N个文本的长度是LN;
正常情况下单词-文本对 出现的总次数为i=1NLi\sum_{i=1}^N L_i

书中还讲到了等价的共现模型(对称模型)
P(w,d)=zP(z)P(dz)P(wz)P(w,d)=\sum _{z}P(z)P(d|z)P(w|z)

  1. 现在进行E步,计算Q函数
    argmaxθL(θ)=argmaxθi=1Mj=1Nn(wi,dj)logk=1KP(zkwi,dj)P(wizk)P(zkdj)P(zkwi,dj)\arg \max _{\theta} L(\theta)=\arg \max _{\theta} \sum_{i=1}^{M} \sum_{j=1}^{N} n\left(w_{i}, d_{j}\right) \log \sum_{k=1}^{K} P\left(z_{k} \mid w_{i}, d_{j}\right) \frac{P\left(w_{i} \mid z_{k}\right) P\left(z_{k} \mid d_{j}\right)}{P\left(z_{k} \mid w_{i}, d_{j}\right)}
    其中log右边为关于z的期望:
    Ez[P(wizk)P(zkdj)P(zkwi,dj)]=k=1KP(zkwi,dj)P(wizk)P(zkdj)P(zkwi,dj)E_z \left[\frac{P\left(w_{i} \mid z_{k}\right) P\left(z_{k} \mid d_{j}\right)}{P\left(z_{k} \mid w_{i}, d_{j}\right)}\right] = \sum_{k=1}^{K} P\left(z_{k} \mid w_{i}, d_{j}\right) \frac{P\left(w_{i} \mid z_{k}\right) P\left(z_{k} \mid d_{j}\right)}{P\left(z_{k} \mid w_{i}, d_{j}\right)}
    所以:
    argmaxθL(θ)=argmaxθi=1Mj=1Nn(wi,dj)logEz[P(wizk)P(zkdj)P(zkwi,dj)]\arg \max _{\theta} L(\theta)=\arg \max _{\theta} \sum_{i=1}^{M} \sum_{j=1}^{N} n\left(w_{i}, d_{j}\right) \log E_z \left[\frac{P\left(w_{i} \mid z_{k}\right) P\left(z_{k} \mid d_{j}\right)}{P\left(z_{k} \mid w_{i}, d_{j}\right)}\right]
    根据jensen不等式,我们可以得到:
    logEz[P(wizk)P(zkdj)P(zkwi,dj)]Ez[logP(wizk)P(zkdj)P(zkwi,dj)]\log E_{z}\left[\frac{P\left(w_{i} \mid z_{k}\right) P\left(z_{k} \mid d_{j}\right)}{P\left(z_{k} \mid w_{i}, d_{j}\right)}\right] \geq E_{z}\left[\log \frac{P\left(w_{i} \mid z_{k}\right) P\left(z_{k} \mid d_{j}\right)}{P\left(z_{k} \mid w_{i}, d_{j}\right)}\right]
    得到L(θ)的下界:
    i=1Mj=1Nn(wi,dj)Ez[logP(wizk)P(zkdj)P(zkwi,dj)]=i=1Mj=1Nn(wi,dj)k=1KP(zkwi,dj)[logP(wizk)P(zkdj)logP(zkwi,dj)]\sum_{i=1}^{M} \sum_{j=1}^{N} n\left(w_{i}, d_{j}\right) E_z \left[\log \frac{P\left(w_{i} \mid z_{k}\right) P\left(z_{k} \mid d_{j}\right)}{P\left(z_{k} \mid w_{i}, d_{j}\right)}\right] \\ =\sum_{i=1}^{M} \sum_{j=1}^{N} n\left(w_{i}, d_{j}\right) \sum_{k=1}^{K} P\left(z_{k} \mid w_{i}, d_{j}\right)\left[\log P\left(w_{i} \mid z_{k}\right) P\left(z_{k} \mid d_{j}\right)-\log P\left(z_{k} \mid w_{i}, d_{j}\right)\right]
    最后,我们将K的累加项拆开,可以得到一项 k=1KP(zkwi,dj)logP(zkwi,dj)\sum_{k=1}^{K} P\left(z_{k} \mid w_{i}, d_{j}\right) \log P\left(z_{k} \mid w_{i}, d_{j}\right) ,这一项在M步中没有作用,可以省去,于是我们可以得到Q函数为:
    Q=i=1Mj=1Nn(wi,dj)k=1KP(zkwi,dj)log[P(wizk)P(zkdj)]Q=\sum_{i=1}^{M} \sum_{j=1}^{N} n\left(w_{i}, d_{j}\right) \sum_{k=1}^{K} P\left(z_{k} \mid w_{i}, d_{j}\right) \log \left[P\left(w_{i} \mid z_{k}\right) P\left(z_{k} \mid d_{j}\right)\right]
    需要优化的参数为 P(zkwi,dj)P(wizk)P(zkdj)P\left(z_{k} \mid w_{i}, d_{j}\right),P\left(w_{i} \mid z_{k}\right), P\left(z_{k} \mid d_{j}\right) 这三项,在Q步中,第一项是变量,后两项是常量,于是可以由贝叶斯公式获得:
    P(zkwi,dj)=P(wizk)P(zkdj)k=1KP(wizk)P(zkdj)P\left(z_{k} \mid w_{i}, d_{j}\right)=\frac{P\left(w_{i} \mid z_{k}\right) P\left(z_{k} \mid d_{j}\right)}{\sum_{k=1}^{K} P\left(w_{i} \mid z_{k}\right) P\left(z_{k} \mid d_{j}\right)}

  2. M步
    在M步中,我们需要优化的是P(wizk)P(zkdj)P\left(w_{i} \mid z_{k}\right), P\left(z_{k} \mid d_{j}\right) 这两项(两项的乘积代表的完全数据,是未知变量),此时 P(zkwi,dj)P\left(z_{k} \mid w_{i}, d_{j}\right)为常量(代表不完全数据,是已知变量),极大化Q函数的M步可以使用拉格朗日乘子法来优化两个参数,即:
    maxQs.t.i=1MP(wizk)=1,k=1,2,,Kk=1KP(zkdj)=1,j=1,2,,N\max Q \\s.t. \quad \begin{array}{l}\sum_{i=1}^{M} P\left(w_{i} \mid z_{k}\right)=1, \quad k=1,2, \cdots, K \\ \sum_{k=1}^{K} P\left(z_{k} \mid d_{j}\right)=1, \quad j=1,2, \cdots, N\end{array}
    根据上述约束条件构造拉格朗日函数:
    Λ=Q+k=1Kτk(1i=1MP(wizk))+j=1Nρj(1k=1KP(zkdj))\Lambda=Q^{\prime}+\sum_{k=1}^{K} \tau_{k}\left(1-\sum_{i=1}^{M} P\left(w_{i} \mid z_{k}\right)\right)+\sum_{j=1}^{N} \rho_{j}\left(1-\sum_{k=1}^{K} P\left(z_{k} \mid d_{j}\right)\right)
    分别对两个参数P(wizk)P(zkdj)P\left(w_{i} \mid z_{k}\right), P\left(z_{k} \mid d_{j}\right)求偏导,并令偏导数为0:
    j=1Nn(wi,dj)P(zkwi,dj)τkP(wizk)=0,i=1,2,,M;k=1,2,,Ki=1Mn(wi,dj)P(zkwi,dj)ρjP(zkdj)=0,j=1,2,,N;k=1,2,,K\begin{array}{l}\sum_{j=1}^{N} n\left(w_{i}, d_{j}\right) P\left(z_{k} \mid w_{i}, d_{j}\right)-\tau_{k} P\left(w_{i} \mid z_{k}\right)=0, \quad i=1,2, \cdots, M ; \quad k=1,2, \cdots, K \\ \sum_{i=1}^{M} n\left(w_{i}, d_{j}\right) P\left(z_{k} \mid w_{i}, d_{j}\right)-\rho_{j} P\left(z_{k} \mid d_{j}\right)=0, \quad j=1,2, \cdots, N ; \quad k=1,2, \cdots, K\end{array}
    求解上面的方程组,就可以得到M步的参数估计:
    P(wizk)=j=1Nn(wi,dj)P(zkwi,dj)m=1Mj=1Nn(wm,dj)P(zkwm,dj)P\left(w_{i} \mid z_{k}\right)=\frac{\sum_{j=1}^{N} n\left(w_{i}, d_{j}\right) P\left(z_{k} \mid w_{i}, d_{j}\right)}{\sum_{m=1}^{M} \sum_{j=1}^{N} n\left(w_{m}, d_{j}\right) P\left(z_{k} \mid w_{m}, d_{j}\right)}
    P(zkdj)=i=1Mn(wi,dj)P(zkwi,dj)n(dj)P\left(z_{k} \mid d_{j}\right)=\frac{\sum_{i=1}^{M} n\left(w_{i}, d_{j}\right) P\left(z_{k} \mid w_{i}, d_{j}\right)}{n\left(d_{j}\right)}
    最后,在E步和M步间不停迭代,直到得到优化后的两个参数
    n(dj)=i=1Mn(wi,dj)n(d_j) = \sum_{i=1}^M n(w_i,d_j)表示文本djd_j中的单词个数,n(wi,dj)n(w_i,d_j)表示单词wiw_i在文本djd_j中出现的次数。

参考文献

[18-1] Hofmann T. Probabilistic Latent Semantic analysis. Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence, 1999: 289-296.
[18-2] Hofmann T. Probabilistic Latent Semantic Indexing. Proceedings of the 22nd Annual International ACM SIGIR Conference in Research and Development in Information Retrieval, 1999.
[18-3] Hofmann T. Unsupervised learning by probabilistic latent semantic analysis. Machine Learning, 2001, 42: 177-196.
[18-4] Ding C, Li T, Peng W. On the equivalence between non-negative matrix factorization and probabilistic latent semantic indexing. Computational Statistics & Data Analysis, 2008, 52(8): 3913-3927.

第 19 章 马尔可夫链蒙特卡罗法

马尔可夫链蒙特卡罗Markov Chain Monte Carlo, MCMC)由两个MC组成,即蒙特卡罗方法Monte Carlo Simulation, MC)和马尔可夫链Markov Chain, MC)。

要弄懂MCMC的原理我们首先得搞清楚蒙特卡罗方法和马尔可夫链的原理。

马尔可夫链在前面的章节有讲到,再结合书中的内容。这里补充下几个知识:
马尔可夫链的遍历定理
书中啰嗦了很多,我的理解是遍历足够多能达到平稳分布的马尔可夫链。并且达到任何一个状态的概率不能为0;(不可约,非周期且正常返)

可逆马尔可夫链(reversible Markov chain):
设有马尔可夫链X={X0,X1,...,Xt,...}X=\{X_0,X_1,...,X_t,...\},状态空间S,转移矩阵P,如果有状态分布π=(π1,π2,...)T\pi = (\pi_1,\pi_2,...)^T,对任意状态i,jSi,j \in S,对任意时刻t满足
P(Xt=iXt1=j)πj=P(Xt1=jXt=i)πi,i,j=1,2,...P(X_t=i|X_{t-1}=j)\pi_j = P(X_{t-1}=j|X_t=i)\pi_i ,\quad i,j =1,2,...
或者简写
pijπj=pjiπi,i,j=1,2,...p_{ij}\pi_j = p_{ji}\pi_i ,\quad i,j =1,2,...
则称此马尔可夫链为可逆马尔可夫链;简写的等式称为细致平衡方程(detailed balance equation),并且满足细致平衡方程的状态分布π\pi就是该马尔可夫链的平稳分布(并不是所有的马尔可夫链都是可逆的)。
可逆马尔可夫链满足遍历定理。

采样法Sampling Method)也称为蒙特卡罗方法(Monte Carlo Method, MC)或统计模拟方法(Statistical Simulation Method)

蒙特卡罗方法诞生于20 世纪 40 年代美国的“曼哈顿计划”,其名字来源于摩纳哥的一个以赌博业闻名的城市蒙特卡罗,象征概率.
蒙特卡罗方法是一种通过随机采样来近似估计一些计算问题数值解(Numerical solution与其对应的是闭式解Closed-form solution或解析解Analytical solution)的方法.

最早的蒙特卡罗方法都是为了求解一些不太好求解的求和或者积分问题。比如积分:θ=abf(x)dx\theta = \int_a^b f(x)dx 或者估计π\pi值或圆的面积(积分)。

我们可以通过蒙特卡罗方法来模拟求解近似值。如何模拟呢?

随机采样指从给定概率密度函数 𝑝(𝑥) 中抽取出符合其概率分布的样本.

随机采样 采样法的难点是如何进行随机采样,即如何让计算机生成满足概率密度函数 𝑝(𝑥) 的样本.我们知道,计算机可以比较容易地随机生成一个在 [0, 1]区间上均布分布的样本 𝜉.如果要随机生成服从某个非均匀分布的样本,就需要一些间接的采样方法.
如果一个分布的概率密度函数为 𝑝(𝑥),其累积分布函数 cdf(𝑥) 为连续的严格增函数,且存在逆函数cdf1(𝑦),𝑦[0,1]cdf^{−1}(𝑦), 𝑦 ∈ [0, 1],那么我们可以利用累积分布函数的逆函数(inverse CDF)来生成服从该随机分布的样本.假设 𝜉 是 [0, 1] 区间上均匀分布的随机变量,则 cdf1(ξ)cdf^{−1}(\xi) 服从概率密度函数为𝑝(𝑥)的分布.但当 𝑝(𝑥) 非常复杂,其累积分布函数的逆函数难以计算,或者不知道 𝑝(𝑥)的精确值,只知道未归一化的分布 ̂𝑝(𝑥)时,就难以直接对𝑝(𝑥)进行采样,往往需要使用一些间接的采样策略,比如拒绝采样、重要性采样、马尔可夫链蒙特卡罗采样等.这些方法一般是先根据一个比较容易采样的分布进行采样,然后通过一些策略来间接得到符合𝑝(𝑥)分布的样本.

rejection sampling, inverse CDF, Box-Muller, Ziggurat algorithm

拒绝采样(Rejection Sampling)是一种间接采样方法,也称为接受-拒绝采样(Acceptance-Rejection Sampling).
假设原始分布𝑝(𝑥)难以直接采样,我们可以引入一个容易采样的分布𝑞(𝑥),一般称为提议分布(Proposal Distribution),然后以某个标准来拒绝一部分的样本使得最终采集的样本服从分布 𝑝(𝑥)。我们需要构建一个提议分布 𝑞(𝑥) 和一个常数 𝑘,使得 𝑘𝑞(𝑥) 可以覆盖函数𝑝(𝑥),即𝑘𝑞(𝑥) ≥ 𝑝(𝑥),

对于每次抽取的样本 xˆ\^{x} 计算接受概率(Acceptance Probability):α(xˆ)=p(xˆ)kq(xˆ)\alpha(\^{x}) = \frac{p(\^{x})}{kq(\^{x})},并以概率α(xˆ)\alpha(\^{x}))来接受样本 xˆ\^{x}

拒绝采样的采样过程

输入: 提议分布𝑞(𝑥),常数𝑘,样本集合𝒱 = ∅;
1 repeat
2   根据𝑞(𝑥)随机生成一个样本 ;̂𝑥
3   计算接受概率𝛼( ̂𝑥);
4   从(0, 1)的均匀分布中随机生成一个值𝑧;
5   if 𝑧 ≤ 𝛼( ̂𝑥) then // 以𝛼(𝑥)̂ 的概率接受𝒙̂
6       𝒱 = 𝒱 ∪ { ̂𝑥};
7   end
8 until 获得 𝑁 个样本(|𝒱| = 𝑁);
输出: 样本集合𝒱

判断一个拒绝采样方法的好坏就是看其采样效率,即总体的接受率.如果函数𝑘𝑞(𝑥)远大于原始分布函数 ̂𝑝(𝑥),拒绝率会比较高,采样效率会非常不理想.但要找到一个和 ̂𝑝(𝑥)比较接近的提议分布往往比较困难.特别是在高维空间中,其采样率会非常低,导致很难应用到实际问题中.

重要性采样(Importance sampling)
如果采样的目的是计算分布𝑝(𝑥)下函数𝑓(𝑥)的期望,那么实际上抽取的样本不需要严格服从分布 𝑝(𝑥).也可以通过另一个分布,即提议分布 𝑞(𝑥),直接采样并估计𝔼𝑝[𝑓(𝑥)].
函数𝑓(𝑥)在分布𝑝(𝑥)下的期望可以写为
Ep[f(x)]=xf(x)p(x)dx=xf(x)p(x)q(x)q(x)dx=xf(x)w(x)q(x)dx=Eq[f(x)w(x)]E_p[f(x)] = \int_x f(x)p(x)dx = \int_x f(x)\frac{p(x)}{q(x)}q(x)dx = \int_x f(x)w(x)q(x)dx = E_q[f(x)w(x)]
其中𝑤(𝑥)称为重要性权重.

重要性采样(Importance Sampling)是通过引入重要性权重,将分布 𝑝(𝑥)下𝑓(𝑥)的期望变为在分布𝑞(𝑥)下𝑓(𝑥)𝑤(𝑥)的期望,从而可以近似为
Ep[f(x)]=Eq[f(x)w(x)]=1Ni=1Nf(xi)w(xi)=1Ni=1Nf(xi)p(xi)q(xi)E_p[f(x)] = E_q[f(x)w(x)] =\frac{1}{N} \sum_{i=1}^N f(x_i)w(x_i) =\frac{1}{N} \sum_{i=1}^N f(x_i)\frac{p(x_i)}{q(x_i)}