Bayesian signal processing

该博客为张颢老师所讲的现代信号处理的贝叶斯信号处理部分的课程个人学习笔记,涵盖贝叶斯公式的含义,贝叶斯学派与频率学派的关系,贝叶斯统计模型的目标函数,蒙特卡洛采样和贝叶斯滤波的内容。

Introduce to bayes method

Bayes and frequency

使用贝叶斯公式并不意味着在使用贝叶斯方法,频率学派也会经常使用贝叶斯公式,早在贝叶斯之前,贝叶斯公式就已经在Laplce等数学家中非常娴熟的应用,贝叶斯方法在思想上面的含义要远远大于数学形式上的贝叶斯公式。

摘抄知乎上Xiangyu Wang的回答：频率学派和贝叶斯学派最大的差别其实产生于对参数空间的认知上。所谓参数空间,就是你关心的那个参数可能的取值范围。频率学派（其实就是当年的Fisher）并不关心参数空间的所有细节,他们相信数据都是在这个空间里的“某个”参数值下产生的（虽然你不知道那个值是啥）,所以他们的方法论一开始就是从“哪个值最有可能是真实值”这个角度出发的。于是就有了最大似然（maximum likelihood）以及置信区间（confidence interval）这样的东西,你从名字就可以看出来他们关心的就是我有多大把握去圈出那个唯一的真实参数。而贝叶斯学派恰恰相反,他们关心参数空间里的每一个值,因为他们觉得我们又没有上帝视角,怎么可能知道哪个值是真的呢？所以参数空间里的每个值都有可能是真实模型使用的值,区别只是概率不同而已。于是他们才会引入先验分布（prior distribution）和后验分布（posterior distribution）这样的概念来设法找出参数空间上的每个值的概率。最好诠释这种差别的例子就是想象如果你的后验分布是双峰的,频率学派的方法会去选这两个峰当中较高的那一个对应的值作为他们的最好猜测,而贝叶斯学派则会同时报告这两个值,并给出对应的概率。

假如我们有一个统计模型$f(x,\theta)$去估计原始数据$x$的分布,频率学派认为$\theta$是固定的,尽管我们现在还没有办法得到,贝叶斯学派认为$\theta$是随机的,是服从一个概率分布的随机变量。

Bayes: The Science of cognitive law

贝叶斯的观点是符合人脑的认知规律的,也是现代认知心理学的基础：在认识一件事物之前,我们会先对这个事物有一定的认知,在与事物接触的过程中,我们的认知会逐渐的发生变化,或加强我们的认知,或减弱我们的认知,我们的认知从而得到了进步,这与马原中所说的认知随着实践螺旋式上升也有异曲同工之妙。贝叶斯公式就是表述了这个认知的过程：

$$ P(\theta \mid X)=\frac{P(X \mid \theta) P(\theta)}{P(X)} $$

贝叶斯学派通过引入三个名词对贝叶斯公式进行解释：

名词	公式
先验分布(Prior)	$P(\theta)$
似然(Likelihood)	$P(X \mid \theta)$
后验分布(Postertor)	$P(\theta \mid X)$

先验分布就是在接触事物数据之前我们已经有的认知,似然是在接触事物的过程中,在给定先验知识的条件下我们对数据的认知,后验分布是我们在数据的基础上对事物的认知,机器学习中的“Learning”这一概念,很大程度上就是贝叶斯方法中的使用后验分布去调整先验分布,当数据量足够大的时候,$P(\theta)$就越来越接近真实的概率模型,另外,分母$P(X)$经常被成为贝叶斯因子。

我们再来观察特殊的情况。如果数据与先验独立,即$X \perp \theta$,有$P(\theta \mid X)=P(\theta)$,即先验等于后验,也就是说在接触数据的过程中,我们并没有学到有用的信息。如果$P(\theta)=0$,我们的先验分布和后验分布会一直等于0,所以在很多机器学习算法中设置初始迭代解的时候,会避开零概率的设置。

Example: Naive Bayes

朴素贝叶斯是一种简单的机器学习分类器,分类任务就是在给定某些特征 $\lbrace F_i \rbrace_{i=1}^N$(Features)的条件下,判断某一事物的分类$C$(Catogories),即我们需要获得：$P(C \mid F_1:F_N)$,数据是给定条件,分类是随机变量,这就是贝叶斯方法中所说的后验分布,我们可以使用先验分布和似然来表示该后验分布：

