笔记

贝叶斯数据分析(1)——入门

贝叶斯数据分析(1)——贝叶斯推断入门

本笔记基于Bayesian Data Analysis Third Edition, Andrew Gelman et. al. 学习编写,由于是英文教材,可能在学习过程中一些翻译或者内容有误,如有问题或错误,可以发送至邮箱[email protected]反馈。在学习本书中,需要有一定数学基础或者数理统计的基础,并且存在很多需要计算的场景,对于每一条定理,从头开始的证明会让你更明白每一步是如何实现的。

贝叶斯推断

概率标记

首先我们需要一些在贝叶斯统计里常用到的标记方法:

  1. p()p(\cdot|\cdot) 表示条件概率密度 (conditional probability density),p()p(\cdot)表示边缘分布 (marginal distribution). 应注意,我们通常使用的‘分布’ 和 ‘密度’ 两个概念通常是互相可以交换的。
  2. 为了混淆,我们会使用 Pr()Pr(\cdot) 代表事件的概率
  3. 如果使用到参数分布,我们需要标注分布的类型。比如 θ\theta 服从均数为 μ\mu, 方差为 σ\sigma 的正态分布,那么我们写作 θN(μ,σ2)\theta \sim N(\mu,\sigma^2) 或者 p(θ)=N(θμ,σ2)p(\theta) = N(\theta| \mu, \sigma^2), 又或者是更为详细的 p(θμ,σ2)=N(θμ,σ2)p(\theta|\mu, \sigma^2) = N(\theta| \mu, \sigma^2)
  4. 我们使用 N(μ,σ2)N(\mu,\sigma^2) 表示随机变量, N(θμ,σ2)N(\theta|\mu, \sigma^2) 作为密度方程。
  5. 协方差: sd(theta)/E(θ)sd(theta)/E(\theta), 几何均数 exp(E[log(θ)])\exp(E[\log(\theta)]), 几何标准差: exp(sd[log(θ)])\exp(sd[\log(\theta)])

贝叶斯法则

这里,我们需要明确几个概念,也是贝叶斯里面最重要的概念,我们一步步来看:

为了简化方便,现在下面所描述的均为单参数模型

我们从联合分布开始,对于联合分布 p(θ,y)p(\theta,y) ,我们可以写成:

p(θ,y)=p(θ)p(yθ)p(\theta,y) = p(\theta)p(y|\theta)

这里的 p(θ)p(\theta)先验分布(prior distribution),而p(yθ)p(y|\theta)则为样本分布(sampling distribution/data distribution)

我们往往想知道的是 p(θy)p(\theta|y),那么我们就可以通过上面的贝叶斯公式生成我们想要的 后验分布(posterior distribution):

p(θy)=p(θ,y)p(y)=p(θ)p(yθ)p(y)p(\theta|y) = \frac{p(\theta,y)}{p(y)} = \frac{p(\theta)p(y|\theta)}{p(y)}

这里的 p(y)=θp(θ)p(yθ)p(y) = \sum_\theta p(\theta)p(y|\theta) 又或者是 p(y)=p(θ)p(yθ)dθp(y) = \int p(\theta)p(y|\theta)d\theta, 这里的 p(y)p(y),我们很容易知道它是一个不依赖于 θ\theta 的函数,只与y有关,所以我们可以认为其是一个常数。于是乎,就有了非标准化后验密度(unnormalized posterior density)

p(θy)p(θ)p(yθ)p(\theta|y) \propto p(\theta)p(y|\theta)

应注意的是,p(yθ)p(y|\theta)是一个与 θ\theta 有关的一个函数,而并非 yy 的函数。这个基本上贯穿了我们整个贝叶斯的学习。

预测

然后正常的,我们会对于一个不知道的观测值进行推断,这时候我们叫做预测推断(predictive inference)

对于一个未知,但可以观测的 yy 的分布是:

p(y)=p(y,θ)dθ=p(θ)p(yθ)dθp(y) = \int p(y,\theta) d\theta = \int p(\theta) p(y|\theta) d\theta

这个也通常是 yy 的边缘分布,但其还有一个信息更加丰富的名字叫做先验预测分布(prior predictive distribution)。这个名字我们需要好好解释一下,先验:是因为它是非条件性的,因为它不基于先前的观测值;预测:是以为它是可观测量的分布。

