统计

发表于 阅读次数： Waline： 本文字数： 2.8k 阅读时长 ≈ 10 分钟

All models are wrong, but some are useful. – George E.P. Box

  • 概率论：在良定义/理想化的模型中、基于严格的理论、考虑一个事件的可能性
  • 统计学：观察、收集、去芜存菁、总结、抽象、假设、建模、预测、推断

点估计

最大似然估计

最大似然估计（Maximum Likelihood Estimation, MLE）：给定一个模型，找到最有可能的参数值。

投一枚硬币 100 次，49 次正面，51 次反面。硬币正面的概率是多少？

正面概率是未知参数 pp，概率质量函数

f(x;p)={p,x=11−p,x=0f(x; p) = \begin{cases} p, & x = 1 \\ 1 - p, & x = 0 \end{cases}

似然函数

L(p)=L(p;x1,…,xn)=∏i=1nf(xi;p)=p49(1−p)51L(p) = L(p; x_1, \dots, x_n) = \prod_{i=1}^n f(x_i; p) = p^{49} (1 - p)^{51}

从而有

L′(p)=(1−p)50p48(49−100p)=0L'(p) = (1-p)^{50} p^{48} (49 - 100p) = 0

得出 p=0.49p = 0.49。

估计量 p^\hat{p} 的最大似然估计值 p^=0.49\hat{p} = 0.49。即哪个参数最有可能使得当前数据出现：arg max⁡p^Pr⁡(X=x∣p=p^)\argmax\limits_{\hat{p}} \Pr(X = x \mid p = \hat{p})。

最大似然估计期望并不总是等于样本均值：

  • 是：伯努利试验、二项分布、泊松分布、正态分布
  • 否：几何分布、指数分布

求最大似然估计即最优化似然函数 L(p)L(p)。

这不总是易解，还有迭代逼近的方法：

  • 梯度下降法
  • 牛顿-拉弗森法
  • 拟牛顿法
  • 期望最大化（Expectation Maximization, EM）算法

矩方法

矩问题：两个概率分布如果所有的矩都相同，那么它们是同一个分布。

根据各阶矩可以确定一个概率分布，即给定参数之后的那个未知分布。

大数定理有

Xˉnk=1n(X1k+⋯+Xnk)→E[Xk]\bar{X}_n^{k} = \dfrac{1}{n}(X_1^{k} + \dots + X_n^{k}) \to \mathbb{E}[X^{k}]

kk-阶样本矩依概率/几乎处处收敛到 kk-阶矩。

假设样本服从某正态分布 N(μ,σ2)\mathcal{N}(\mu, \sigma^2)。

有 E[X]=μ\mathbb{E}[X] = \mu 与 E[X2]=μ2+σ2\mathbb{E}[X^2] = \mu^2 + \sigma^2。

样本矩有

{m1=1n∑i=1nXim2=1n∑i=1nXi2\left\lbrace\begin{aligned} m_1 = \dfrac{1}{n} \sum_{i=1}^n X_i \\ m_2 = \dfrac{1}{n} \sum_{i=1}^n X_i^2 \end{aligned}\right.

令 m1=μ,m2=μ2+σ2m_1 = \mu, m_2 = \mu^2 + \sigma^2 解得

{μ^=m1=Xˉnσ^2=m2−m12=1n∑i=1n(Xi−Xˉn)2=Sn2\left\lbrace\begin{aligned} \hat{\mu} &= m_1 = \bar{X}_n\\ \hat{\sigma}^2 &= m_2 - m_1^2 = \dfrac{1}{n} \sum_{i=1}^n (X_i - \bar{X}_n)^2 = S_n^2 \end{aligned}\right.

假设样本服从某泊松分布 Pois⁡(λ)\operatorname{Pois}(\lambda)。

有 E[X]=λ\mathbb{E}[X] = \lambda 与 E[X2]=λ2+λ\mathbb{E}[X^2] = \lambda^2 + \lambda。

同样解得

{λ^1=Xˉnλˉ2=Sn2−Xˉn−14\left\lbrace\begin{aligned} \hat{\lambda}_1 &= \bar{X}_n\\ \bar{\lambda}_2 &= \sqrt{S_n^2 - \bar{X}_n} - \dfrac{1}{4} \end{aligned}\right.

说明矩估计答案并不唯一。一般情况下采用低阶矩的结果，即 λ^=λ^1=Xˉn\hat{\lambda} = \hat{\lambda}_1 = \bar{X}_n。

贝叶斯估计

对于投针试验，MLE 结果是 3.13.1，在已知先验知识 π≈3.14\pi \approx 3.14 的情况下：

  • 频率学派：限制可选的 π∈[3.14,3.15]\pi \in [3.14, 3.15]，MLE 修正为 π=3.14\pi = 3.14；
  • 贝叶斯学派：假设 π\pi 服从与 3.143.14 有关的先验概率分布 Pr⁡(π=π^)\Pr(\pi = \hat{\pi})。
    • 通过后验概率分布 Pr⁡(π=π^∣X=x)\Pr(\pi = \hat{\pi} \mid X = x) 选择最好的估计值 π^\hat{\pi}；
    • 给定数据 xx，已知 Pr⁡(X=x∣π=π^)\Pr(X = x \mid \pi = \hat{\pi}) 和 Pr⁡(X=x)\Pr(X = x)；
    • 贝叶斯定理有 Pr⁡(π=π^∣X=x)=Pr⁡(X=x∣π=π^)Pr⁡(π=π^)Pr⁡(X=x)\Pr(\pi = \hat{\pi} \mid X = x) = \dfrac{\Pr(X = x \mid \pi = \hat{\pi}) \Pr(\pi = \hat{\pi})}{\Pr(X = x)}。

选择先验分布：

  • 应反映现实情况
    • 实践中选择比较任意
    • 选择适应面较广的分布，通过调参具体确定
  • 便于计算

共轭（conjugacy）分布*：P(π)P(\pi) 与 P(π∣x)P(\pi \mid x) 同分布族，P(π)P(\pi) 是 P(x∣π)P(x \mid \pi) 的共轭先验。

由于 P(π∣x)∝P(x∣π)⋅P(π)P(\pi \mid x) \propto P(x \mid \pi) \cdot P(\pi)，选择 π^\hat{\pi} 只需比较 P(x∣π)⋅P(π)P(x \mid \pi) \cdot P(\pi)。

计算后验预测分布 P(x~∣x)=∫πP(x~∣π)⋅P(π∣x) ⁣dπP(\tilde{x} \mid x) = \displaystyle \int_{\pi} P(\tilde{x} \mid \pi) \cdot P(\pi \mid x) \d \pi。

假设样本服从某伯努利分布，参数 p=θp = \theta，令 Y=∑iXi∼Bin⁡(n,p)Y = \sum_i X_i \sim \operatorname{Bin}(n, p)。

则伯努利分布与二项分布的共轭先验分布是贝塔分布。

贝塔分布

Beta⁡(x,y)=∫01tx−1(1−t)y−1 ⁣dt∼2πxx−1/2yy−1/2(x+y)x+y−1/2\operatorname{Beta}(x, y) = \int_0^1 t^{x-1} (1-t)^{y-1} \d t \sim \sqrt{2 \pi} \dfrac{x^{x- 1 / 2}y^{y - 1 / 2}}{(x + y)^{x + y - 1 / 2}}

假设后验分布 P(θ)=θα−1(1−θ)β−1Beta⁡(α,β)P(\theta) = \dfrac{\theta^{\alpha-1}(1-\theta)^{\beta-1}}{\operatorname{Beta}(\alpha, \beta)}，α,β\alpha, \beta 是灵活调节的参数。

有

P(y)=∫01P(y∣θ)P(θ) ⁣dθ=(ny)Beta⁡(α+y,n+β−y)Beta⁡(α,β)\begin{aligned} P(y) &= \int_0^1 P(y \mid \theta) P(\theta) \d \theta \\ &= \dbinom{n}{y} \dfrac{\operatorname{Beta}(\alpha + y, n + \beta - y)}{\operatorname{Beta}(\alpha, \beta)} \end{aligned}

于是

P(θ∣y)=P(y,θ)P(y)=θy+α−1(1−θ)n−y+β−1Beta⁡(α+y,n+β−y)\begin{aligned} P(\theta \mid y) &= \dfrac{P(y, \theta)}{P(y)}\\ &= \dfrac{\theta^{y + \alpha - 1} (1 - \theta)^{n - y + \beta - 1}}{\operatorname{Beta}(\alpha + y, n + \beta - y)} \end{aligned}

从而有

P(θ∣y)∝P(y∣θ)P(θ)=θy+α−1(1−θ)n−y+β−1Beta⁡(α,β)\begin{aligned} P(\theta \mid y) &\propto P(y \mid \theta) P(\theta)\\ &= \dfrac{\theta^{y + \alpha - 1} (1 - \theta)^{n - y + \beta - 1}}{\operatorname{Beta}(\alpha, \beta)} \end{aligned}

给定当前数据，最有可能的 π^\hat{\pi}：arg max⁡π^P(π=π^∣X=x)\argmax\limits_{\hat{\pi}} P(\pi = \hat{\pi} \mid X = x)。

有最大后验估计（Maximum A Posteriori, MAP）：

θ^=y+α−1n+α+β−2\hat{\theta} = \dfrac{y + \alpha - 1}{n + \alpha + \beta - 2}

后验预测分布

P(y~∣y)=∫01P(y~∣θ)P(θ∣y) ⁣dθ=y+αn+α+β\begin{aligned} P(\tilde{y} \mid y) &= \int_0^1 P(\tilde{y} \mid \theta) P(\theta \mid y) \d \theta\\ &= \dfrac{y + \alpha}{n + \alpha + \beta} \end{aligned}

比较估计量

  • 最大似然估计（MLE）
    • 哪个参数最有可能使得当前数据出现
    • 直观，适合大样本，有时不易计算
  • 矩估计（MOM）
    • 大数定律：样本矩会收敛到总体矩
    • 容易计算，有时不太准确，甚至自相矛盾
  • 最大后验估计（MAP）
    • 给定当前数据，最有可能的那个参数
    • 引入先验知识的辅助，适合小样本，经常难以计算

相合/一致（Consistency）

大数定理保证了 Xˉn→E[X]\bar{X}_n \to \mathbb{E}[X]。对于估计量能否有数据越多，估计越准？

相合（Consistency）：θ^n→θ\hat{\theta}_n \to \theta

  • （弱）相合：θ^→Pθ\hat{\theta} \xrightarrow[]{P} \theta
  • 强相合*：θ^→a.s.θ\hat{\theta} \xrightarrow[]{\text{a.s.}} \theta
  • rr 阶矩相合*：∣θ^−θ^∣r→0|\hat{\theta} - \hat{\theta}|^r \to 0

无偏（Unbiasedness）

假设样本 X1,…,XnX_1, \dots, X_n 服从某分布，估计其方差。

二阶样本中心矩 Sn2=∑i(Xi−Xˉ)2nS_n^2 = \displaystyle \sum_i \dfrac{(X_i - \bar{X})^2}{n}，大数定理有 Sn2→E[(X−μ)2]S_n^2 \to \mathbb{E}[(X - \mu)^2]。

有

E[Sn2]=∑iE[Xi2]−2E[XˉnXi]+E[Xˉn2]n=E[X2]−E[Xˉn2]=E[X2]−E[∑iXi2]+∑i≠jE[XiXj]n2=n−1n(E[X2]−E[X]2)=n−1nσ2\begin{aligned} \mathbb{E}[S_n^2] &= \sum_i \dfrac{\mathbb{E}[X_i^2] - 2\mathbb{E}[\bar{X}_n X_i] + \mathbb{E}[\bar{X}_n^2]}{n}\\ &= \mathbb{E}[X^2] - \mathbb{E}[\bar{X}_n^2]\\ &= \mathbb{E}[X^2] - \dfrac{\mathbb{E}[\sum_i X_i^2] + \sum_{i \ne j}\mathbb{E}[X_i X_{j}]}{n^2}\\ &= \dfrac{n-1}{n}(\mathbb{E}[X^2] - \mathbb{E}[X]^2)\\ &= \dfrac{n-1}{n}\sigma^2 \end{aligned}

第二个等号是因为

E[XˉnXi]=1n∑jE[XiXj]=1n2∑i,jE[XiXj]=E[∑iXin∑jXjn]=E[Xˉn2]\begin{aligned} \mathbb{E}[\bar{X}_n X_i] &= \dfrac{1}{n} \sum_{j} \mathbb{E}[X_i X_j]\\ &= \dfrac{1}{n^2} \sum_{i, j} \mathbb{E}[X_i X_{j}]\\ &= \mathbb{E}\left[ \sum_i \dfrac{X_i}{n} \sum_{j} \dfrac{X_{j}}{n} \right] \\ &= \mathbb{E}[\bar{X}_n^2] \end{aligned}

Sn2S_n^2 称为渐近无偏估计量（未修正的样本方差）。有 E[Sn2]→σ2\mathbb{E}[S_n^2] \to \sigma^2。

样本方差的无偏估计量

S2=nn−1Sn2=∑i(Xi−Xˉ)2n−1S^2 = \dfrac{n}{n-1} S_n^2 = \dfrac{\sum_i (X_i - \bar{X})^2}{n-1}

有 E[S2]=σ2\mathbb{E}[S^2] = \sigma^2。

偏差（bias）为 E[θ^]−θ\mathbb{E}[\hat{\theta}] - \theta。

平均绝对误差（average absolute deviation）为 E[∣θ^−θ∣]\mathbb{E}[|\hat{\theta} - \theta|]。

均方误差（mean squared error, MSE）为 E[(θ^−θ)2]\mathbb{E}[(\hat{\theta} - \theta)^2]。

贝叶斯估计：

  • 最大后验估计（MAP）：π^=arg max⁡π^P(π=π^∣X=x)\hat{\pi} = \argmax\limits_{\hat{\pi}} P(\pi = \hat{\pi} \mid X = x)
  • 后验中位数估计（Posterior median）：平均绝对误差最小 π^=arg min⁡π^E[∣π−π^∣∣X=x]\hat{\pi} = \argmin\limits_{\hat{\pi}} \mathbb{E}[|\pi - \hat{\pi}| \mid X = x]；
  • 最小均方估计（Least Mean Squares, LMS）：最小均方误差估计（MMSE）π^=arg min⁡π^E[(π−π^)2∣X=x]\hat{\pi} = \argmin\limits_{\hat{\pi}} \mathbb{E}[(\pi - \hat{\pi})^2 \mid X = x]。

有效性（Efficiency）

估计量有方差，方差更小的无偏估计更「有效」（efficient）。

最小方差无偏估计（Minimum Variance Unbiased Estimator, MVUE）：方差最小的无偏估计。

若无论 θ\theta 的取值，θ^(X;θ)\hat{\theta}(X; \theta) 都是 MVUE，则它是一致最小无偏估计（Uniformly Minimum Variance Unbiased Estimator, UMVUE）。UMVUE 是唯一的。

Cramér–Rao bound* & Fisher information*

Cramér-Ra bound

假设样本 X1,…,XnX_1, \dots, X_n 服从某分布，pdf 为 f(x;θ)f(x; \theta)。若 θ^(X1,…,Xn)\hat{\theta}(X_1, \dots, X_n) 是 θ\theta 的无偏估计量，则

Var⁡[θ^]⩾1n⋅I(θ)\operatorname{Var} [\hat{\theta}] \ge \dfrac{1}{n \cdot I(\theta)}

其中 I(θ)I(\theta) 是该概率分布关于未知参数 θ\theta 的费希尔讯息数（Fisher Information）：

I(θ)=E[(∂∂θlog⁡f(X;θ))2]=∫−∞∞(∂∂θlog⁡f(x;θ))2f(x;θ) ⁣dx\begin{aligned} I(\theta) &= \mathbb{E}\left[ \left( \dfrac{\partial}{\partial \theta} \log f(X; \theta) \right)^2 \right]\\ &= \int_{-\infty}^{\infty} \left( \dfrac{\partial}{\partial \theta} \log f(x; \theta) \right)^2 f(x; \theta) \d x \end{aligned}

Cramér-Rao bound 不一定可达。

若 θ^\hat{\theta} 是 θ\theta 的无偏估计：

  • 有效估计（efficient estimator）：Var⁡[θ^]=1nI(θ)\operatorname{Var} [\hat{\theta}] = \dfrac{1}{n I(\theta)}，即效率为 11；
  • θ^\hat{\theta} 的效率（efficiency）为 en(θ^)=1nI(θ)Var⁡[θ^]∈[0,1]e_n(\hat{\theta}) = \dfrac{1}{n I(\theta) \operatorname{Var}[\hat{\theta}]} \in [0, 1]；
  • 渐近有效（asymptotically efficient）：en(θ^)→1e_n(\hat{\theta}) \to 1。

充分统计量（sufficient statistic）*

假设样本 X1,…,XnX_1, \dots, X_n 服从某伯努利分布 p=θp = \theta，令 Y=∑iXi∼Bin⁡(n,θ)Y = \sum_i X_i \sim \operatorname{Bin}(n, \theta)。

已知 YY，在获得别的信息，能更好地估计 θ\theta 吗？

若给定统计量 T(X1,…,Xn)T(X_1, \dots, X_n)，样本 (X1,…,Xn)(X_1, \dots, X_n) 与参数 θ\theta 无关，则统计量 TT 是参数 θ\theta 的充分统计量。

充分统计量 TT 包含了样本中所有关于参数 θ\theta 的信息。

Fisher-Neyman 因子分解定理

令 ff 是似然函数/pdf/pmf，统计量 TT 是参数 θ\theta 的充分统计量当且仅当存在非负函数 h,gh, g 使得

f(x;θ)=h(x)⋅g(θ,T(x))f(x; \theta) = h(x) \cdot g(\theta, T(x))

概率分布可以分解为两个函数的乘积，其中一个与 θ\theta 无关，另一个仅通过 TT 与样本 xx 产生关联。

充分统计量不唯一。如果统计量 S(X)S(X) 是 θ\theta 的充分统计量，且对于任意 θ\theta 的充分统计量 T(X)T(X) 都存在函数 gg 使得 S(X)=g(T(X))S(X) = g(T(X))，则 S(X)S(X) 是 θ\theta 的最小充分统计量（minimal sufficiency）。

最小充分统计量通常存在，但不总是存在。

Rao-Blackwell-Kolmogorov 定理

已知参数 θ\theta 的充分统计量 TT，对于任意估计量 θ^1(X)\hat{\theta}_1(X)，利用充分统计量的新估计量 θ^2=E[θ^1(X)∣T(X)]\hat{\theta}_2 = \mathbb{E}[ \hat{\theta}_1(X) \mid T(X)]，有

E[(θ^2(X)−θ)2]⩽E[(θ^1(X)−θ)2]\mathbb{E}\left[ (\hat{\theta}_2(X) - \theta)^2 \right] \le \mathbb{E}\left[ (\hat{\theta}_1(X) - \theta)^2 \right]

利用充分统计量可以得到均方误差更小的估计量。

Lehmann-Scheffé 定理

若 θ^\hat{\theta} 是 θ\theta 的无偏估计量，TT 是完备（complete）充分统计量，则 E[θ^∣T]\mathbb{E}[\hat{\theta} \mid T] 是唯一的一致最小方差无偏估计（UMVUE）。

不严格地说，完备统计量只包含了样本中关于目标参数的信息，不含其他信息。

区间估计（interval estimation）

假设样本 X1,…,XnX_1, \dots, X_n 服从某伯努利分布 p=θp = \theta，令 Y=∑iXi∼Bin⁡(n,θ)Y = \sum_i X_i \sim \operatorname{Bin}(n, \theta)。

霍夫丁不等式有

Pr⁡(∣Y−np∣⩾t)⩽2exp⁡(−2t2n)\Pr(|Y - np| \ge t) \le 2 \exp\left( -\dfrac{2t^2}{n} \right)

取一定的 tt（Θ(n)\Theta(\sqrt{n})），可使 Pr⁡[p=Yn±t]=0.95\Pr\left[p = \dfrac{Y}{n} \pm t\right] = 0.95，即以 95%95\% 的概率，p=Y/n±O(n)p = Y/n \pm O(\sqrt{n})。

更严格地说，应该是有 95%95\% 的概率，该置信区间包含真实参数 pp。而非 pp 在该区间的概率为 95%95\%，因为 pp 是一个固定值。

Yn±O(n)\dfrac{Y}{n} \pm O(\sqrt{n}) 是一个随机的区间，取决于随机样本 X1,…,XnX_1, \dots, X_n。

从后验分布中选择一段区间（可信区间，credible interval）使得

Pr⁡[θ∈[a,b]∣X=x]=1−α\Pr[\theta \in [a, b] \mid X = x] = 1 - \alpha

选法*：

  • 选择最短的区间（highest density interval, HDI）
  • 基于分位数（quantile-based interval, QBI）：如中位数，左右各取 1−α2\dfrac{1-\alpha}{2}。

（原则上）给定置信度 1−α1 - \alpha，任意统计量 a(X1,… ), b(X1,… )a(X_1, \dots),\, b(X_1, \dots)，使得

Pr⁡[a⩽θ⩽b]⩾1−α\Pr[a \le \theta \le b] \ge 1 - \alpha

枢轴变量（pivot/pivotal quantity）法：假设样本 X1,…,XnX_1, \dots, X_n 服从 N(μ,1)\mathcal{N}(\mu, 1)：

  • 统计量 Xˉn∼N(μ,1n)\bar{X}_n \sim \mathcal{N}\left(\mu, \frac{1}{n}\right)，枢轴变量 Xˉn−μ∼N(0,1n)\bar{X}_n - \mu \sim \mathcal{N}\left(0, \frac{1}{n}\right)
  • 寻找分位数 [c,d][c, d] 使得 Pr⁡[Xˉn−μ∈[c,d]]⩾1−α\Pr[\bar{X}_n - \mu \in [c, d]] \ge 1 - \alpha
  • 置信水平 1−α1 - \alpha 的置信区间为 [Xˉn−d,Xˉn−c][\bar{X}_n - d, \bar{X}_n - c]