What is Bayesian?

贝叶斯学派认为要估计的参数并不是一个具体值, 它也是一个随机变量. 而我们的目标就是估计这个随机变量的分布.

首先在估计之前, 我们会先形成一个对这个参数的先验(prior): 记为$\pi(\theta)$. 对于给定的数据$D$, 我们可以认为这一组数据是某个Draw: $\theta\sim \pi(\theta)$下得出的. 因此在这个Draw下得出这组数据的概率就是: $P(D|\theta)=\Pi_{i=1}^{n}P(D_{i},\theta)$. Partiallig out $\theta$, 有: $P(D)=\int P(D|\theta)\pi(\theta)d\theta$. 最后据贝叶斯公式, 可得后验分布(posterior): $\pi(\theta|x)=\frac{P(D|\theta)\pi(\theta)}{P(D)}\propto P(D|\theta)\pi(\theta)$. 能只关注分子是因为分母是个常数, 与$\theta$无关..

Application

常见的应用就是预测和后验预期

Prediction

我用不到..

MCMC

对于一些事件, 其后验分布可能会很复杂, 完全无法写出解析形式, 那么我们可以用MCMC方法来生成一系列服从后验分布的数据, 以此来推测后验分布.

MC (蒙特卡洛模拟法)

如果知道了某个分布形式, 那么可以通过抽象的方式获得其期望.

步骤很简单: 对于$X\sim P(X)$, 生成$X={x^{(1)},x^{(2)},x^{(3)},…x^{(R)}}$, for $x^{(r)}\sim P(X)$. 如果我们想知道某个函数$f(X)$在$P(X)$下的均值, 那么当样本量够大的时候, 有: $\hat{E}[f(X)] = \frac{1}{R} \sum_{r=1}^{R} f\bigl(x^{(r)}\bigr)\to_{p} \int f(x)P(x)dx$.

MC 显而易见地有一些性质:1. MC无偏; 2. 其方差$=O(\frac{1}{R})$.

Application: MC估计$\pi$的大小

