参数估计、似然与贝叶斯更新

拟合参数不只是找曲线上的最大值。必须说明曲线代表什么数据机制、哪些量未知，以及所报告结果采用什么损失或不确定性解释。本课推导 Bernoulli 与高斯最大似然，区分偏差和方差，完整归一化 beta–Bernoulli 更新，并从噪声和先验推导正则化系数。似然、抽样分布与后验是三个不同数学对象。

检索准备：模块 16的导数、模块 17的密度积分、模块 23的矩，以及模块 24的抽样极限。实验仅使用 Python 标准库。即使观测离散，连续参数与高斯例子仍明确需要微积分。

统计模型是一组概率规律 P_θ，以参数 θ∈Θ 索引。模型正确指定时，数据服从其中一个规律。Bernoulli 的 Θ=[0,1]，θ=p 为成功率；已知正方差的高斯位置模型中，Θ 为实数，θ=μ。如果方差也未知，参数是 (μ,σ²)，σ²>0。对 μ 为常数的因子，此时可能依赖未知参数。

频率学派抽样分析把 θ 视为固定未知量，样本及统计量随机。估计量 T(X₁,…,X_n) 是看数据之前指定的规则，其实现输出 t 是估计值。估计量有随 θ 改变的抽样分布，估计值则只是一次结果。E_θ[T]=θ 不能改写为每次 t 都等于真值；无偏性从不这样保证。

贝叶斯分析给参数先验，通过联合模型和数据更新。给定参数时各观测可能独立，积分掉共同未知参数后却可能依赖。因此独立声明必须说明条件层次。后验均值也是数据规则，可在每个固定 θ 下计算频率学派抽样风险，而不将该风险等同于后验概率。

可识别性（identifiability）要求不同参数对应不同可观测规律。若 θ≠θ′ 却 P_θ=P_θ′，再多数据也无法仅凭这个规律识别标签。X∼N(θ²,1)、θ∈R 中，θ 与 −θ 具有相同分布，符号不可识别。将 Θ 限制为非负可使参数到均值的映射单射，但零附近的不确定性仍需认真处理。

抽样方案属于模型。IID Bernoulli 似然需要共同率与给定率下的独立试验。复制行或变化率不因数据为零一就满足该乘积模型。依赖未观测结果的选择也会改变观测样本的似然。模型只是现实近似时，拟合参数可能代表模型内最佳近似，而非真实总体机制；后续回归通过残差与校准检查错设证据。

估计报告应写可观测变量、抽样单位、参数空间、依赖条件、规则和目标。这些内容决定端点是否允许、导数方程是否足够以及偏差由哪个期望定义。无约束公式给出空间外的值无效；合法边界最优也不能仅因没有内部导数零点就被丢弃。

数据 x 已实现时，L(θ;x)=p_θ(x)；离散情形用质量函数，连续情形用相对于共同参考测度的密度。表达式可与概率模型相同，但似然固定 x、改变 θ；概率规律固定 θ、遍历 x。似然在 Θ 上不会自动积分为 1，将其按某参数坐标归一化也是额外建模选择。

给定参数独立时 L=∏p_θ(x_i)，对数似然 ℓ=Σlog p_θ(x_i)。在正似然处，对数严格递增，所以最大点相同；参数下不可能的数据给 L=0、ℓ=−∞。对数将乘法变加法并避免很多下溢，却不修复错误独立条件。端点应使用明确扩展值，不能让零乘 log 零意外成为 NaN。

顺序 Bernoulli 样本有 k 成功、n−k 失败，L=p^k(1−p)^(n−k)。只记录 K=k 则乘上 binom(n,k)。固定 n,k 时该因子不依赖 p，所以两者的 p 最大点以及同先验下的后验形状相同。但事件不同，数据概率也不同。一个因子可在参数拟合中消去，却仍在计算概率或比较模型时重要。

省略常数以前必须说明优化或归一化哪个自变量。已知 σ 的高斯均值拟合中，前因子对 μ 恒定；拟合 σ 时，σ^(−n) 不可丢。后验中 θ 无关因子可从分子和证据消去，但仍需有限有效归一化。模型比较中的模型相关因子，也可能影响结果，即使它不改变每个模型内的参数最优点。