$$ P(C \mid F_1:F_N)=\frac{P(F_1:F_N \mid C) P(C)}{P(F_1:F_N)} $$

对不同特征求给定条件的联合分布是比较困难的,naive bayes做了一个假设,在给定分类条件下,不同特征之间相互独立,即我们可以重写上述公式为

$$ P(C \mid F_1:F_N)=\frac{\prod_{i=1}^N P(F_i \mid C) P(C)}{P(F_1:F_N)} $$

分母可以通过统计样本得出,但是并不重要,因为所有的后验分布共享一个分母,当分类器判别样本时,都具有相同的权重。

Example:Use bayes to process Gaussian distribution

考虑正态似然与正态先验：$X_1,\ldots,X_N\mid\mu\overset{\mathrm{iid}}{\sim}\mathcal{N}(\mu,\sigma_x^2)$，且先验为 $\mu\sim\mathcal{N}(\mu_\theta,\sigma_\mu^2)$。由贝叶斯公式，后验密度满足

$$ p(\mu\mid x_{1:N})\propto p(x_{1:N}\mid\mu)p(\mu). $$

忽略与 $\mu$ 无关的项，后验的对数密度为

$$ \begin{aligned} \log p(\mu\mid x_{1:N}) &=-\frac{1}{2}\left(\frac{N}{\sigma_x^2}+\frac{1}{\sigma_\mu^2}\right)\mu^2 +\left(\frac{\sum_{i=1}^N x_i}{\sigma_x^2}+\frac{\mu_\theta}{\sigma_\mu^2}\right)\mu+C\\ &=-\frac{1}{2v_N}(\mu-m_N)^2+C', \end{aligned} $$

其中

$$ v_N=\left(\frac{N}{\sigma_x^2}+\frac{1}{\sigma_\mu^2}\right)^{-1},\qquad m_N=v_N\left(\frac{\sum_{i=1}^N x_i}{\sigma_x^2}+\frac{\mu_\theta}{\sigma_\mu^2}\right). $$

因此后验分布仍为正态分布：

$$ \mu\mid x_{1:N}\sim\mathcal{N}(m_N,v_N). $$

令 $\bar X=N^{-1}\sum_{i=1}^N x_i$，后验均值也可写成数据均值与先验均值的加权平均：

$$ m_N=\frac{N/\sigma_x^2}{N/\sigma_x^2+1/\sigma_\mu^2}\,\bar X +\frac{1/\sigma_\mu^2}{N/\sigma_x^2+1/\sigma_\mu^2}\,\mu_\theta. $$

样本量增大时，数据在后验均值中的权重增大；似然方差增大时，先验均值的权重相对增大。

Statistical decision theory

贝叶斯方法的核心观点就是利用后验分布去学习统计模型$\hat{\theta}(x)$,去对先验分布进行优化,那么就需要我们定义一种度量$d(\hat{\theta}(x),\theta)$,从而通过最小化这种度量的方式完成先验分布的优化,$d(\hat{\theta}(x),\theta)$是随机的,所以我们希望消除他的随机性求平均距离,即我们的目标为

$$ \hat{\theta}(x) = \underset{\hat{\theta}}{\arg \min }\mathbb{E}_{x,\theta}[d(\hat{\theta}(x),\theta) ] $$

平均距离$\mathbb{E}_{x,\theta}[d(\hat{\theta}(x),\theta) ]$ 一般被成为风险(Risk),该目标函数也被称为风险极小化,频率学派和贝叶斯学派对该目标函数有不同方式的处理：

$$ \begin{aligned} E[d(\hat{\theta}(x), \theta)]&=\iint d(\theta, \hat{\theta}(x)) P_{x, \theta}\left(x, \theta\right) d x d\theta\\ &=\iint d(\theta, \hat{\theta}(x)) P_{x\mid \theta}\left(x\mid \theta\right) d xP(\theta) d\theta\\ &=\iint d(\theta, \hat{\theta}(x)) P_{x\mid \theta}\left(\theta \mid x\right) d \theta P(x) dx\\ \end{aligned} $$

