Introduction

对概率的诠释有两大学派，一种是频率派另一种是贝叶斯派。后面我们对观测集采用下面记号：

XN×p=(x1,x2,⋯ ,xN)T,xi=(xi1,xi2,⋯ ,xip)TX_{N\times p}=(x_{1},x_{2},\cdots,x_{N})^{T},x_{i}=(x_{i1},x_{i2},\cdots,x_{ip})^{T}

这个记号表示有 NN 个样本，每个样本都是 pp 维向量。其中每个观测都是由 p(x∣θ)p(x|\theta) 生成的。

频率派的观点

p(x∣θ)p(x|\theta)中的 θ\theta 是一个常量。对于 NN 个观测来说观测集的概率为 p(X∣θ)=iid∏i=1Np(xi∣θ))p(X|\theta)\mathop{=}\limits _{iid}\prod\limits _{i=1}^{N}p(x_{i}|\theta)) 。为了求 θ\theta 的大小，我们采用最大对数似然MLE的方法：

θMLE=argmaxθlog⁡p(X∣θ)=iidargmaxθ∑i=1Nlog⁡p(xi∣θ)\theta_{MLE}=\mathop{argmax}\limits _{\theta}\log p(X|\theta)\mathop{=}\limits _{iid}\mathop{argmax}\limits _{\theta}\sum\limits _{i=1}^{N}\log p(x_{i}|\theta)

贝叶斯派的观点

贝叶斯派认为 p(x∣θ)p(x|\theta) 中的 θ\theta 不是一个常量。这个 θ\theta 满足一个预设的先验的分布 θ∼p(θ)\theta\sim p(\theta) 。于是根据贝叶斯定理依赖观测集参数的后验可以写成：

p(θ∣X)=p(X∣θ)⋅p(θ)p(X)=p(X∣θ)⋅p(θ)∫θp(X∣θ)⋅p(θ)dθp(\theta|X)=\frac{p(X|\theta)\cdot p(\theta)}{p(X)}=\frac{p(X|\theta)\cdot p(\theta)}{\int\limits _{\theta}p(X|\theta)\cdot p(\theta)d\theta}

为了求 θ\theta 的值，我们要最大化这个参数后验MAP：

θMAP=argmaxθp(θ∣X)=argmaxθp(X∣θ)⋅p(θ)\theta_{MAP}=\mathop{argmax}\limits _{\theta}p(\theta|X)=\mathop{argmax}\limits _{\theta}p(X|\theta)\cdot p(\theta)

其中第二个等号是由于分母和 θ\theta 没有关系。求解这个 θ\theta 值后计算p(X∣θ)⋅p(θ)∫θp(X∣θ)⋅p(θ)dθ\frac{p(X|\theta)\cdot p(\theta)}{\int\limits _{\theta}p(X|\theta)\cdot p(\theta)d\theta} ，就得到了参数的后验概率。其中 p(X∣θ)p(X|\theta) 叫似然，是我们的模型分布。得到了参数的后验分布后，我们可以将这个分布用于预测贝叶斯预测：

p(xnew∣X)=∫θp(xnew∣θ)⋅p(θ∣X)dθp(x_{new}|X)=\int\limits _{\theta}p(x_{new}|\theta)\cdot p(\theta|X)d\theta

其中积分中的被乘数是模型，乘数是后验分布。

小结

频率派和贝叶斯派分别给出了一系列的机器学习算法。频率派的观点导出了一系列的统计机器学习算法而贝叶斯派导出了概率图理论。在应用频率派的 MLE 方法时最优化理论占有重要地位。而贝叶斯派的算法无论是后验概率的建模还是应用这个后验进行推断时积分占有重要地位。因此采样积分方法如 MCMC 有很多应用。

MathBasics

高斯分布

一维情况 MLE

高斯分布在机器学习中占有举足轻重的作用。在 MLE 方法中：