有限候选表可用 log-sum-exp：log Σexp a_j=m+log Σexp(a_j−m)，m=max a_j 有限。这稳定算术，却不自动创造候选后验；还须指定先验质量与候选模型。连续证据为积分，任意网格和还需求积权重及正确近似分析。

连续数据的似然值是密度，不是精确实数点的概率。密度可大于 1，数值单位随数据变换改变。可逆数据变换需要 Jacobian，若其不依赖 θ，则相对拟合不变。参数变换则不同：似然按坐标复合，但先验、后验密度还要 Jacobian。第六节据此说明连续 MAP 依坐标而变。

0<k<n 时，ℓ=k log p+(n−k)log(1−p)，0<p<1。导数 k/p−(n−k)/(1−p)=0，解得 p_hat=k/n。二阶导 −k/p²−(n−k)/(1−p)²<0，证明严格凹与唯一内部最大点。两种结果均出现时端点似然为零。完整推导要检查曲率和允许边界，不能仅因找到导数根就认定全局最优。

k=0、n>0 时 L=(1−p)^n，在 [0,1] 的最大点为 0；k=n 时最大点为 1。在开空间 (0,1) 上，只有趋近缺失端点的上确界，没有最大点。无数据 n=0 时各参数似然均为 1，不能由数据选唯一值。相同代数式在不同参数空间可有不同存在结论。

IID 高斯、已知 σ²>0 时，完整 logL=−n log(√(2π)σ)−Σ(x_i−μ)²/(2σ²)。拟合 μ 等价于最小化残差平方和。恒等式 Σ(x_i−μ)²=Σ(x_i−x_bar)²+n(μ−x_bar)² 由中心残差和为零得到。第二项非负，只在 μ=x_bar 时零，故 n>0 的全局均值 MLE 是样本平均；导数 n(x_bar−μ)/σ² 也给相同解。

方差 v 也未知且中心残差和 R>0 时，代入 μ=x_bar，目标为 −(n/2)log v−R/(2v) 加常数。导数 −n/(2v)+R/(2v²) 给 v_hat=R/n；其下递增、其上递减，两个方差极端都趋 −∞，所以是全局最大。MLE 分母 n，与下一节无偏方差的 n−1 目的不同。

如果高斯观测全相同，R=0。令 μ 等于该值、v→0⁺，似然无上界，因此正方差空间没有有限 MLE。代入 v_hat=0 越出模型。数值方差下限可以是实际限制，但必须说明它改变问题。简单模型中的唯一 MLE 也不能推出所有混合、非线性或高维模型都有唯一解。

这里平方损失来自高斯观测噪声。IID 已知尺度 b>0 的 Laplace 位置噪声，似然正比 exp(−Σ|x_i−μ|/b)，因此样本中位数最小化绝对偏差。偶数观测可以有整个中位区间。选择损失既是统计假设，也是计算选择。

固定标量 θ、估计量 T 有限二阶矩时，偏差为 E_θT−θ，方差为 E_θ[(T−E_θT)²]，MSE 为 E_θ[(T−θ)²]。将误差拆成中心项与偏差，平方取期望，交叉项因中心项均值零而消失，故 MSE=Var+Bias²。它们是在指定真参数下重复样本的性质，不由一次绝对误差决定，排序也可能随 θ 改变。

IID Bernoulli 的 K/n 无偏，方差和 MSE 都是 p(1−p)/n。固定 Beta(a,b) 后验均值规则 T_B=(K+a)/(n+a+b) 的频率期望为 (np+a)/(n+a+b)，偏差 [a−(a+b)p]/(n+a+b)，方差 np(1−p)/(n+a+b)²。收缩降低方差，却引入参数相关偏差。是否改善 MSE 要看总和，不能仅看方差或“贝叶斯”名称。