当我们知道观测 yy 的数据后,我们就可以对不知道的可观测的 y~\tilde{y} 进行预测,那么对于这个 y~\tilde{y} 的分布,我们叫做后验预测分布(posterior predictive distribution)。一样的,我们也需要进行解释:后验,是因为它是基于已经观测到的值的条件化量,而预测是因为它是对可观察到的 y~\tilde{y} 的预测:

p(y~y)=p(y~,θy)dθ=p(y~θ,y)p(θy)dθ=p(y~θ)p(θy)dθ\begin{aligned} p(\tilde y|y) &= \int p(\tilde y, \theta|y)d\theta \\ &=\int p(\tilde y|\theta,y)p(\theta|y) d\theta \\ &=\int p(\tilde y|\theta)p(\theta|y) d\theta \end{aligned}

对于第二个等号,后验预测分布式条件预测在 θ\theta 的后验分布上的均值,而最后一步式因为我们假设 y~yθ\tilde y \perp y|\theta

似然(likelihood)

当我们对于一个选择好的概率模型使用贝叶斯就意味着我们的 yy 只通过 p(yθ)p(y|\theta) 来影响我们的后验推断,当我们认为其是一个固定 yy 的,关于 θ\theta 的函数,其也被称作似然方程(likelihood function) ,在贝叶斯推断里,它也偶尔被称作 似然法则(likelihood principle). 其表示在给定观测数据的情况下,关于模型参数的所有合理推断,仅依赖于数据的似然函数

似然与比值比

在一个给定模型下的后验密度的评价 θ1\theta_1θ2\theta_2 的比叫做后验比值 (posterior odds):

p(θ1y)p(θ2y)=p(θ1)p(yθ1)/p(y)p(θ2)p(yθ2)/p(y)=p(θ1)p(yθ1)p(θ2)p(yθ2)\frac{p(\theta_1|y)}{p(\theta_2|y)} = \frac{p(\theta_1)p(y|\theta_1)/p(y)}{p(\theta_2)p(y|\theta_2)/p(y)} = \frac{p(\theta_1)p(y|\theta_1)}{p(\theta_2)p(y|\theta_2)}

换句话,后验比 = 先验比 × 似然比

概率作为不确定性的衡量

在贝叶斯统计里,概率不仅用来描述重复随机实验的结果,也被用作衡量未知状态的基本标尺。贝叶斯方法关心的是:在已有数据和背景知识下,我们对于某个尚未观测、不可直接观测,或者尚未确定的状态到底知道多少。

从数学上看,概率被定义在一组可能结果上,且满足非负性、互斥事件可加性,以及所有互斥可能结果的概率总和为 1。贝叶斯统计在这个概率演算的基础上进一步强调:对于任何未知量,我们当前知识状态都可以用一个概率分布来描述。因此,后验分布不只是一个计算结果,它表示的是在结合数据之后,我们对未知参数或未知状态的不确定性。

为什么概率是一个合理的方式去量化不确定性呢?

  1. 类推(Analogy):物理随机性会导致不确定性,所以用随机事件的语言来描述不确定性是自然的。

  2. 公理和规范方法:这一思路通常和决策理论有关。我们可以把统计推断放到带有收益和损失的决策问题中,然后要求决策者的偏好满足一些合理公理,例如排序、传递性等。在这些规范性条件下,不确定性应当用概率表示。这个理由说明概率有很强的逻辑吸引力,但它更多是规范性的论证,并不能单独证明所有实际问题都必须这样建模。

  3. 投注的一致性:我们也可以用愿意接受的投注价格来定义主观概率。对于事件 EE,如果你愿意支付 pp 元,换取事件 EE 发生时得到 1 元的回报,那么 pp 就可以被看作你赋予事件 EE 的概率。如果 EE 发生,你净赚 1p1-p;如果 EE 不发生,你损失 pp。所谓一致性要求是:你给所有事件分配的概率不能让别人通过一组投注安排从你这里获得确定收益。可以证明,如果概率分配满足这种一致性要求,那么它必须满足概率论的基本公理。

主观性和客观性

所有使用概率的统计方法在某种意义上都是主观的,因为它们都依赖于对现实世界的数学理想化。贝叶斯方法常常因为使用先验分布而被认为更主观,但事实上,似然函数的选择、模型形式的设定、哪些数据应当纳入分析,同样都需要科学判断。例如,线性回归模型中的线性假设、误差分布假设,也可能和某个先验假设一样值得被质疑。