θ=(μ,Σ)=(μ,σ2),θMLE=argmaxθlog⁡p(X∣θ)=iidargmaxθ∑i=1Nlog⁡p(xi∣θ)\theta=(\mu,\Sigma)=(\mu,\sigma^{2}),\theta_{MLE}=\mathop{argmax}\limits _{\theta}\log p(X|\theta)\mathop{=}\limits _{iid}\mathop{argmax}\limits _{\theta}\sum\limits _{i=1}^{N}\log p(x_{i}|\theta)

一般地，高斯分布的概率密度函数PDF写为：

p(x∣μ,Σ)=1(2π)p/2∣Σ∣1/2e−12(x−μ)TΣ−1(x−μ)p(x|\mu,\Sigma)=\frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}}e^{-\frac{1}{2}(x-\mu)^{T}\Sigma^{-1}(x-\mu)}

带入 MLE 中我们考虑一维的情况

log⁡p(X∣θ)=∑i=1Nlog⁡p(xi∣θ)=∑i=1Nlog⁡12πσexp⁡(−(xi−μ)2/2σ2)\log p(X|\theta)=\sum\limits _{i=1}^{N}\log p(x_{i}|\theta)=\sum\limits _{i=1}^{N}\log\frac{1}{\sqrt{2\pi}\sigma}\exp(-(x_{i}-\mu)^{2}/2\sigma^{2})

首先对 μ\mu 的极值可以得到 ：

μMLE=argmaxμlog⁡p(X∣θ)=argmaxμ∑i=1N(xi−μ)2\mu_{MLE}=\mathop{argmax}\limits _{\mu}\log p(X|\theta)=\mathop{argmax}\limits _{\mu}\sum\limits _{i=1}^{N}(x_{i}-\mu)^{2}

于是：

∂∂μ∑i=1N(xi−μ)2=0⟶μMLE=1N∑i=1Nxi\frac{\partial}{\partial\mu}\sum\limits _{i=1}^{N}(x_{i}-\mu)^{2}=0\longrightarrow\mu_{MLE}=\frac{1}{N}\sum\limits _{i=1}^{N}x_{i}

其次对 θ\theta 中的另一个参数 σ\sigma ，有：

σMLE=argmaxσlog⁡p(X∣θ)=argmaxσ∑i=1N[−log⁡σ−12σ2(xi−μ)2]=argminσ∑i=1N[log⁡σ+12σ2(xi−μ)2]\begin{align} \sigma_{MLE}=\mathop{argmax}\limits _{\sigma}\log p(X|\theta)&=\mathop{argmax}\limits _{\sigma}\sum\limits _{i=1}^{N}[-\log\sigma-\frac{1}{2\sigma^{2}}(x_{i}-\mu)^{2}]\nonumber\\ &=\mathop{argmin}\limits _{\sigma}\sum\limits _{i=1}^{N}[\log\sigma+\frac{1}{2\sigma^{2}}(x_{i}-\mu)^{2}] \end{align}

于是：

∂∂σ∑i=1N[log⁡σ+12σ2(xi−μ)2]=0⟶σMLE2=1N∑i=1N(xi−μ)2\frac{\partial}{\partial\sigma}\sum\limits _{i=1}^{N}[\log\sigma+\frac{1}{2\sigma^{2}}(x_{i}-\mu)^{2}]=0\longrightarrow\sigma_{MLE}^{2}=\frac{1}{N}\sum\limits _{i=1}^{N}(x_{i}-\mu)^{2}

值得注意的是，上面的推导中，首先对 μ\mu 求 MLE， 然后利用这个结果求 σMLE\sigma_{MLE} ，因此可以预期的是对数据集求期望时 ED[μMLE]\mathbb{E}_{\mathcal{D}}[\mu_{MLE}] 是无偏差的：