Def: $P(x,y)=\left\{\begin{matrix}1, &x\in(0,1),y\in(0,1) \ 0, &others\end{matrix}\right.$, 有: $\pi = 4\int\int 1[x^2+y^2<1]P(x,y)dxdy \approx 4\tfrac1R \sum_{r=1}^{R}1[x^{(r)2}+y^{(r)2}<1]$.

Inverse CDF

当我们想通过MC来generate data的时候, 我们有一系列方法. 其中比较常用的就是Inverse CDF了:

令$F(X)$为CDF, 如果其严格单调且可逆, 则$F^{-1}(u)\in [0,1]$. 我们在$[0,1]$上生成一个均匀分布的随机数$U$, 有: $X=F^{-1}(U)$, 则我们将生成X的任务变成了生成$u\sim U(0,1)$. 假设生成的是$u={u^{(1)},u^{(2)},…,u^{(R)}}$, 则根据$X=F^{-1}(U)$ 可以计算得到一系列 ${x^{(1)},…,x^{(R)}}$. 近似服从$f(X)$.

Transformation

一些分布可以用其他分布来模拟: 常见的有: log-normal, Beta-Gamma, Truncated Normal.

Log Normal: 对于$X\sim N(\mu,\sigma^2)$, 指数化得到$Y\sim Lognormal(\mu,\sigma^2)$. 相当于先生成一系列X, 然后再取其指数. 这样得到的数据全都是正的.

Beta-Gamma: 如果$X\sim \Gamma(\alpha,1)$ and $Y\sim \Gamma(\beta,1)$ 且X,Y 相互独立, 那么$Z=\frac{X}{X+Y}\sim B (\alpha,\beta)$. 这样能直接简单的用两个Gamma分布变量模拟Beta分布变量. (pdf of Z: $f(z)=\frac{\Gamma(X+Y)}{\Gamma(X)\Gamma(Y)}z^{\alpha-1}(1-z)^{\beta-1}$)

Truncated Normal: 对于$X\sim N(\mu,\delta^2)$, 只保留$[a,b]$部分的正态分布, 有:pdf: $f_{T}(x)=\frac{\phi(x)}{\Phi(b)-\Phi(a)}$, $x\in [a,b]$. 截断正态分布的生成比较麻烦一些, 见下: (还是直接用Robert 算法吧)

Sum of RV

一些随机变量的和为另一随机变量的分布.

常见的Normal-Chi Square-F dist就不说了, 说个不常见的: exponential – Gamma

Accept-Reject Algorithm

这一套算法将$f_{X}$分布(目标分布)通过另一个分布(提议分布)$f_{Y}$来模拟:

Def $X\sim f_{X}$ and $Y\sim f_{Y}$, and support of Y contains that of X(指Y可取到的值比X更广), 有: def $M=\sup_{x}\frac{f_{X}(x)}{f_{Y}(x)}$, (看似玄乎, 实际就是要求$f_{X}(x)=O(f_{Y}(x))$, 这个要求相当于$f_{X}(x)\le Mf_{Y}(x)$).

则: 从$f_{Y}(x)$中生成一个x, 然后从$u\sim U(0,1)$中生成一个u, 如果$u\le \frac{f_{X}}{Mf_{Y}}$, 那么就接受x, 否则就抛弃x. 最终被接受的样本服从$f_{X}$.

接受率在30%-40%为佳, 否则会浪费算力. 此外, 这种方法在高维分布中效果有限.

改进方法: Adaptive Rejection Sampling, Adaptive Rejection Metropolis Sampling

Importance Sampling

如果我们想估计期望值, 那么在ARA中拒绝一些x会浪费算力, 那么我们想办法使用这些被拒绝的x:

比如我们对$h(x)$感兴趣, 那么对于$\int h(x)f_{X}(x)dx$, 有:

对于$\int h(x)\frac{f_{X}(x)}{f_{Y}(x)}f_{Y}(x)dx$, 可以直接MC: $\approx \frac{1}{R}\sum h(x^{r})\frac{f_{X}(x^{r})}{f_{Y}(x^{r})}$. 这里是把$f_{Y}$视作前面这个组合函数的分布了, 秒啊.

如果想求期望, 则有: $\hat{E}[f(x)] = \frac{\sum_{r=1}^{R} f(x^{(r)})\,w(x^{(r)})}{\sum_{r=1}^{R} w(x^{(r)})}$, for $w(x)=\frac{f_{X}(x)}{f_{Y}(x)}$ (期望必须归一化吗?)

核心思想是：通过从一个更容易采样的辅助分布（称为提议分布或重要性分布）中采样，再利用适当的权重来校正样本与目标分布之间的差异，从而近似估计目标分布下的期望或积分。

Slice Sampling

额.. 好抽象,等下再看

MC (马尔可夫链)

它是一系列随机变量: $X^{(0)},X^{(1)},X^{(2)},…X^{(n)}$, s.t. $P(X^{(n+1)}=y|X^{(j)}=x^{(j)},j=0,1,..,n)$. 当前阶段的结果仅依靠上一阶段的结果, 和之前的所有结果无关. 标注P为transition kernel, $S={x_{1},x_{2},…,x_{n}}$ 为state space.

我们可以将$P$表示为一个矩阵, 每个元素就是从某状态到下一状态的概率. 那么n轮之后结果的概率分布就是$P^n$. 对于一个固定的kernel, 结果的分布会在n足够大的时候趋于稳定. 而kernel会变化的时候不好说.

Markov chain的一系列性质就不赘述了. 对于我这种马工来说用不上嘻嘻.

Gibbs Sampler

暂时跳过. 反正可以直接调用bayesm包中的函数秒杀…

MCMC (马尔可夫链蒙特卡洛模拟法)

通过构造一个马尔可夫链, 使其平稳分布恰好为目标分布. 要求该kernel满足:

1. Detailed Balance: $\pi(x)P(x\to y)=\pi(y)P(y\to x)$.
2. Ergoicity: 从$x_{i}, \forall i$ 出发, $\Pr(x_{i}\to x_{j},\forall j)>0$.

利用该马尔可夫链生成样本${X_{1},..,X_{i}}$. 当样本生成的足够多了之后, 他自然而然会服从目标分布. 然后利用这些生成的样本就可以估计目标参数.

Metropdis-Hasting

从$X_{t}=x$出发, 据$q(x\to y)$ (提议分布)来生成y, 再用接受分布$\alpha(x\to y)$来决定是否保留y. 如果接受,则$X_{t+1}=y$, 否则$X_{t+1}=x$. 则总转移概率为: $P(x\to y)=q(x\to y)\alpha(x\to y)$, for $\alpha(x\to y)=\min (1, \frac{\pi(y)q(y\to x)}{\pi (x)q(x\to y)})$.

因为$x^{(i+1)}=y\cdot 1[u\le \alpha(x^{(i)},y)]+x^{(i)}1[u>\alpha(x^{(i)},y)]$. 有transition kernel: $P(x,y)=\left\{\begin{matrix}q(x^{(i)}|y)\alpha(x^{(i)},y),
&y\not = x \ 1-\sum_{v\not = x}q(x^{(i)}|v)\alpha(x^{(i)},y),
&y=x \end{matrix}\right.$

如果假设$\pi(x)q(y|x)>\pi(y)q(x|y)\leftrightarrow \alpha(x,y)<1\,\, \&\,\,\alpha(y,x)=1$. we have: $\pi(x)p(x,y)=\sim=\pi(y)p(y,x)\to $ Detailed Balance.

Variants

1. Random Walk Sampler
2. Independent Sampler

有空再介绍

Dealing with Missing Data (Data Augmentation)

我们可以使用这种方法来处理missing data: 记总共的data为: $y=\{y_{missing},y_{observed}\}$, 同时已知$f(y|\theta)$

那么有: $f(\theta|y_{obs})\propto f(y_{obs}|\theta)\pi(\theta)=\int f(y_{obs},y_{missing}|\theta)\pi(\theta)dy_{missing}\pi(\theta)$.

则通过Gibbs生成$y_{missing}: y_{missing}|y_{obs},\theta$. 可以用生成的data来得到后验分布: $\theta|y_{missing},y_{obs}$.

这相当于将$y_{missing}$当作未知变量: $\pi(\theta,y_{missing}|y_{obs})\propto f(y_{missing},y_{obs}|\theta)\pi(\theta)=f(y_{missing}|y_{obs},\theta)f(y_{obs}|\theta)\pi(\theta).$