IID 有限方差 X_i 的 R=Σ(X_i−x_bar)²=Σ(X_i−μ)²−n(x_bar−μ)²。取期望得 E[R]=nσ²−n(σ²/n)=(n−1)σ²。因此 n>1 的 S²=R/(n−1) 无偏。共同均值、方差以及交叉协方差零是推导条件；相关行可使修正失效。该期望恒等式不需高斯，但特定精确抽样分布可能需要。

(1,2,4) 的均值 7/3、R=14/3，方差 MLE 14/9、无偏方差 7/3，后者大 n/(n−1) 倍。无偏是针对方差，不是其平方根；非线性开方一般改变偏差。n−1 也不保证比所有有偏替代有较小 MSE，应先声明评价准则。

一致性指 T_n 依概率收敛到 θ。Bernoulli 比例由弱法则一致；固定正 a,b 时，T_B−K/n=[a−(a+b)K/n]/(n+a+b)，因 K/n∈[0,1] 而一致趋零，所以该后验均值规则也一致。有限样本偏差与一致性可共存。若先验强度随 n 成比例增加，可能保留非零极限收缩，不能照搬此结论。

无偏也不推出一致：一直只返回第一个 Bernoulli 观测，虽无偏，但方差 p(1−p) 不随 n 缩小。消失偏差与消失方差则可通过平方误差概率界证明一致。不能未经条件和证明将样本均值的性质转给所有 MLE 算法；一般渐近理论还涉及可识别、正则性及边界。

风险比较应固定数据预算与估计规则。看完同一结果才选择规则，构成一个新数据依赖估计量，其期望不自动等于原来任一规则公式。实验 2 在生成前指定规则与合成率；下一课的模型选择方案使这种区别可操作。

Beta(a,b) 先验在 0<p<1 的密度为 p^(a−1)(1−p)^(b−1)/B(a,b)，a,b>0，B 为该分子在 [0,1] 的积分。正参数使可能的端点奇异仍可积。正整数参数时，重复分部积分给 B(a,b)=(a−1)!(b−1)!/(a+b−1)!。先验表达对率的模型选择，不是额外观察试验。

先验乘似然，幂指数相加，得到 p^(a+k−1)(1−p)^(b+n−k−1)。其积分 B(a+k,b+n−k) 归一化 Beta(a+k,b+n−k) 后验。顺序样本证据为更新 B 与原 B 的比，计数证据另含 binom(n,k)。该因子在后验中消去，但在观测计数概率中保留。条件化需要证据正且有限。

令 A=a+k、B=b+n−k。积分比 B(A+1,B)/B(A,B) 给均值 A/(A+B)；二阶矩 A(A+1)/[(A+B)(A+B+1)] 给方差 AB/[(A+B)²(A+B+1)]。这是此数据、此先验下 p 的条件不确定性，区别于固定 p 下对新数据平均的规则 (K+a)/(n+a+b) 抽样方差。

给定 p 后与旧试验独立的未来 Bernoulli 成功概率，后验预测为 ∫pπ(p|x)dp=E[p|x]=A/(A+B)，而非代入众数。多个未来试验给定 p 独立，积分后却有协方差 Var(p|x)>0，因为共同未知率。因此未来总计数一般不是后验均值率的独立二项，而是 beta-binomial 混合。

等尾 95% 可信区间取后验 CDF 的 .025 与 .975 分位，在所选模型先验下包含 .95 后验概率。实验 1 通过 CDF 质量核验数值反演。它不自动对每个固定 p 给 95% 频率覆盖；模块 26 将区分置信区间。均匀先验仅在一个特定坐标中均匀，小数据的结论可对先验与模型错设敏感。

r 个未来成功数 S 的条件均值为 rp、方差 rp(1−p)。全方差给 Var(S|x)=rE[p(1−p)|x]+r²Var(p|x)。写 m=A/(A+B)，化为 r m(1−m)+r(r−1)Var(p|x)。第一项是代入二项方差，第二项是共享未知率的额外波动；r=1 时第二项零，故单试验预测简单但联合律不同。Beta(10,4) 的十次未来成功均值 50/7，方差约 3.265306，而代入二项为 2.040816。它是条件预测，不是新训练数据间后验均值的波动。