ED[μMLE]=ED[1N∑i=1Nxi]=1N∑i=1NED[xi]=μ\mathbb{E}_{\mathcal{D}}[\mu_{MLE}]=\mathbb{E}_{\mathcal{D}}[\frac{1}{N}\sum\limits _{i=1}^{N}x_{i}]=\frac{1}{N}\sum\limits _{i=1}^{N}\mathbb{E}_{\mathcal{D}}[x_{i}]=\mu

但是当对 σMLE\sigma_{MLE} 求 期望的时候由于使用了单个数据集的 μMLE\mu_{MLE}，因此对所有数据集求期望的时候我们会发现 σMLE\sigma_{MLE} 是 有偏的：

ED[σMLE2]=ED[1N∑i=1N(xi−μMLE)2]=ED[1N∑i=1N(xi2−2xiμMLE+μMLE2)=ED[1N∑i=1Nxi2−μMLE2]=ED[1N∑i=1Nxi2−μ2+μ2−μMLE2]=ED[1N∑i=1Nxi2−μ2]−ED[μMLE2−μ2]=σ2−(ED[μMLE2]−μ2)=σ2−(ED[μMLE2]−ED2[μMLE])=σ2−Var[μMLE]=σ2−Var[1N∑i=1Nxi]=σ2−1N2∑i=1NVar[xi]=N−1Nσ2\begin{align} \mathbb{E}_{\mathcal{D}}[\sigma_{MLE}^{2}]&=\mathbb{E}_{\mathcal{D}}[\frac{1}{N}\sum\limits _{i=1}^{N}(x_{i}-\mu_{MLE})^{2}]=\mathbb{E}_{\mathcal{D}}[\frac{1}{N}\sum\limits _{i=1}^{N}(x_{i}^{2}-2x_{i}\mu_{MLE}+\mu_{MLE}^{2})\nonumber \\&=\mathbb{E}_{\mathcal{D}}[\frac{1}{N}\sum\limits _{i=1}^{N}x_{i}^{2}-\mu_{MLE}^{2}]=\mathbb{E}_{\mathcal{D}}[\frac{1}{N}\sum\limits _{i=1}^{N}x_{i}^{2}-\mu^{2}+\mu^{2}-\mu_{MLE}^{2}]\nonumber\\ &= \mathbb{E}_{\mathcal{D}}[\frac{1}{N}\sum\limits _{i=1}^{N}x_{i}^{2}-\mu^{2}]-\mathbb{E}_{\mathcal{D}}[\mu_{MLE}^{2}-\mu^{2}]=\sigma^{2}-(\mathbb{E}_{\mathcal{D}}[\mu_{MLE}^{2}]-\mu^{2})\nonumber\\&=\sigma^{2}-(\mathbb{E}_{\mathcal{D}}[\mu_{MLE}^{2}]-\mathbb{E}_{\mathcal{D}}^{2}[\mu_{MLE}])=\sigma^{2}-Var[\mu_{MLE}]\nonumber\\&=\sigma^{2}-Var[\frac{1}{N}\sum\limits _{i=1}^{N}x_{i}]=\sigma^{2}-\frac{1}{N^{2}}\sum\limits _{i=1}^{N}Var[x_{i}]=\frac{N-1}{N}\sigma^{2} \end{align}

所以：

σ^2=1N−1∑i=1N(xi−μ)2\hat{\sigma}^{2}=\frac{1}{N-1}\sum\limits _{i=1}^{N}(x_{i}-\mu)^{2}

多维情况

多维高斯分布表达式为：

p(x∣μ,Σ)=1(2π)p/2∣Σ∣1/2e−12(x−μ)TΣ−1(x−μ)p(x|\mu,\Sigma)=\frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}}e^{-\frac{1}{2}(x-\mu)^{T}\Sigma^{-1}(x-\mu)}

