贝叶斯估计

贝叶斯估计通过先验分布与数据似然得到后验分布，再根据指定的损失函数选择估计或行动。后验分布描述参数的不确定性；从中取均值、众数或分位点，是另一个需要明确目标的步骤。不同损失可以从同一后验得到不同的最优决定。

例如两次合成测量为2和4，测量噪声方差已知为四。若未知总体均值μ事先服从 $N(0,1)$，则后验均值为一，后验方差为2/3。样本均值是三，结果向先验均值零移动，但并非把先验和两次数值简单平均；各部分的精度决定权重。

从后验风险到估计规则

设参数先验为π，观测的条件密度或质量为 $f(y\mid \theta )$，证据m(y)有限且正，从而后验 $\pi (\theta \mid y)=f(y\mid \theta )\pi (\theta )/m(y)$ 有定义。给定行动空间A与损失 $L(\theta ,a)\ge 0$，后验风险是 $$\rho (a\mid y)=𝔼[L(\theta ,a)\mid y].$$ 贝叶斯行动选取其最小点 $a^{\ast }(y)\in \mathrm{\ast }{argmin}_{a\in A}\rho (a\mid y)$，若最小点可测且能取得。理由是将先验平均风险按数据分组，利用非负积分换序，得到 $${𝔼}_{\theta ,Y}L(\theta ,a(Y))=\int \rho (a(y)\mid y)m(y)\mathrm{d}y.$$ 对每个数据y分别使被积的后验风险最小，就使整体先验平均风险最小。这个平均是在给定先验下作的，不表示对每一个固定真实参数都支配其他估计。

若 $A=ℝ$、损失为平方误差、后验二阶矩有限，记 $m=𝔼(\theta \mid y)$，则展开平方得 $$\rho (a\mid y)=𝔼[(\theta -m+m-a{)}^2\mid y]=\mathrm{Var}(\theta \mid y)+(a-m{)}^2.$$ 交叉项为零，因此唯一最优行动是后验均值m。若损失改为绝对误差，并且后验满足 $𝔼(|\theta |\mid y)<\mathrm{\infty }$，后验中位数是最优；对有连续分布函数G的情形，有限风险的导数为 $2G(a)-1$，其符号在中位数两侧改变。若行动只允许整数或其他约束集合，不能直接照搬无约束均值，须在允许行动中最小化风险。

正态均值的归一化后验

假定给定μ后 $Y_1,\dots ,Y_n$ 独立同分布于 $N(\mu ,{\sigma }^2)$，σ的平方为已知正数；先验为 $\mu \sim N(m_0,v_0)$，v0为正方差。似然乘先验的指数中，与μ有关的二次式为 $$\frac{1}{2{\sigma }^2}\sum _i(y_i-\mu {)}^2+\frac{(\mu -m_0{)}^2}{2v_0}=\frac{{\mu }^2}{2}(\frac{n}{{\sigma }^2}+\frac{1}{v_0})-\mu (\frac{\sum _i{y}_i}{{\sigma }^2}+\frac{m_0}{v_0})+\text{常数}.$$ 令 $$v_n={(v_0^{-1}+n{\sigma }^{-2})}^{-1},m_n=v_n(\frac{m_0}{v_0}+\frac{\sum _i{y}_i}{{\sigma }^2}),$$ 配方后得到 $(\mu -m_n{)}^2/(2v_n)$ 加与μ无关的项。积分高斯核给出归一化密度 $(2\pi v_n{)}^{-1/2}\exp [-(\mu -m_n{)}^2/(2v_n)]$，所以后验恰为 $N(m_n,v_n)$。方差没有未知常数，平方损失的最优估计为mn。

对于n大于零，mn也可写成两个均值按精度的加权：先验精度为 $1/v_0$，样本均值精度为 $n/{\sigma }^2$。更精确的信息权重更大。这里σ的平方是单次观测方差，不能误写为样本均值方差，也不能把v0当成标准差。

主例 $m_0=0,v_0=1,{\sigma }^2=4,n=2,\sum y_i=6$，故后验精度 $1+2/4=3/2$，$v_n=2/3$、$m_n=(2/3)(6/4)=1$。平方损失风险在行动0、1、2处依次为 $5/3,2/3,5/3$；在样本均值三处为14/3。后验均值一确实使当前后验平方风险最小，这不声称当固定真实μ很远离零时该先验一定有利。

后验比先验更集中；曲线归一化为一，横轴是未知均值μ。