可见,频率学派尽可能的削弱$\theta$的随机性,因为他们认为$\theta$是确定的（所以在频率学派中有极大化似然这种方法）,但是频率学派种所需要的似然往往是不可以直接计算的,因此有了经验风险极小化(Empirical Risk Minimization,ERM)理论,该理论一般将样本均值取代理论均值,即

$$ E[d(\hat{\theta}(x), \theta)]=\frac{1}{N}\sum_{k=1}^Nd(\theta_k, \hat{\theta}_k(x)) $$

贝叶斯学派则是尽可能削弱数据的随机性,认为已经发生的实验就是确定的,将问题的重心转化到后验分布的计算上,在后面我们会看到,使用贝叶斯方法处理统计模型的关键就是得到后验分布以及如何处理后验分布。

Choice of measurement

1. MSE

平方损失对应的贝叶斯风险为

$$ R(\hat\theta\mid x)=\int(\theta-\hat\theta)^2p(\theta\mid x)\,d\theta. $$

使风险最小的估计是后验均值：

$$ \hat\theta_{\mathrm{MSE}}=\mathbb{E}[\theta\mid x] =\int\theta p(\theta\mid x)\,d\theta. $$

2. MAE

绝对损失对应的风险为

$$ R(\hat\theta\mid x)=\int|\theta-\hat\theta|p(\theta\mid x)\,d\theta. $$

其最小值由后验中位数取得。对连续后验分布，最优点满足

$$ \int_{-\infty}^{\hat\theta}p(\theta\mid x)\,d\theta =\int_{\hat\theta}^{+\infty}p(\theta\mid x)\,d\theta=\frac12. $$

3. MAP

最大后验估计选择后验密度的众数：

$$ \hat\theta_{\mathrm{MAP}}\in\operatorname*{arg\,max}_{\theta}p(\theta\mid x). $$

它也可由小邻域的零一损失得到。令 $d_\delta(\theta,\hat\theta)=\mathbf{1}\{|\theta-\hat\theta|>\delta\}$，则当 $\delta$ 足够小时，最小化后验风险等价于最大化 $p(\hat\theta\mid x)$。

Example:Use bayes to process Bernoulli distribution

假设我们有先验分布$\theta$,似然$X\sim B(n,\theta)$,目标计算$\theta$的后验分布以及在均方意义下的最优估计。

由Bayes公式,有

$$ \begin{aligned} P(\theta \mid x) & =\frac{P(x \mid \theta) P(\theta)}{P(x)} \\ & =\frac{P(x \mid \theta) P(\theta)}{\int P(x, \theta) d \theta} \end{aligned} $$

接下来,我们考虑当先验分布服从标准均匀分布的情况下后验概率的计算,即计算

$$ P(\theta \mid x)=\frac{\theta^x(1-\theta)^{n-x}}{\int_0^1 \theta^x(1-\theta)^{n-x} d \theta} $$

在计算后验概率之前,我们需要引入Beta分布来辅助我们的贝叶斯因子的计算,Beta分布的pdf为

$$ \begin{aligned} f(x ; \alpha, \beta) & =\frac{x^{\alpha-1}(1-x)^{\beta-1}}{\int_0^1 u^{\alpha-1}(1-u)^{\beta-1} d u} \\ & =\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha) \Gamma(\beta)} x^{\alpha-1}(1-x)^{\beta-1} \\ & =\frac{1}{\mathrm{~B}(\alpha, \beta)} x^{\alpha-1}(1-x)^{\beta-1} \end{aligned} $$

其中,$\alpha,\beta>0$,$\Gamma(x)$是Gamma函数,有性质

$$ \Gamma(x)=\int_0^{+\infty} t^{x-1} e^{-t} d t(x>0)=(x-1)! $$

Beta函数$B(\alpha, \beta)$的定义为$B(\alpha, \beta)=\int_0^1 t^{\alpha-1}(1-t)^{\beta-1} \mathbf{d} t$,现在我们来证明$\mathrm{~B}(\alpha, \beta)=\frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha+\beta)}$:

观察Beta函数的定义,类似卷积的定义,注意到：

$$ \begin{aligned} B(\alpha, \beta)&=\int_0^1 t^{\alpha-1}(x-t)^{\beta-1} \mathbf{d} t\mid_{x=1} \\&=f(x)*g(x)\mid_{x=1} \end{aligned} $$

