基于线性回归的CRLB解释

4 minute read

Published: May 25, 2026

克拉美罗下界（Cramér-Rao Lower Bound, CRLB）是系统辨识里的重要概念，用于衡量参数估计的下界精度。

一、最大似然估计

Fisher信息由频率学派最大似然估计（Maximum Likelihood Estimation, MLE）而来。基于MLE进行参数估计：

\[L(\theta|x_1,x_2,\cdots,x_n) = \prod_{i=1}^{n}f(x_{i}\mid\theta)\] \[\hat \theta_{\text{MLE}} = \arg\max_{\theta} L(\theta|x_1,x_2,\cdots,x_n)\]

为方便计算，取对数形式，可得：

\[\log L(\theta|x_1,x_2,\cdots,x_n) = \log L(\theta|X) = \sum_{i=1}^{n} \log f(x_{i}\mid\theta) \\ \hat \theta_{\text{MLE}} = \arg\max_{\theta} \log L(\theta|x_1,x_2,\cdots,x_n) = \arg\max_{\theta} \sum_{i=1}^{n} \log f(x_{i}\mid\theta)\]

如果概率密度/质量函数关于$\theta$可导（该前提一般都满足），则MLE估计量满足以下条件：

\[\frac{\partial}{\partial \theta} \log L(\theta|x_1,x_2,\cdots,x_n) = \frac{\partial}{\partial \theta} \log L(\theta|X) = \frac{\partial}{\partial \theta} \ell(\theta|X) = 0\]

对于一般的线性回归模型：

\[Y = X^T \theta + \epsilon\]

如果误差项 $\epsilon$ 服从正态分布，则最小二乘（Ordinary Least Squares, OLS）所得参数估计结果与MLE相同（见最后附录），即：

\[\begin{aligned} \hat{\theta}_\text{OLS} &= \arg\min_{\theta} \sum_{i=1}^{n} (y_i - f(x_i, \theta))^2 \\ &= \arg\min_{\theta} \|Y - X^T\theta\|^2 \end{aligned}\] \[\frac{\partial}{\partial \theta} \|Y - X^T\theta\|^2 = -2X^T(Y - X^T\theta) = 0\]

则有：

\[\hat{\theta}_\text{OLS} = (X^TX)^{-1}X^TY\]

二、Fisher信息

令 $S(\theta

X)$ 为得分函数：

\[S(\theta|X) = \frac{\partial}{\partial \theta} \ell(\theta|X) = \frac{\partial}{\partial \theta} \log L(\theta|X)\]

对于一组给定的观测样本 $X$，其MLE估计量 $\hat{\theta}_X$ 满足：

\[S(\hat{\theta}_X|X) = 0\]

如果 $\theta^*$ 为真实参数，则其得分函数的期望为零：

\[\mathbb{E}_X[S(\theta^*|X)] = 0\]

而Fisher信息定义为得分函数的方差：

\[I(\theta) = \mathbb{E}_X \left[ \left( \frac{\partial}{\partial \theta} \log L(\theta|X) \right)^2 \right] = \mathbb{E}_X \left[ S(\theta|X)^2 \right]\]

在正则条件下

\[I(\theta) = -\mathbb{E}_X \left[ \frac{\partial^2}{\partial \theta^2} \log L(\theta|X) \right]\]

$I(\theta)$ 越大，表明样本对参数 $\theta$ 估计提供的信息越多，估计方差越小。由此引出Cramér-Rao下界（CRLB）：

任何正则无偏估计量的方差都受到Cramér-Rao下界的限制：

\[\mathrm{Var}(\hat{\theta}) \geq \frac{1}{I(\theta)}\]

而在大样本极限下，MLE的方差趋近于Cramér-Rao下界：

\[\sqrt{n}(\hat{\theta}_\text{MLE} - \theta) \xrightarrow{d} \mathcal{N}\left(0, \frac{1}{I(\theta)}\right)\]

三、Matlab中对CRLB的计算

Matlab系统辨识模型一般形式为：

\[y(t) = G(q, \theta)u(t) + H(q, \theta)e(t)\]

其中，$e(t)$ 为方差 $\lambda = \sigma^2$ 的白噪声。

接下来，Matlab采用观测Fisher信息矩阵来计算CRLB：

\[\hat{I}(\hat{\theta}) = - \frac{\partial^2 \ell(\hat{\theta})}{\partial \theta^2}|_{\theta=\hat{\theta}}\]

在Gauss-Newton或Levenberg-Marquardt优化算法中，该二阶导数Hessian矩阵被近似为雅可比矩阵的乘积，从而简化计算。