其中 x,μ∈Rp,Σ∈Rp×px,\mu\in\mathbb{R}^{p},\Sigma\in\mathbb{R}^{p\times p} ，Σ\Sigma 为协方差矩阵，一般而言也是半正定矩阵。这里我们只考虑正定矩阵。首先我们处理指数上的数字，指数上的数字可以记为 xx 和 μ\mu 之间的马氏距离。对于对称的协方差矩阵可进行特征值分解，Σ=UΛUT=(u1,u2,⋯ ,up)diag(λi)(u1,u2,⋯ ,up)T=∑i=1puiλiuiT\Sigma=U\Lambda U^{T}=(u_{1},u_{2},\cdots,u_{p})diag(\lambda_{i})(u_{1},u_{2},\cdots,u_{p})^{T}=\sum\limits _{i=1}^{p}u_{i}\lambda_{i}u_{i}^{T} ，于是：

Σ−1=∑i=1pui1λiuiT\Sigma^{-1}=\sum\limits _{i=1}^{p}u_{i}\frac{1}{\lambda_{i}}u_{i}^{T}

Δ=(x−μ)TΣ−1(x−μ)=∑i=1p(x−μ)Tui1λiuiT(x−μ)=∑i=1pyi2λi\Delta=(x-\mu)^{T}\Sigma^{-1}(x-\mu)=\sum\limits _{i=1}^{p}(x-\mu)^{T}u_{i}\frac{1}{\lambda_{i}}u_{i}^{T}(x-\mu)=\sum\limits _{i=1}^{p}\frac{y_{i}^{2}}{\lambda_{i}}

我们注意到 yiy_{i} 是 x−μx-\mu 在特征向量 uiu_{i} 上的投影长度，因此上式子就是 Δ\Delta 取不同值时的同心椭圆。

下面我们看多维高斯模型在实际应用时的两个问题

  1. 参数 Σ,μ\Sigma,\mu 的自由度为 O(p2)O(p^{2}) 对于维度很高的数据其自由度太高。解决方案：高自由度的来源是 Σ\Sigma 有 p(p+1)2\frac{p(p+1)}{2} 个自由参数，可以假设其是对角矩阵，甚至在各向同性假设中假设其对角线上的元素都相同。前一种的算法有 Factor Analysis，后一种有概率 PCA(p-PCA) 。

  2. 第二个问题是单个高斯分布是单峰的，对有多个峰的数据分布不能得到好的结果。解决方案：高斯混合GMM 模型。

下面对多维高斯分布的常用定理进行介绍。

我们记 x=(x1,x2,⋯ ,xp)T=(xa,m×1,xb,n×1)T,μ=(μa,m×1,μb,n×1),Σ=(ΣaaΣabΣbaΣbb)x=(x_1, x_2,\cdots,x_p)^T=(x_{a,m\times 1}, x_{b,n\times1})^T,\mu=(\mu_{a,m\times1}, \mu_{b,n\times1}),\Sigma=\begin{pmatrix}\Sigma_{aa}&\Sigma_{ab}\\\Sigma_{ba}&\Sigma_{bb}\end{pmatrix}，已知 x∼N(μ,Σ)x\sim\mathcal{N}(\mu,\Sigma)。

首先是一个高斯分布的定理：

定理：已知 x∼N(μ,Σ),y∼Ax+bx\sim\mathcal{N}(\mu,\Sigma), y\sim Ax+b，那么 y∼N(Aμ+b,AΣAT)y\sim\mathcal{N}(A\mu+b, A\Sigma A^T)。

证明：E[y]=E[Ax+b]=AE[x]+b=Aμ+b\mathbb{E}[y]=\mathbb{E}[Ax+b]=A\mathbb{E}[x]+b=A\mu+b，Var[y]=Var[Ax+b]=Var[Ax]=A⋅Var[x]⋅ATVar[y]=Var[Ax+b]=Var[Ax]=A\cdot Var[x]\cdot A^T。