探索器限定整数先验参数，以便不依赖外部统计库就透明计算 beta 密度与 CDF。两个面板分别绘制相对似然与归一化后验密度，纵轴独立且明确标记。相对似然高度 1 不是后验概率。MLE、后验均值、密度众数与区间分别报告，也处理全成功、全失败和均匀后验。

真正新的条件独立批次按顺序更新，与一次合并相同，因为计数相加。对同一批更新两次会平方似然、虚构证据：八/二数据从 Beta(10,4) 再用一次得到 Beta(18,6)，并非新数据的合理改进。若先验从同一数据估出，也需披露和相应分析，不能默称独立固定先验。

最大后验（MAP）选择后验质量众数，或指定连续坐标的密度众数。密度情形最大化 Lπ 等价于最小化负 logL 加负 log 先验，允许省参数无关常数。连续点概率为零，密度众数不是有正点概率的参数。平坦或多峰后验可有非唯一 MAP，也不必接近均值或适合预测。

Beta(A,B)、A,B>1 的 log 密度导数给众数 (A−1)/(A+B−2)。A=1、B>1 时，连续扩展到 [0,1] 的最大密度端点为零；B=1、A>1 时为一；两者皆一时均匀，全部点并列。参数小于一可有端点密度发散，应明确众数约定。探索器限整数形状至少一，避免隐藏奇异情况。

IID y_i∼N(μ,σ²)、已知 σ²>0，先验 μ∼N(0,τ²)、τ²>0，负 log 后验为 Σ(y_i−μ)²/(2σ²)+μ²/(2τ²) 加 μ 无关项。乘 2σ² 得 SSE+λμ²，λ=σ²/τ²。求导给 μ_MAP=n y_bar/(n+λ)。若使用平均平方误差，对应系数为 λ/n；平均数据项却不变系数，会改变先验等价目标。

独立零均值高斯系数先验产生平方欧氏范数惩罚，独立 Laplace 先验产生绝对值惩罚。这个解释同时依赖系数先验与观测噪声：高斯噪声给平方残差，Laplace 噪声给绝对残差。惩罚截距、改变特征单位或噪声尺度也改变概率解释。后续回归处理设计矩阵，此标量推导已固定系数约定。

连续 MAP 依参数坐标。若 p 后验为 Beta(A,B)，η=log[p/(1−p)]，η 密度含 dp/dη=p(1−p)，所以按 p 表达的形状正比 p^A(1−p)^B，η 众数对应 p=A/(A+B)。原 p 密度众数却是 (A−1)/(A+B−2)。正确变换的分布与事件概率一致，密度高度和众数依坐标。某坐标平坦不能直接称为中立先验。

MAP 正则化是建模解释，不保证不确定性校准或泛化。点优化器丢掉大量后验信息，验证集调出的系数需说明选择协议。σ 也拟合时要恢复归一化项；未知超参数改变先验时，其归一化也可能影响拟合。只保留熟悉损失而不检查变化的自变量，可悄悄改变模型。

实验 1 · 似然与归一化 beta 更新

下载 lab1_likelihood_and_beta_update.py

"""CPU-only: likelihood estimates and an integer-parameter beta posterior."""
import math


def beta_cdf(x, a, b):
    """For positive integer a,b, I_x(a,b) is a binomial upper tail."""
    if x <= 0:
        return 0.0
    if x >= 1:
        return 1.0
    n = a + b - 1
    return math.fsum(math.comb(n, j) * x**j * (1-x)**(n-j)
                     for j in range(a, n + 1))


def beta_quantile(u, a, b):
    lo, hi = 0.0, 1.0
    for _ in range(70):
        mid = (lo + hi) / 2
        if beta_cdf(mid, a, b) < u:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2