最后，计算持续激励大样本量条件下的CRLB和参数估计的协方差矩阵：

\[\text{Cov}(\hat{\theta}) \sim \hat{I}(\hat{\theta})^{-1}\]

四、基于线性回归的算例

这里使用一个简单的2x2线性回归模型作为算例，以阐述Fisher信息矩阵和CRLB的计算过程。该线性回归的方程为：

\[y = \theta_0 + \theta_1 x + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \Sigma)\]

其中，

\[\begin{aligned} \theta_0 &= [1, 2]^T \\ \theta_1 &= \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix} \\ \Sigma &= \begin{bmatrix} 0.1 & 0 \\ 0 & 0.1 \end{bmatrix} \end{aligned}\]

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

def gen_samples(N, sigma_1, sigma_2, theta_0, theta_1):
    Sigma = np.array([[sigma_1, 0], [0, sigma_2]])
    epsilon = np.random.multivariate_normal([0, 0], Sigma, N)

    x = np.random.rand(N, 2)
    y = theta_0 + x @ theta_1.T + epsilon
    return x, y


def plot_figs(x, y):
    N = x.shape[0]
    fig = plt.figure(figsize=(4, 4))
    subfigs = fig.subfigures(1, 2, wspace=0.1)

    axs_left = subfigs[0].subplots(2, 1)
    subfigs[0].suptitle("input x")
    axs_left[0].plot(range(N), x[:, 0], alpha=1, linewidth=0.6)
    axs_left[0].set_ylabel("x1")
    axs_left[0].grid(True, alpha=0.3)
    axs_left[1].plot(range(N), x[:, 1], alpha=1, linewidth=0.6)
    axs_left[1].set_ylabel("x2")
    axs_left[1].set_xlabel("Sample index")
    axs_left[1].grid(True, alpha=0.3)
    
    axs_right = subfigs[1].subplots(2, 1)
    subfigs[1].suptitle("output y")
    axs_right[0].plot(range(N), y[:, 0], alpha=1, linewidth=0.6, color="C1")
    axs_right[0].set_ylabel("y1")
    axs_right[0].grid(True, alpha=0.3)
    axs_right[1].plot(range(N), y[:, 1], alpha=1, linewidth=0.6, color="C1")
    axs_right[1].set_ylabel("y2")
    axs_right[1].set_xlabel("Sample index")
    axs_right[1].grid(True, alpha=0.3)

    plt.show()


def est_ols_param(x, y):
    N = x.shape[0]

    # 构造设计矩阵 X_design: [1, x1, x2] for each sample
    X_design = np.hstack([np.ones((N, 1)), x])

    # OLS估计: theta_hat = (X^T X)^{-1} X^T y
    theta_hat = np.linalg.inv(X_design.T @ X_design) @ X_design.T @ y

    # 提取估计的参数
    theta_0_hat = theta_hat[0, :]
    theta_1_hat = theta_hat[1:, :].T

    # 计算Fisher信息矩阵和CRLB
    # 计算残差
    y_pred = X_design @ theta_hat
    residuals = y - y_pred
    
    # 估计噪声协方差矩阵 Sigma_hat
    Sigma_hat = (residuals.T @ residuals) / N
    
    # 计算 Sigma 的逆矩阵
    Sigma_inv = np.linalg.inv(Sigma_hat)
    
    # 计算 X^T X
    XTX = X_design.T @ X_design
    
    # 计算Fisher信息矩阵
    # Fisher_info = (X^T X) ⊗ Sigma^{-1}
    Fisher_info = np.kron(XTX, Sigma_inv)
    
    # 计算CRLB (Cramér-Rao Lower Bound)
    CRLB = np.linalg.inv(Fisher_info)
    
    # 提取参数标准差（对角线元素的平方根）
    param_std = np.sqrt(np.diag(CRLB))

    return theta_0_hat, theta_1_hat, param_std

生成样本

theta_0 = np.array([1, 2])
theta_1 = np.array([[1, 2], [3, 4]])

N = 1000
sigma_1 = 0.1
sigma_2 = 0.1

# 生成样本
x, y = gen_samples(N, sigma_1, sigma_2, theta_0, theta_1)

画图：

# 画图
plot_figs(x, y)

基于OLS进行参数估计：

theta_0_hat, theta_1_hat, param_std = est_ols_param(x, y)

print(f"theta_0_hat = {theta_0_hat}")
print(f"theta_1_hat =\n{theta_1_hat}")
print(f"param_std = {param_std}")