下面利用这个定理得到 p(xa),p(xb),p(xa∣xb),p(xb∣xa)p(x_a),p(x_b),p(x_a|x_b),p(x_b|x_a) 这四个量。

  1. xa=(Im×mOm×n))(xaxb)x_a=\begin{pmatrix}\mathbb{I}_{m\times m}&\mathbb{O}_{m\times n})\end{pmatrix}\begin{pmatrix}x_a\\x_b\end{pmatrix}，代入定理中得到：

    E[xa]=(IO)(μaμb)=μaVar[xa]=(IO)(ΣaaΣabΣbaΣbb)(IO)=Σaa\mathbb{E}[x_a]=\begin{pmatrix}\mathbb{I}&\mathbb{O}\end{pmatrix}\begin{pmatrix}\mu_a\\\mu_b\end{pmatrix}=\mu_a\\ Var[x_a]=\begin{pmatrix}\mathbb{I}&\mathbb{O}\end{pmatrix}\begin{pmatrix}\Sigma_{aa}&\Sigma_{ab}\\\Sigma_{ba}&\Sigma_{bb}\end{pmatrix}\begin{pmatrix}\mathbb{I}\\\mathbb{O}\end{pmatrix}=\Sigma_{aa}

    所以 xa∼N(μa,Σaa)x_a\sim\mathcal{N}(\mu_a,\Sigma_{aa})。

  2. 同样的，xb∼N(μb,Σbb)x_b\sim\mathcal{N}(\mu_b,\Sigma_{bb})。

  3. 对于两个条件概率，我们引入三个量：

    xb⋅a=xb−ΣbaΣaa−1xaμb⋅a=μb−ΣbaΣaa−1μaΣbb⋅a=Σbb−ΣbaΣaa−1Σabx_{b\cdot a}=x_b-\Sigma_{ba}\Sigma_{aa}^{-1}x_a\\ \mu_{b\cdot a}=\mu_b-\Sigma_{ba}\Sigma_{aa}^{-1}\mu_a\\ \Sigma_{bb\cdot a}=\Sigma_{bb}-\Sigma_{ba}\Sigma_{aa}^{-1}\Sigma_{ab}

    特别的，最后一个式子叫做 Σbb\Sigma_{bb} 的 Schur Complementary。可以看到：

    xb⋅a=(−ΣbaΣaa−1In×n)(xaxb)x_{b\cdot a}=\begin{pmatrix}-\Sigma_{ba}\Sigma_{aa}^{-1}&\mathbb{I}_{n\times n}\end{pmatrix}\begin{pmatrix}x_a\\x_b\end{pmatrix}

    所以：

    E[xb⋅a]=(−ΣbaΣaa−1In×n)(μaμb)=μb⋅aVar[xb⋅a]=(−ΣbaΣaa−1In×n)(ΣaaΣabΣbaΣbb)(−Σaa−1ΣbaTIn×n)=Σbb⋅a\mathbb{E}[x_{b\cdot a}]=\begin{pmatrix}-\Sigma_{ba}\Sigma_{aa}^{-1}&\mathbb{I}_{n\times n}\end{pmatrix}\begin{pmatrix}\mu_a\\\mu_b\end{pmatrix}=\mu_{b\cdot a}\\ Var[x_{b\cdot a}]=\begin{pmatrix}-\Sigma_{ba}\Sigma_{aa}^{-1}&\mathbb{I}_{n\times n}\end{pmatrix}\begin{pmatrix}\Sigma_{aa}&\Sigma_{ab}\\\Sigma_{ba}&\Sigma_{bb}\end{pmatrix}\begin{pmatrix}-\Sigma_{aa}^{-1}\Sigma_{ba}^T\\\mathbb{I}_{n\times n}\end{pmatrix}=\Sigma_{bb\cdot a}

    利用这三个量可以得到 xb=xb⋅a+ΣbaΣaa−1xax_b=x_{b\cdot a}+\Sigma_{ba}\Sigma_{aa}^{-1}x_a。因此：

    E[xb∣xa]=μb⋅a+ΣbaΣaa−1xa\mathbb{E}[x_b|x_a]=\mu_{b\cdot a}+\Sigma_{ba}\Sigma_{aa}^{-1}x_a

    Var[xb∣xa]=Σbb⋅aVar[x_b|x_a]=\Sigma_{bb\cdot a}

    这里同样用到了定理。

  4. 同样：

    xa⋅b=xa−ΣabΣbb−1xbμa⋅b=μa−ΣabΣbb−1μbΣaa⋅b=Σaa−ΣabΣbb−1Σbax_{a\cdot b}=x_a-\Sigma_{ab}\Sigma_{bb}^{-1}x_b\\ \mu_{a\cdot b}=\mu_a-\Sigma_{ab}\Sigma_{bb}^{-1}\mu_b\\ \Sigma_{aa\cdot b}=\Sigma_{aa}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba}

    所以：

    E[xa∣xb]=μa⋅b+ΣabΣbb−1xb\mathbb{E}[x_a|x_b]=\mu_{a\cdot b}+\Sigma_{ab}\Sigma_{bb}^{-1}x_b

    Var[xa∣xb]=Σaa⋅bVar[x_a|x_b]=\Sigma_{aa\cdot b}