if __name__ == "__main__":
    n, k, a, b = 10, 8, 2, 2
    post_a, post_b = a+k, b+n-k
    mle, mean = k/n, post_a/(post_a+post_b)
    mode = (post_a-1)/(post_a+post_b-2)
    print("Ordered Bernoulli likelihood: p^8(1-p)^2; count factor cancels for p fitting")
    for p in (.2, .5, .7, .8, .9):
        print(f"p={p:.1f}: likelihood={p**k*(1-p)**(n-k):.9f}")
    low, high = (beta_quantile(u, post_a, post_b) for u in (.025, .975))
    assert abs(beta_cdf(high, post_a, post_b)-beta_cdf(low, post_a, post_b)-.95) < 1e-12
    print(f"MLE={mle:.6f}; posterior Beta({post_a},{post_b})")
    print(f"posterior mean / next-trial predictive={mean:.6f}; posterior mode={mode:.6f}")
    print(f"equal-tail 95% credible interval=[{low:.6f}, {high:.6f}] (numerical inversion)")
    xs = [1., 2., 4.]
    mu = math.fsum(xs)/len(xs)
    rss = math.fsum((x-mu)**2 for x in xs)
    print(f"Gaussian data {xs}: mean MLE={mu:.6f}; variance MLE={rss/3:.6f}; unbiased variance={rss/2:.6f}")
    print("All failures: Bernoulli MLE=0; Beta(2,12) mean=1/7, mode=1/12")

Ordered Bernoulli likelihood: p^8(1-p)^2; count factor cancels for p fitting
p=0.2: likelihood=0.000001638
p=0.5: likelihood=0.000976562
p=0.7: likelihood=0.005188321
p=0.8: likelihood=0.006710886
p=0.9: likelihood=0.004304672
MLE=0.800000; posterior Beta(10,4)
posterior mean / next-trial predictive=0.714286; posterior mode=0.750000
equal-tail 95% credible interval=[0.461868, 0.909080] (numerical inversion)
Gaussian data [1.0, 2.0, 4.0]: mean MLE=2.333333; variance MLE=1.555556; unbiased variance=2.333333
All failures: Bernoulli MLE=0; Beta(2,12) mean=1/7, mode=1/12

先推三个 Bernoulli 摘要，再核验区间后验质量，解释高斯方差分母与全失败端点。

实验 2 · 重复抽样偏差、方差和 MSE

下载 lab2_bias_variance_and_mse.py

"""Repeated independent datasets compare specified frequentist risks."""
import math
import random


def summarize(values, truth):
    avg = math.fsum(values)/len(values)
    var = math.fsum((v-avg)**2 for v in values)/len(values)
    mse = math.fsum((v-truth)**2 for v in values)/len(values)
    assert abs(mse-var-(avg-truth)**2) < 1e-12
    return avg-truth, var, mse


if __name__ == "__main__":
    rng = random.Random(25025)
    n, replications = 10, 6000
    print(f"seed=25025; independent synthetic datasets={replications}; n={n}; Beta(2,2) fixed before draws")
    for p in (.3, .95):
        counts = [sum(rng.random() < p for _ in range(n)) for _ in range(replications)]
        estimates = {"MLE": [k/n for k in counts], "posterior mean": [(k+2)/(n+4) for k in counts]}
        for name, values in estimates.items():
            bias, var, mse = summarize(values, p)
            model_bias = 0. if name == "MLE" else (2-4*p)/(n+4)
            model_var = p*(1-p)/n if name == "MLE" else n*p*(1-p)/(n+4)**2
            print(f"p={p:.2f}, {name}: empirical bias={bias:.6f}, variance={var:.6f}, MSE={mse:.6f}")
            print(f"  exact sampling bias={model_bias:.6f}, variance={model_var:.6f}, MSE={model_var+model_bias**2:.6f}")
    print("A posterior mean can be assessed by repeated-sample risk; its Bayesian interpretation does not guarantee lower risk at every fixed p.")

seed=25025; independent synthetic datasets=6000; n=10; Beta(2,2) fixed before draws
p=0.30, MLE: empirical bias=0.000417, variance=0.021185, MSE=0.021185
  exact sampling bias=0.000000, variance=0.021000, MSE=0.021000
p=0.30, posterior mean: empirical bias=0.057440, variance=0.010809, MSE=0.014108
  exact sampling bias=0.057143, variance=0.010714, MSE=0.013980
p=0.95, MLE: empirical bias=-0.000250, variance=0.004723, MSE=0.004723
  exact sampling bias=0.000000, variance=0.004750, MSE=0.004750