图中行动0、1、2的风险分别为5/3、2/3、5/3，最低值仍为正，反映参数剩余不确定性。

预测仍要保留参数的不确定性

新观测与已有观测条件于μ相互独立，且 $Y_{\mathrm{n}\mathrm{e}\mathrm{w}}\mid \mu \sim N(\mu ,4)$。积分掉后验μ，可以表示为 $Y_{\mathrm{n}\mathrm{e}\mathrm{w}}=m_n+\sqrt{v_n}{Z}_1+\sigma Z_2$，其中Z1、Z2是独立标准正态变量。因此 $$Y_{\mathrm{n}\mathrm{e}\mathrm{w}}\mid y\sim N(m_n,{\sigma }^2+v_n)=N(1,14/3).$$ 总预测方差等于噪声方差四，加上尚未消除的均值不确定性2/3。若只把μ估计为一后用 $N(1,4)$，会遗漏后一部分。给μ的后验区间与给新观测的预测区间具有不同宽度和对象，不应混称。

变式。 再获得一次测量三。把原后验 $N(1,2/3)$ 作为新先验，精度增加1/4，得到新精度7/4、方差4/7，均值 $(1/(2/3)+3/4)/(7/4)=9/7$。直接用三次测量总和九更新，也得到 $(9/4)/(7/4)=9/7$。这个一致性来自同一条件独立模型下似然相乘，不是每增加一条记录就任意给旧结论相同权重。

Python 3

import math

def normal_mean_update(observations, prior_mean, prior_variance, noise_variance):
    data = tuple(observations)
    if (not math.isfinite(prior_mean)
            or any(not math.isfinite(v) or v <= 0 for v in (prior_variance, noise_variance))
            or any(not math.isfinite(v) for v in data)):
        raise ValueError("finite observations and positive finite variances required")
    try:
        precision = 1 / prior_variance + len(data) / noise_variance
        if not math.isfinite(precision) or precision <= 0:
            raise ArithmeticError("posterior precision exceeds floating range")
        variance = 1 / precision
        weighted_sum = prior_mean / prior_variance + math.fsum(data) / noise_variance
        mean = variance * weighted_sum
        predictive_variance = variance + noise_variance
    except (OverflowError, ZeroDivisionError) as error:
        raise ArithmeticError("posterior calculation exceeds floating range") from error
    if (not math.isfinite(weighted_sum) or not math.isfinite(mean)
            or any(not math.isfinite(v) or v <= 0 for v in (variance, predictive_variance))):
        raise ArithmeticError("posterior result exceeds floating range")
    return mean, variance, predictive_variance

m, v, predictive_v = normal_mean_update([2, 4], 0, 1, 4)
assert abs(m - 1) < 1e-12 and abs(v - 2 / 3) < 1e-12
assert abs(predictive_v - 14 / 3) < 1e-12
sequential = normal_mean_update([3], m, v, 4)
joint = normal_mean_update([2, 4, 3], 0, 1, 4)
assert all(abs(a - b) < 1e-12 for a, b in zip(sequential, joint))
assert abs(joint[0] - 9 / 7) < 1e-12 and abs(joint[1] - 4 / 7) < 1e-12
assert normal_mean_update([], 0, 1, 4) == (0, 1, 5)
print(m, v, predictive_v, joint[:2])

代码接受有限测量值和有限正方差，返回后验均值、后验方差、下一次同噪声观测的预测方差；不修改输入。若派生精度、加权和或结果超出浮点可表示范围，或方差下溢成零，程序显式抛出 ArithmeticError。可表示的空数据更新返回原先验及相应预测方差。时间与数据条数成正比；为复制输入使用线性空间，可在数据很大时改为只累积条数和总和的常数额外空间实现。

损失、先验与存在条件

贝叶斯估计不是一律求后验众数。连续后验下，精确猜中一个实数的零一损失几乎处处等于一，不能直接用它证明连续MAP最优；若要用小区间命中或其他容差损失，须明确该损失与坐标。对后验均值使用平方损失时，还必须有相应矩存在。例如某些重尾后验没有有限均值或二阶矩，配方证明的风险分解就不能不加检查地使用。

本例高斯先验和已知方差保证所用后验、矩与风险均有限。未知噪声方差、相关测量、异常值机制和先验取值改变后，不能继续直接使用同一更新公式。选择先验、检验模型与优化给定后验风险，是需要分别完成的工作。