在上述公式中,我们定义了两个幂函数$f(x): =x^{\alpha-1},g(x):=x^{\beta-1},$这两个幂函数的支撑集应该是[0,1],这样才能满足积分的上下限要求,现在对两个幂级数的卷积进行Laplce变换：

$$ \mathcal{L}\left\{x^{\alpha-1} * x^{\beta-1}\right\}=\mathcal{L}\left\{x^{\alpha-1}\right\} \cdot \mathcal{L}\left\{x^{\beta-1}\right\} $$

而幂级数的Laplace变换可以和Gamma函数联系起来,即

$$ \begin{aligned} \mathcal{L}\left\{t^a\right\}&=\int_0^{\infty} t^a e^{-s t} d t \\ &=\frac{1}{s^{a+1}} \int_0^{\infty}(s t)^a e^{-(s t)}(s d t) \\ &=\frac{1}{s^{a+1}} \int_0^{\infty} \tau^{(a+1)-1} e^{-\tau} d \tau \\ &=\frac{\Gamma(a+1)}{s^{a+1}} \end{aligned} $$

继续我们上述的Laplace变换,有

$$ \begin{aligned} \mathcal{L}\left\{x^{\alpha-1}\right\} \cdot \mathcal{L}\left\{x^{\beta-1}\right\}&=\frac{\Gamma(\alpha) \Gamma(\beta)}{s^{\alpha+\beta}} \\ &=\frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha+\beta)}\frac{\Gamma(\alpha+\beta-1+1)}{s^{(\alpha+\beta-1)+1}}\\ &=\frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha+\beta)} \mathcal{L}\left\{x^{\alpha+\beta-1}\right\} \end{aligned} $$

这时令$x=1$,便可以得到$\mathrm{~B}(\alpha, \beta)=\frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha+\beta)}$,证毕。

现在我们使用Beta函数来计算后验分布,注意,如果令$\alpha = x+1,\beta = n-x+1$,后验分布的贝叶斯因子为$B(\alpha,\beta)=B(x+1,n-x+1)$,那么有

$$ \begin{aligned} \int_0^1 \theta^x(1-\theta)^{n-x} d \theta & =B(x+1, n-x+1) \\ & =\frac{\Gamma(x+1) \Gamma(n-x+1)}{\Gamma(n+2)} \\ & =\frac{x!(n-x)!}{(n+1)!} \end{aligned} $$

则我们计算出来的后验分布最终形式为

$$ P(\theta, x)=\frac{(n+1)!}{x!(n-x)!} \theta^x(1-\theta)^{n-x} $$

可见该后验分布实际上是服从Beta分布的。根据MSE意义下的最优估计,有

$$ \begin{aligned} E[\theta \mid x] & =\int_0^1 \theta \frac{(n+1)!}{x!(n-x)!} \theta^x(1-\theta)^{n-x} d \theta \\ & =\frac{(n+1)!}{x!(n-x)!} \int_0^1 \theta^{x+1}(1-\theta)^{n-x} d \theta \\ & =\frac{(n+1)!}{x!(n-x)!} B(x+2, n-x+1) \\ & =\frac{(n+1)!}{n!(n-x)!} \frac{(x+1)!(n-x)!}{(n+2)!} \\ & =\frac{x+1}{n+2} \end{aligned} $$

至此,我们求得了后验分布和MSE意义下的模型最优估计。如果再进一步的处理该条件期望,我们可以得到线型凸组合的形式：

$$ E[\theta \mid x]=\frac{n}{n+2}\frac{x}{n}+\frac{2}{n+2}\frac{1}{2} $$

$\frac{x}{n}$反应了似然意义上的某种均值,$\frac{1}{2}$为先验分布的均值,这也显示了在伯努利实验的情况下,最优估计的过程是不断线性凸组合的过程。

Example:Conjugate prior

在上个例子中,我们一开始假设先验分布服从标准均匀分布,但是我们根据伯努利的似然求出后验分布之后,发现后验分布并不服从标准均匀分布,而是服从Beta分布,这说明我们的先验分布并没有反映事物的客观性质,在统计学中,如果我们的先验分布和后验分布的分布类型一致,我们就称该先验分布为共轭先验。

