基于Scipy最小二乘法拟合耦合方程参数及解决报错问题

1. 修正残差函数（解决TypeError） #

scipy.optimize.leastsq要求输入的残差函数必须接收参数数组并返回残差向量，而非直接计算平方和的数值。你的错误在于将foo定义为一个浮点数（平方和结果），而不是可调用的函数。

正确的残差函数应该：

• 接收参数数组K（包含K11、K12、K21、K22）
• 计算模型f和g的预测值
• 返回实验值与预测值的残差（将两组残差合并为一维数组）

代码实现：

import numpy as np
from scipy.optimize import leastsq

# 定义原模型函数
def f(L, K11, K12):
    return (K11*L + 2*K11*K12*L**2) / (2*(1 + K11*L + K11*K12*L**2))

def g(L, K21, K22):
    return (K21*L + 2*K21*K22*L**2) / (2*(1 + K21*L + K21*K22*L**2))

# 正确的残差函数
def residuals(K, L, Y1_exp, Y2_exp):
    K11, K12, K21, K22 = K
    # 计算模型预测值
    f_pred = f(L, K11, K12)
    g_pred = g(L, K21, K22)
    # 返回合并后的残差向量
    return np.concatenate([Y1_exp - f_pred, Y2_exp - g_pred])

2. 执行最小二乘法拟合 #

假设你已经有实验数据L（自变量）、Y1_exp（Y1实验值）、Y2_exp（Y2实验值），以及初始参数K0：

# 示例初始参数（需根据你的数据调整）
K0 = np.array([0.1, 0.05, 0.12, 0.06])

# 执行拟合
fit_result = leastsq(residuals, K0, args=(L, Y1_exp, Y2_exp))
K_fit = fit_result[0]  # 拟合得到的参数数组 [K11_fit, K12_fit, K21_fit, K22_fit]

# 计算参数的标准误差（基于协方差矩阵）
cov_x = fit_result[1]
if cov_x is not None:
    # 残差方差估计
    residual_variance = (residuals(K_fit, L, Y1_exp, Y2_exp)**2).sum() / (len(Y1_exp) + len(Y2_exp) - len(K_fit))
    # 缩放协方差矩阵
    cov_scaled = cov_x * residual_variance
    # 提取参数标准误差
    std_errors_leastsq = np.sqrt(np.diag(cov_scaled))
else:
    std_errors_leastsq = None
    print("无法计算协方差矩阵，建议使用Bootstrap方法")

3. Bootstrap方法计算参数误差 #

Bootstrap通过重复采样实验数据并重新拟合，得到参数的分布，进而估计误差：

def bootstrap_param_errors(L, Y1_exp, Y2_exp, K0, n_iter=1000):
    n_samples = len(L)
    param_samples = []
    
    for _ in range(n_iter):
        # 随机采样数据索引（带替换）
        indices = np.random.choice(n_samples, size=n_samples, replace=True)
        L_boot = L[indices]
        Y1_boot = Y1_exp[indices]
        Y2_boot = Y2_exp[indices]
        
        # 对Bootstrap样本拟合
        try:
            boot_result = leastsq(residuals, K0, args=(L_boot, Y1_boot, Y2_boot))
            param_samples.append(boot_result[0])
        except:
            # 跳过拟合失败的样本
            continue
    
    # 转换为数组
    param_samples = np.array(param_samples)
    # 计算每个参数的标准差（作为误差估计）
    boot_std_errors = np.std(param_samples, axis=0)
    
    return boot_std_errors, param_samples

# 执行Bootstrap
boot_std_errors, param_samples = bootstrap_param_errors(L, Y1_exp, Y2_exp, K0, n_iter=1000)

结果输出 #

print("拟合参数：")
print(f"K11: {K_fit[0]:.4f} (Leastsq误差: {std_errors_leastsq[0]:.4f} | Bootstrap误差: {boot_std_errors[0]:.4f})")
print(f"K12: {K_fit[1]:.4f} (Leastsq误差: {std_errors_leastsq[1]:.4f} | Bootstrap误差: {boot_std_errors[1]:.4f})")
print(f"K21: {K_fit[2]:.4f} (Leastsq误差: {std_errors_leastsq[2]:.4f} | Bootstrap误差: {boot_std_errors[2]:.4f})")
print(f"K22: {K_fit[3]:.4f} (Leastsq误差: {std_errors_leastsq[3]:.4f} | Bootstrap误差: {boot_std_errors[3]:.4f})")

关键说明 #

• 残差函数的核心要求：必须返回一维残差数组，而非平方和数值，这是解决TypeError的关键。
• Bootstrap优势：无需假设残差服从正态分布，对非线性模型的误差估计更可靠，尤其是当数据量较小时。
• 初始参数选择：初始参数K0需尽量接近真实值，否则leastsq可能收敛到局部最优解，可通过绘制实验数据与模型曲线的初始拟合来调整。

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