线性最小二乘法中矩阵不确定度估计

1. 引言

在数据处理和参数估计中，线性最小二乘法是一种常用的方法。当我们使用线性最小二乘法估计参数矩阵时，了解估计结果的不确定度对于评估模型的可靠性和精度非常重要。本技术文档将详细介绍如何在 MATLAB 中计算参数矩阵的估计值及其不确定度，包括相关的理论推导和代码示例。

2. 问题描述

给定两个矩阵 $\mathbf{A}$ 和 $\mathbf{B}$ ，它们的关系可以表示为：

$\mathbf{B} = \mathbf{R}_a \mathbf{A} + \mathbf{E}$

其中：

$\mathbf{B} \in \mathbb{R}^{3 \times N}$ ：观测矩阵。
$\mathbf{A} \in \mathbb{R}^{3 \times N}$ ：已知矩阵。
$\mathbf{R}_a \in \mathbb{R}^{3 \times 3}$ ：待估计的参数矩阵。
$\mathbf{E} \in \mathbb{R}^{3 \times N}$ ：残差矩阵或误差项。
$N$ ：观测值的数量。

我们的目标是估计参数矩阵 $\mathbf{R}_a$ ，并计算其不确定度。

3. 线性最小二乘估计

3.1 参数矩阵的估计

使用线性最小二乘法，参数矩阵的估计值为：

$\hat{\mathbf{R}}_a = \mathbf{B} \mathbf{A}^\top (\mathbf{A} \mathbf{A}^\top)^{-1}$

3.2 MATLAB 实现

% 假设已知矩阵 A 和 B
% A 和 B 的维度为 3 x N% 计算参数矩阵的估计值
Ra_hat = B * A' * inv(A * A');

4. 不确定度估计

4.1 计算残差矩阵

首先，计算残差矩阵 $\mathbf{E}$ ：

$\mathbf{E} = \mathbf{B} - \hat{\mathbf{R}}_a \mathbf{A}$

MATLAB 实现

% 计算残差矩阵
E = B - Ra_hat * A;

4.2 估计残差方差

假设误差项是独立同分布的，具有相同的方差 $\sigma^2$ ，则可以估计残差的方差：

$\hat{\sigma}^2 = \frac{1}{N - p} \sum_{i=1}^N \|\mathbf{e}_i\|^2$

其中：

$\mathbf{e}_i$ 是残差矩阵 $\mathbf{E}$ 的第 $i$ 列。
$p$ 是估计参数的数量，对于 $\times 3$ 的矩阵， $p = 9$ 。

MATLAB 实现

% 计算残差的平方和
residuals_squared = sum(E.^2, 'all');% 计算参数数量 p 和观测数量 N
p = numel(Ra_hat);
N = size(A, 2);% 估计残差方差
sigma2 = residuals_squared / (N - p);

4.3 计算参数矩阵的协方差矩阵

参数矩阵的协方差矩阵为：

$\text{Cov}(\hat{\mathbf{R}}_a) = \hat{\sigma}^2 (\mathbf{A} \mathbf{A}^\top)^{-1}$

MATLAB 实现

