Python中拟合多输入多输出带协方差矩阵的模型实现

核心思路 #

普通的scipy.optimize.curve_fit无法直接处理输入输出的相关性与协方差矩阵，必须自定义联合负对数似然函数，把输入误差的传播、输出观测协方差全部纳入拟合逻辑。lmfit的Minimizer类对这类复杂拟合的支持更友好，也能直接输出参数协方差；Scipy则需手动处理优化与协方差估计。

步骤1：模拟带协方差的数据集 #

先生成符合需求的测试数据（与你的场景匹配）：

import numpy as np
import lmfit
from scipy.optimize import minimize

# 真实模型参数
true_a, true_b, true_c = 2.0, 1.5, 3.0
n_points = 50

# 生成无噪声的真实输入
x_true = np.random.uniform(-5, 5, n_points)
y_true = np.random.uniform(-5, 5, n_points)

# 输入协方差（含相关性）
x_err, y_err, x_y_corr = 0.2, 0.3, 0.6
input_cov = np.array([[x_err**2, x_y_corr*x_err*y_err],
                      [x_y_corr*x_err*y_err, y_err**2]])
# 带相关噪声的观测输入
x_obs, y_obs = np.random.multivariate_normal([x_true, y_true], input_cov, n_points).T

# 生成真实输出
u_true = true_a * x_true + true_b * y_true
v_true = true_b * x_true + true_c * y_true

# 输出协方差（含相关性）
u_err, v_err, u_v_corr = 0.4, 0.5, 0.7
output_cov = np.array([[u_err**2, u_v_corr*u_err*v_err],
                       [u_v_corr*u_err*v_err, v_err**2]])
# 带相关噪声的观测输出
u_obs, v_obs = np.random.multivariate_normal([u_true, v_true], output_cov, n_points).T

步骤2：基于lmfit的拟合实现 #

lmfit支持参数约束、自动协方差计算，是首选方案：

def vector_field_model(params, x, y):
    """定义2D向量场模型"""
    a = params['a']
    b = params['b']
    c = params['c']
    u = a * x + b * y
    v = b * x + c * y
    return u, v

def neg_log_likelihood(params, x_obs, y_obs, u_obs, v_obs, input_cov, output_cov):
    """联合负对数似然函数：纳入输入输出协方差与相关性"""
    total_nll = 0.0
    a, b, c = params['a'].value, params['b'].value, params['c'].value
    
    for idx in range(len(x_obs)):
        # 模型预测值
        u_pred, v_pred = vector_field_model(params, x_obs[idx], y_obs[idx])
        
        # 输入误差传播到输出的协方差（雅克比矩阵变换）
        jacobian = np.array([[a, b], [b, c]])
        model_cov_from_input = jacobian @ input_cov @ jacobian.T
        
        # 总输出协方差：输入误差传播项 + 观测输出协方差
        total_output_cov = model_cov_from_input + output_cov
        inv_cov = np.linalg.inv(total_output_cov)
        det_cov = np.linalg.det(total_output_cov)
        
        # 残差向量
        res = np.array([u_obs[idx] - u_pred, v_obs[idx] - v_pred])
        
        # 单个数据点的负对数似然
        nll = 0.5 * (np.log(det_cov) + res.T @ inv_cov @ res + 2*np.log(2*np.pi))
        total_nll += nll
    
    return total_nll

# 初始化参数
params = lmfit.Parameters()
params.add('a', value=1.0, min=0, max=5)
params.add('b', value=1.0, min=0, max=5)
params.add('c', value=2.0, min=0, max=5)

# 执行拟合
minimizer = lmfit.Minimizer(neg_log_likelihood, params, 
                            args=(x_obs, y_obs, u_obs, v_obs, input_cov, output_cov))
result = minimizer.minimize()

# 打印结果（含参数协方差）
lmfit.report_fit(result)

# 提取最优参数与协方差矩阵
best_params = {k: v.value for k, v in result.params.items()}
param_cov_matrix = result.covar

步骤3：基于Scipy的替代实现 #

如果偏好Scipy，可手动定义目标函数并优化：

def scipy_nll(params, x_obs, y_obs, u_obs, v_obs, input_cov, output_cov):
    a, b, c = params
    total_nll = 0.0
    
    for idx in range(len(x_obs)):
        u_pred = a * x_obs[idx] + b * y_obs[idx]
        v_pred = b * x_obs[idx] + c * y_obs[idx]
        
        jacobian = np.array([[a, b], [b, c]])
        model_cov_from_input = jacobian @ input_cov @ jacobian.T
        total_output_cov = model_cov_from_input + output_cov
        
        inv_cov = np.linalg.inv(total_output_cov)
        det_cov = np.linalg.det(total_output_cov)
        res = np.array([u_obs[idx] - u_pred, v_obs[idx] - v_pred])
        
        nll = 0.5 * (np.log(det_cov) + res.T @ inv_cov @ res + 2*np.log(2*np.pi))
        total_nll += nll
    
    return total_nll

# 初始参数猜测
initial_guess = [1.0, 1.0, 2.0]

# 执行优化（选择支持Hessian估计的方法）
scipy_result = minimize(scipy_nll, initial_guess, 
                        args=(x_obs, y_obs, u_obs, v_obs, input_cov, output_cov),
                        method='L-BFGS-B')

# 打印最优参数与协方差矩阵
print("最优参数:", scipy_result.x)
print("参数协方差矩阵:")
print(np.linalg.inv(scipy_result.hess_inv.todense()))

关键说明 #

• 输入误差通过雅克比矩阵传播为模型输出的不确定性，再与输出观测协方差合并，得到总误差协方差，这是准确拟合的核心。
• lmfit自动处理参数边界、协方差计算，比Scipy更省心；Scipy需手动指定优化方法并计算Hessian逆来获取参数协方差。
• 该方案能同时拟合影响多输出的参数（如示例中的b同时影响u和v，c影响v），解决单输出拟合无法获取c的问题。

内容的提问来源于stack exchange，提问作者Swike