p=0.95, posterior mean: empirical bias=-0.128750, variance=0.002410, MSE=0.018986
  exact sampling bias=-0.128571, variance=0.002423, MSE=0.018954
A posterior mean can be assessed by repeated-sample risk; its Bayesian interpretation does not guarantee lower risk at every fixed p.

先算两个真率下的精确风险。有限模拟估计这些风险，不能证明统一排序；解释后验均值为何同时也是频率估计规则。

实验 3 · 数值与模型修复

下载 lab3_underflow_normalisation_and_prior.py

"""Diagnose probability underflow, parameter-dependent factors, and reused data."""
import math


def bernoulli_log_likelihood(p, k, n):
    if not 0 <= p <= 1 or not 0 <= k <= n:
        raise ValueError("invalid model or count")
    if (p == 0 and k > 0) or (p == 1 and n-k > 0):
        return -math.inf
    return (k*math.log(p) if k else 0.) + ((n-k)*math.log1p(-p) if n-k else 0.)


if __name__ == "__main__":
    n, k = 10000, 8000
    ll = [bernoulli_log_likelihood(p, k, n) for p in (.79, .8)]
    print(f"n={n}, k={k}: raw likelihoods={[math.exp(x) for x in ll]}")
    print(f"log likelihoods={ll[0]:.6f}, {ll[1]:.6f}; relative L(.79)/L(.8)={math.exp(ll[0]-ll[1]):.6f}")
    assert ll[1] > ll[0] and all(math.exp(x) == 0 for x in ll)
    for sigma in (1., 2.):
        full = -2*math.log(sigma)-math.log(2*math.pi)
        print(f"two zero residuals, sigma={sigma:.0f}: Gaussian full logL={full:.6f}; omitted-scale expression=0")
    a, b, k, n = 2, 2, 8, 10
    correct = (a+k)/(a+b+n)
    doubled = (a+2*k)/(a+b+2*n)
    print(f"One update: Beta(10,4), mean={correct:.6f}")
    print(f"Same data reused as new evidence: Beta(18,6), mean={doubled:.6f} (invalid independent-evidence claim)")
    sigma2, tau2, n, observed_mean = 1., 4., 4, 3.
    lam = sigma2/tau2
    map_mean = n*observed_mean/(n+lam)
    print(f"Gaussian mean, prior N(0,4), known noise variance1: MAP={map_mean:.6f}; sum-SSE penalty coefficient={lam:.6f}")
    print(f"For average-SSE convention the matching coefficient is lambda/n={lam/n:.6f}.")

n=10000, k=8000: raw likelihoods=[0.0, 0.0]
log likelihoods=-5007.074165, -5004.024235; relative L(.79)/L(.8)=0.047362
two zero residuals, sigma=1: Gaussian full logL=-1.837877; omitted-scale expression=0
two zero residuals, sigma=2: Gaussian full logL=-3.224171; omitted-scale expression=0
One update: Beta(10,4), mean=0.714286
Same data reused as new evidence: Beta(18,6), mean=0.750000 (invalid independent-evidence claim)
Gaussian mean, prior N(0,4), known noise variance1: MAP=2.823529; sum-SSE penalty coefficient=0.250000
For average-SSE convention the matching coefficient is lambda/n=0.062500.

修复下溢、遗漏尺度归一化、复用证据及总和/平均系数。指出每个优化需要保留哪些参数相关因子。

练习 1–12 必做，选学 13–14 在十二小时以外增加 35 分钟。

读 MIT 经典估计、贝叶斯推断及 Mathematics for Machine Learning 第八章。初等估计与 beta 更新均在本课推导。

推两个 MLE、方差分母、MSE 分解、beta 归一化与预测，以及高斯先验惩罚尺度。不依赖图解释可识别、端点存在与三种不确定性对象。

结业任务：比较 .8、5/7、.75 的不同目标，再修复重复更新与平均损失后不变的惩罚系数。

进入下一课的标准：能够在明确模型下推导并解释点估计。下一课添加置信覆盖、检验和实验协议。返回课程总览。