对比试验：样本量对估计标准差的影响

接下来，分别对比不同样本量 $N = 100, 1000, 10000$ 时所得参数估计结果和标准差：

# 完成代码,将结果整理为表格
results = []

for N in [100, 1000, 10000]:
    x, y = gen_samples(N, sigma_1=1, sigma_2=1, theta_0=np.array([1, 2]), theta_1=np.array([[1, 2], [3, 4]]))
    theta_0_hat, theta_1_hat, param_std = est_ols_param(x, y)
    results.append({
        "N": N,
        "θ₀₁": f"{theta_0_hat[0]:.4g}",
        "θ₀₂": f"{theta_0_hat[1]:.4g}",
        "θ₁₁": f"{theta_1_hat[0, 0]:.4g}",
        "θ₁₂": f"{theta_1_hat[0, 1]:.4g}",
        "θ₂₁": f"{theta_1_hat[1, 0]:.4g}",
        "θ₂₂": f"{theta_1_hat[1, 1]:.4g}",
        "std(θ₀₁)": f"{param_std[0]:.4g}",
        "std(θ₀₂)": f"{param_std[1]:.4g}",
        "std(θ₁₁)": f"{param_std[2]:.4g}",
        "std(θ₁₂)": f"{param_std[3]:.4g}",
        "std(θ₂₁)": f"{param_std[4]:.4g}",
        "std(θ₂₂)": f"{param_std[5]:.4g}"
    })

df_results = pd.DataFrame(results)
df_results

可见，增加样本量能够显著降低参数估计的标准差，提高估计的精度。

对比试验：噪声对估计标准差的影响

对比 $y_1$ 通道上加入不同强度噪声 $\sigma_1 = 0.001, 0.1, 1$ 对参数估计标准差的影响：

# 完成代码,将结果整理为表格
results_noise = []

for sigma in [0.001, 0.1, 1]:
    x, y = gen_samples(N=10000, sigma_1=sigma, sigma_2=sigma_2, theta_0=np.array([1, 2]), theta_1=np.array([[1, 2], [3, 4]]))
    theta_0_hat, theta_1_hat, param_std = est_ols_param(x, y)
    results_noise.append({
        "σ": sigma,
        "θ₀₁": f"{theta_0_hat[0]:.4g}",
        "θ₀₂": f"{theta_0_hat[1]:.4g}",
        "θ₁₁": f"{theta_1_hat[0, 0]:.4g}",
        "θ₁₂": f"{theta_1_hat[0, 1]:.4g}",
        "θ₂₁": f"{theta_1_hat[1, 0]:.4g}",
        "θ₂₂": f"{theta_1_hat[1, 1]:.4g}",
        "std(θ₀₁)": f"{param_std[0]:.4g}",
        "std(θ₀₂)": f"{param_std[1]:.4g}",
        "std(θ₁₁)": f"{param_std[2]:.4g}",
        "std(θ₁₂)": f"{param_std[3]:.4g}",
        "std(θ₂₁)": f"{param_std[4]:.4g}",
        "std(θ₂₂)": f"{param_std[5]:.4g}"
    })

df_results_noise = pd.DataFrame(results_noise)
df_results_noise

从上述结果可见，当样本量充足时，输出 $y_1$ 通道上的噪声标准差越大，与该通道有关的参数估计结果的标准差也越大，而其他通道估计结果几乎不受影响。

附录

假设 $\epsilon \sim N(0, \sigma^2)$，则 $Y_i

X_i \sim N(X_i^T \theta, \sigma^2)$，对应的对数似然函数为：

\[\ell(\theta) = -\frac{n}{2} \log(2\pi\sigma^2) - \frac{1}{2\sigma^2} \sum_{i=1}^n (Y_i - X_i^T \theta)^2\]

在固定的 $\sigma^2$ 下，最大化对数似然函数等价于最小化残差平方和，即：

\[\arg\min_{\theta} \sum_{i=1}^n (Y_i - X_i^T \theta)^2\]

因此，对于线性模型，如果误差项服从正态分布，则OLS估计量等价于MLE估计量。

Share on

WeChat

罗磊

基于线性回归的CRLB解释

一、最大似然估计

二、Fisher信息

三、Matlab中对CRLB的计算

四、基于线性回归的算例

对比试验：样本量对估计标准差的影响

对比试验：噪声对估计标准差的影响

附录

Share on

You May Also Enjoy

马尔可夫链定阶

工业过程数据稳态判断

小样本马尔可夫链独立性检验

多尺度时序因果分析