% 计算协方差矩阵
Cov_Ra = sigma2 * inv(A * A');

4.4 提取标准不确定度（标准差）

参数矩阵中每个元素的标准不确定度为协方差矩阵对角线元素的平方根：

$\text{不确定度} = \sqrt{\text{diag}(\text{Cov}(\hat{\mathbf{R}}_a))}$

MATLAB 实现

% 提取标准不确定度
standard_uncertainty = sqrt(diag(Cov_Ra));

5. 理论解释

5.1 克罗内克积的省略

在一般情况下，当我们将矩阵参数展开为向量时，协方差矩阵的计算需要使用克罗内克积：

$\text{Cov}(\text{vec}(\hat{\mathbf{R}}_a)) = \hat{\sigma}^2 (\mathbf{A} \mathbf{A}^\top)^{-1} \otimes \mathbf{I}_3$

其中：

$\text{vec}(\hat{\mathbf{R}}_a)$ 是将矩阵 $\hat{\mathbf{R}}_a$ 展开为列向量。
$\otimes$ 表示克罗内克积。
$\mathbf{I}_3$ 是 $\times 3$ 的单位矩阵。

然而，在我们的计算中，由于矩阵维度匹配，且直接使用矩阵形式计算，克罗内克积的作用被隐含地包含在矩阵运算中，因此可以省略。

5.2 协方差矩阵的含义

计算得到的协方差矩阵 $\text{Cov}(\hat{\mathbf{R}}_a)$ 包含了参数矩阵中元素的不确定度信息。对角线元素表示各参数的方差，非对角线元素表示参数之间的协方差。

6. 完整的 MATLAB 代码示例

% 假设已知矩阵 A 和 B
% A 和 B 的维度为 3 x N% 1. 计算参数矩阵的估计值
Ra_hat = B * A' * inv(A * A');% 2. 计算残差矩阵
E = B - Ra_hat * A;% 3. 计算残差的平方和
residuals_squared = sum(E.^2, 'all');% 4. 计算参数数量 p 和观测数量 N
p = numel(Ra_hat);
N = size(A, 2);% 5. 估计残差方差
sigma2 = residuals_squared / (N - p);% 6. 计算协方差矩阵
Cov_Ra = sigma2 * inv(A * A');% 7. 提取标准不确定度
standard_uncertainty = sqrt(diag(Cov_Ra));% 8. 显示结果
disp('估计的参数矩阵 Ra_hat：');
disp(Ra_hat);disp('参数矩阵 Ra_hat 的标准不确定度：');
disp(standard_uncertainty);

7. 示例数据与结果

7.1 示例数据

为了演示上述方法，我们使用一些示例数据。

% 生成示例数据
N = 100; % 观测数量
A = randn(3, N); % 已知矩阵 A
true_Ra = [1.2, 0.5, -0.3; -0.4, 0.8, 0.6; 0.2, -0.7, 1.1]; % 真实的参数矩阵
noise = 0.05 * randn(3, N); % 噪声
B = true_Ra * A + noise; % 观测矩阵 B% 使用上述步骤计算估计值和不确定度

7.2 计算结果

运行完整的 MATLAB 代码后，我们可以得到估计的参数矩阵和标准不确定度。

% 输出示例结果（具体数值根据运行结果而定）
估计的参数矩阵 Ra_hat：1.1985    0.5012   -0.2987-0.4023    0.8021    0.59980.1989   -0.7005    1.1007参数矩阵 Ra_hat 的标准不确定度：0.00500.00500.0050

8. 结论

通过以上步骤，我们成功地使用线性最小二乘法估计了参数矩阵 $\mathbf{R}_a$ ，并计算了其不确定度。关键步骤包括计算残差、估计残差方差、计算协方差矩阵以及提取标准不确定度。

在实际应用中，了解参数估计的不确定度对于评估模型的可靠性和进一步的数据分析具有重要意义。

9. 注意事项

假设条件：以上方法假设误差项是独立同分布的，且具有相同的方差。
数据质量：观测数据的质量直接影响参数估计的准确性和不确定度。
矩阵求逆：在计算 $(\mathbf{A} \mathbf{A}^\top)^{-1}$ 时，需要确保矩阵 $\mathbf{A} \mathbf{A}^\top$ 是可逆的，即矩阵 $\mathbf{A}$ 的行满秩。

10. 参考文献

Björck, Å. (1996). Numerical Methods for Least Squares Problems. SIAM.
Menke, W. (2012). Geophysical Data Analysis: Discrete Inverse Theory (3rd ed.). Academic Press.
Strang, G. (2009). Introduction to Linear Algebra (4th ed.). Wellesley-Cambridge Press.

附录：符号说明