下面利用上边四个量，求解线性模型：

已知：p(x)=N(μ,Λ−1),p(y∣x)=N(Ax+b,L−1)p(x)=\mathcal{N}(\mu,\Lambda^{-1}),p(y|x)=\mathcal{N}(Ax+b,L^{-1})，求解：p(y),p(x∣y)p(y),p(x|y)。

解：令 y=Ax+b+ϵ,ϵ∼N(0,L−1)y=Ax+b+\epsilon,\epsilon\sim\mathcal{N}(0,L^{-1})，所以 E[y]=E[Ax+b+ϵ]=Aμ+b\mathbb{E}[y]=\mathbb{E}[Ax+b+\epsilon]=A\mu+b，Var[y]=AΛ−1AT+L−1Var[y]=A \Lambda^{-1}A^T+L^{-1}，因此：

p(y)=N(Aμ+b,L−1+AΛ−1AT)p(y)=\mathcal{N}(A\mu+b,L^{-1}+A\Lambda^{-1}A^T)

引入 z=(xy)z=\begin{pmatrix}x\\y\end{pmatrix}，我们可以得到 Cov[x,y]=E[(x−E[x])(y−E[y])T]Cov[x,y]=\mathbb{E}[(x-\mathbb{E}[x])(y-\mathbb{E}[y])^T]。对于这个协方差可以直接计算：

Cov(x,y)=E[(x−μ)(Ax−Aμ+ϵ)T]=E[(x−μ)(x−μ)TAT]=Var[x]AT=Λ−1AT\begin{align} Cov(x,y)&=\mathbb{E}[(x-\mu)(Ax-A\mu+\epsilon)^T]=\mathbb{E}[(x-\mu)(x-\mu)^TA^T]=Var[x]A^T=\Lambda^{-1}A^T \end{align}

注意到协方差矩阵的对称性，所以 p(z)=N(μAμ+b),(Λ−1Λ−1ATAΛ−1L−1+AΛ−1AT))p(z)=\mathcal{N}\begin{pmatrix}\mu\\A\mu+b\end{pmatrix},\begin{pmatrix}\Lambda^{-1}&\Lambda^{-1}A^T\\A\Lambda^{-1}&L^{-1}+A\Lambda^{-1}A^T\end{pmatrix})。根据之前的公式，我们可以得到：

E[x∣y]=μ+Λ−1AT(L−1+AΛ−1AT)−1(y−Aμ−b)\mathbb{E}[x|y]=\mu+\Lambda^{-1}A^T(L^{-1}+A\Lambda^{-1}A^T)^{-1}(y-A\mu-b)

Var[x∣y]=Λ−1−Λ−1AT(L−1+AΛ−1AT)−1AΛ−1Var[x|y]=\Lambda^{-1}-\Lambda^{-1}A^T(L^{-1}+A\Lambda^{-1}A^T)^{-1}A\Lambda^{-1}