当我们拥有大量可交换的重复观测时,就可以从数据中估计概率分布的某些特征,从而让分析更接近“客观”。如果整个实验可以重复多次,甚至先验分布中的参数也可以通过数据来估计。不过,无论数据有多少,模型形式、数据选择和模型检验方式仍然需要我们在研究的过程中进行判断。

常用的概率理论知识

链式法则:

p(u,v,w)=p(uv,w)p(vw)p(w)p(u,v,w) = p(u|v,w)p(v|w)p(w)

隐藏假设:

p(θ,yH)=p(θH)p(yθ,H)p(\theta,y|H) = p(\theta|H)p(y|\theta,H)

​ 其中, HH 代表了一系列的假设或者猜想用于定义这个模型,比如我们有时候会用的协变量或者解释变量。

均值和方差:

E(u)=up(u)du,var(u)=(uE(u))2p(u)duE(u) = \int up(u)du, \qquad var(u) = \int(u-E(u))^2p(u)du

对于向量参数 uu,期望的表达式一样,其协方差矩阵为:

var(u)=(uE(u))(uE(u))Tp(u)duvar(u) = \int (u-E(u))(u-E(u))^Tp(u)du

我们应该要注意的是,E(u)E(u)var(u)var(u) 本质上不是关于 uu 的函数,而是关于分布方程的函数,即 p(u)p(u)

条件分布的均值和方差

均值:

离散形式下:

E(u)=E(E(uv))E(u) = E(E(u|v))

内层的期望是在基于条件 vv 的情况下 uu 的均值,外面那层的期望横跨 vv 的均值。对于 u.vu.v 的联合分布,其期望可以写作

E(u)=up(u,v)dudv=up(uv)p(v)dudv=up(uv)dup(v)dv=E(uv)p(v)dvE(u) = \int \int up(u,v) dudv = \int \int u p(u|v) p(v) dudv = \int\int up(u|v) du p(v)dv = \int E(u|v) p(v) dv

对于方差:

var(u)=E(var(uv))+var(E(uv))=E(E(u2v)(E(uv))2)+E((E(uv))2)(E(E(uv)))2=E(u2)E((E(uv))2)+E((E(uv))2)(E(u))2=E(u2)(E(u))2=var(u)\begin{aligned} var(u) &= E(var(u|v)) +var(E(u|v)) \\ &= E(E(u^2|v)-(E(u|v))^2) + E((E(u|v))^2)-(E(E(u|v)))^2 \\ &= E(u^2) - E((E(u|v))^2) + E((E(u|v))^2) - (E(u))^2 \\ &= E(u^2) - (E(u))^2 = var(u) \end{aligned}

证明:

var(x)=E[(XE(X))2]=E[X22XE(X)+E(X)2]=E(X2)2E(X)E(E(X))+E(X)2=E(X2)2E(X)2+E(X)2=E(X2)E(X)2\begin{aligned} var(x) &= E[(X - E(X))^2] \\ &= E[X^2 - 2XE(X)+E(X)^2] \\ &= E(X^2) - 2E(X)E(E(X)) + E(X)^2 \\ &= E(X^2) - 2E(X)^2 + E(X)^2 \\ &= E(X^2) - E(X)^2 \end{aligned}

变量转换

变量转换要解决的问题是:如果我们已经知道 uu 的分布 pu(u)p_u(u),并且定义一个新变量

v=f(u)v=f(u)

那么 vv 的分布 pv(v)p_v(v) 应该怎么写?我们很自然想到可以通过变换的方式得到:

离散变量

如果 uu 是离散变量,且 ff 是一一映射:

pv(v)=pu(f1(v))p_v(v)=p_u(f^{-1}(v))

如果 ff 是多对一映射,则需要把所有能映射到同一个 vvuu 加总:

pv(v)=u:f(u)=vpu(u)p_v(v)=\sum_{u:f(u)=v}p_u(u)

比如 u{1,0,1}u\in\{-1,0,1\}v=u2v=u^2,则:

pv(0)=pu(0),pv(1)=pu(1)+pu(1)p_v(0)=p_u(0), \qquad p_v(1)=p_u(-1)+p_u(1)

连续变量

对于连续变量,我们不能像离散变量一样只做简单代入,因为密度会受到尺度变化的影响。设 v=f(u)v=f(u) 是一一变换,且反函数为

u=f1(v)=g(v)u=f^{-1}(v)=g(v)

则其可以通过如下方式变换:

pv(v)=Jpu(f1(v))p_v(v)=|J|p_u(f^{-1}(v))

其中 JJ 是反变换 u=f1(v)u=f^{-1}(v) 关于 vv 的 Jacobian 行列式。

在一维时:

pv(v)=pu(g(v))dg(v)dvp_v(v)=p_u(g(v))\left|\frac{dg(v)}{dv}\right|

或者写成正向变换的形式:

我会觉得这个形式更加常用和方便计算

pv(v)=pu(u)dudv=pu(u)dvdu,u=f1(v)p_v(v)=p_u(u)\left|\frac{du}{dv}\right| =\frac{p_u(u)}{\left|\frac{dv}{du}\right|}, \qquad u=f^{-1}(v)

我们可以这么想,因为密度不等于概率,而是在单位长度上的概率。变量变换以后,单位长度被拉伸或压缩,所以要乘上尺度校正项,更像是我们在做一个坐标轴变换的事情:

pu(u)du=pv(v)dvp_u(u)du=p_v(v)dv

因此:

pv(v)=pu(u)dudvp_v(v)=p_u(u)\left|\frac{du}{dv}\right|

Jacobian 矩阵

如果 u=(u1,,uk)Tu=(u_1,\dots,u_k)^Tv=(v1,,vk)Tv=(v_1,\dots,v_k)^T,且 u=g(v)=f1(v)u=g(v)=f^{-1}(v),则:

Jg(v)=uv=(u1v1u1vkukv1ukvk)J_g(v)= \frac{\partial u}{\partial v} = \begin{pmatrix} \frac{\partial u_1}{\partial v_1} & \cdots & \frac{\partial u_1}{\partial v_k}\\ \vdots & \ddots & \vdots\\ \frac{\partial u_k}{\partial v_1} & \cdots & \frac{\partial u_k}{\partial v_k} \end{pmatrix}

连续型向量变量的变换公式为:

pv(v)=pu(g(v))detJg(v)p_v(v)=p_u(g(v))\left|\det J_g(v)\right|

如果使用正向变换 v=f(u)v=f(u) 的 Jacobian:

Jf(u)=vuJ_f(u)=\frac{\partial v}{\partial u}

那么:

detJg(v)=1detJf(u),u=g(v)\left|\det J_g(v)\right| = \frac{1}{\left|\det J_f(u)\right|}, \qquad u=g(v)

所以等价地:

pv(v)=pu(u)detJf(u),u=f1(v)p_v(v) = \frac{p_u(u)}{\left|\det J_f(u)\right|}, \qquad u=f^{-1}(v)

多对一的连续变换同样要把所有反函数分支加总:

pv(v)=uj:f(uj)=vpu(uj)dujdvp_v(v)=\sum_{u_j:f(u_j)=v}p_u(u_j) \left|\frac{du_j}{dv}\right|

一个常见例子是 v=u2v=u^2。若 uu 是连续变量,则当 v>0v>0 时有两个反函数分支:

u1=v,u2=vu_1=\sqrt v,\qquad u_2=-\sqrt v

并且:

du1dv=du2dv=12v\left|\frac{du_1}{dv}\right| = \left|\frac{du_2}{dv}\right| = \frac{1}{2\sqrt v}

因此:

pv(v)=pu(v)12v+pu(v)12v,v>0p_v(v) = p_u(\sqrt v)\frac{1}{2\sqrt v} + p_u(-\sqrt v)\frac{1}{2\sqrt v}, \qquad v>0

常用的参数变换

对数转换

对于正数参数 u(0,)u\in(0,\infty),常用对数变换:

v=logu,u=evv=\log u,\qquad u=e^v

Jacobian矩阵为:

dudv=ev\left|\frac{du}{dv}\right|=e^v

所以:

pv(v)=pu(ev)evp_v(v)=p_u(e^v)e^v

反过来,如果已知 vv 的密度,u=evu=e^v,则:

pu(u)=pv(logu)1up_u(u)=p_v(\log u)\frac{1}{u}
logit变换

对于概率参数 u(0,1)u\in(0,1),常用 logit 变换:

logit(u)=log(u1u)\operatorname{logit}(u)=\log\left(\frac{u}{1-u}\right)

其反变换为:

logit1(v)=ev1+ev\operatorname{logit}^{-1}(v)=\frac{e^v}{1+e^v}

这时:

u=ev1+ev,dudv=u(1u)=ev(1+ev)2u=\frac{e^v}{1+e^v},\qquad \frac{du}{dv}=u(1-u)=\frac{e^v}{(1+e^v)^2}

因此如果 uu 的密度为 pu(u)p_u(u),则:

pv(v)=pu(ev1+ev)ev(1+ev)2p_v(v)=p_u\left(\frac{e^v}{1+e^v}\right) \frac{e^v}{(1+e^v)^2}
probit变换

另一个常见的变换是 probit 变换:

v=Φ1(u),u=Φ(v)v=\Phi^{-1}(u),\qquad u=\Phi(v)

其中 Φ\Phi 是标准正态分布的 CDF。因为:

dudv=ϕ(v)\frac{du}{dv}=\phi(v)

所以:

pv(v)=pu(Φ(v))ϕ(v)p_v(v)=p_u(\Phi(v))\phi(v)

逆累积分布函数抽样(Inverse CDF)

模拟的基本想法是:如果能从 UUniform(0,1)U\sim \operatorname{Uniform}(0,1) 生成随机数,就可以通过分布函数的反函数生成其他分布的随机数。

一维随机变量 vv 的累积分布函数为:

F(v)=Pr(vv)={vvp(v),p 是离散分布vp(v)dv,p 是连续分布F(v^*)=\Pr(v\le v^*)= \begin{cases} \sum_{v\le v^*}p(v), & p \text{ 是离散分布}\\ \int_{-\infty}^{v^*}p(v)dv, & p \text{ 是连续分布} \end{cases}

逆 CDF 抽样的步骤:

  1. 生成
UUniform(0,1)U\sim \operatorname{Uniform}(0,1)
V=F1(U)V=F^{-1}(U)

其中 VV 服从目标分布 p(v)p(v)

更严格地说,尤其是在离散分布里,FF 可能不是一一映射,因此常用广义逆函数:

F1(u)=inf{v:F(v)u}F^{-1}(u)=\inf\{v:F(v)\ge u\}

连续情形下可以用一行证明:

Pr(F1(U)v)=Pr(UF(v))=F(v)\Pr(F^{-1}(U)\le v)= \Pr(U\le F(v))=F(v)

所以 F1(U)F^{-1}(U) 的 CDF 正好是 FF

离散分布的抽样

设离散变量可能取值为 v1,,vKv_1,\dots,v_K,概率为:

Pr(V=vk)=qk,k=1Kqk=1\Pr(V=v_k)=q_k,\qquad \sum_{k=1}^Kq_k=1

定义累积概率:

ck=j=1kqj,c0=0c_k=\sum_{j=1}^kq_j,\qquad c_0=0

若:

ck1<Uckc_{k-1}<U\le c_k

则取:

V=vkV=v_k

即离散分布版本的 V=F1(U)V=F^{-1}(U)

指数分布

若:

VExponential(λ)V\sim \operatorname{Exponential}(\lambda)

则其 CDF 为:

F(v)=1eλv,v>0F(v)=1-e^{-\lambda v},\qquad v>0

U=F(v)U=F(v)

U=1eλvU=1-e^{-\lambda v}

解得:

v=log(1U)λv=-\frac{\log(1-U)}{\lambda}

由于 1U1-U 仍然服从 Uniform(0,1)\operatorname{Uniform}(0,1),所以可直接写为:

V=logUλV=-\frac{\log U}{\lambda}

正态分布

标准正态分布没有简单的初等反函数,但软件可以计算 Φ1\Phi^{-1},则:

Z=Φ1(U)N(0,1)Z=\Phi^{-1}(U)\sim N(0,1)

进一步:

X=μ+σΦ1(U)N(μ,σ2)X=\mu+\sigma \Phi^{-1}(U)\sim N(\mu,\sigma^2)

分位数和模拟样本

从分布 p(ω)p(\omega) 中生成 SS 个样本:

ω(1),,ω(S)\omega^{(1)},\dots,\omega^{(S)}

可以用样本经验分布近似真实分布。比如估计 95% 分位数时,将样本排序:

ω(1:S)ω(2:S)ω(S:S)\omega^{(1:S)}\le \omega^{(2:S)}\le \cdots \le \omega^{(S:S)}

取近似:

q0.95ω(0.95S:S)q_{0.95}\approx \omega^{(\lceil0.95S\rceil:S)}

对于大部分的情况,S=1000S=1000 通常已经足够得到比较稳定的 95% 分位数估计