在一些比较简单的概率模型中,共轭先验是已知的并且可以提前制定的,例如当似然服从伯努利分布时,共轭先验为Beta分布,但是更多的(绝大多数的)概率模型的共轭先验是难以得到的。我们沿用上一个例子的情景,令先验分布为Beta分布,即

$$ P(\theta)=\frac{(a+b+1)!}{(a-1)!(b-1)!}\theta^{a-1}(1-\theta)^{b-1} $$

使用与上一个例子同样的方法去求后验分布和条件期望,最终可以求得条件期望为

$$ E[\theta \mid x] = \frac{x+a}{n+a+b} $$

同样的,该条件期望也可以写成线型图组合的形式,即

$$ E[\theta \mid x]=\frac{a+b}{n+a+b}\frac{a}{a+b}+\frac{n}{n+a+b}\frac{x}{n} $$

依然是某种实验中的均值与先验分布均值的线性凸组合,从我们在上一个章节中的Gauss分布似然到这个章节的伯努利分布似然,最后都得到了线性凸组合的学习形式,这可能并不是一种偶然,这些简单的例子给我们启发：在高维空间中,复杂分布的学习过程也可能是迭代的线性凸组合。

Sampling

使用Bayes方法的核心在于后验分布的处理，但是后验分布的形式往往是极其复杂的，极大多数的后验分布都是难以去积分甚至不可积的，统计中常用的方法就是生成符合后验分布的样本，从而利用样本代替解析表达的分布进行计算，这种生成样本的方法又被称作采样。

模拟是指把某一现实的或抽象的系统的某种特征或部分状态，用另一系统（称为模拟模型）来代替或近似。为了解决某问题，把它变成一个概率模型的求解问题，然后产生符合模型的大量随机数，对产生的随机数进行分析从而求解问题，这种方法叫做随机模拟方法，又称为蒙特卡洛(Monte Carlo)方法。

Inverse distribution function

逆分布采样的描述很简单，目标生成随机变量$X\sim F(X)$的样本($F(x)$为累计分布函数)，只需要生成随机样本$Y\sim F^{-1}(U)$即可，$Y$即想要得到的分布，$U$为服从标准正态分布的随机变量。

这样来看，生成随机变量的样本的过程就很简单，第一步，生成均匀分布的样本(可以利用线性同余法等方法生成)，第二步，求得目标分布累积分布函数的逆，第三步，代入累积分布函数的逆求得目标样本。

现在简单的证明一下随机变量$Y$与随机变量$X$服从同一个分布，求$Y$的累积分布函数：

$$ \begin{aligned} F_Y(y)& =P(y\leq Y)\\ & =P(F^{-1}(u)\leq Y)\\ & =P(u\leq F(Y))\\ & =F(Y) \end{aligned} $$

可见，$Y$的累积分布函数实际上就是$X$的累积分布函数，所以它们两个是同一个随机变量，如此，对于易于求得逆累积分布的随机变量，我们便可以通过这种方法求得符合该分布的样本，例如指数分布、泊松分布。

Example:Sampling for Gaussian Distribution Guassian分布的逆分布函数并不好求，虽然我们可以通过中心极限定理生成Gaussian分布，但那种方法并不属于逆分布生成，Box-Muller算法是一种经典的Gaussian分布生成算法，使用了极坐标变换来处理Gaussian分布难以求逆的问题，接下来我们介绍这种方法。

假设我们目标求得两个独立同分布的标准高斯随机变量，$X_1,X_2\overset{\mathrm{iid}}{\sim} N(0,1)$，则它们的联合概率密度为

$$ \begin{aligned} f_{X Y}(x, y) & =\frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{x^2}{2 \sigma^2}\right) \frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{y^2}{2 \sigma^2}\right) \\ & =\frac{1}{2 \pi } \exp \left(-\frac{x^2+y^2}{2 }\right) \end{aligned} $$

对该联合概率密度进行极坐标变换，有

$$ \begin{aligned} f_{X Y}(x, y) &=f_{XY}(\rho \cos\theta,\rho \sin\theta)\\ &=\frac{1}{2 \pi } \exp \left(-\frac{\rho^2}{2 }\right)\rho \end{aligned} $$

Accept-reject sampling

Importance sampling

Markov chain Monte Carlo(MCMC)

这三个章节的内容可以看另一篇博客.