Scipy.optimize.differential_evolution能否计算估计参数的SD、SE及参数不确定度

scipy.optimize.differential_evolution属于启发式全局优化算法，本身没有内置协方差矩阵计算逻辑，你可以通过以下3种常用方法计算参数不确定度：

1. 拔靴法（Bootstrap，优先推荐） #

该方法属于非参数方法，不需要对残差分布、模型线性性做假设，结果可靠性最高：

• 准备原始观测数据集（自变量x、观测值y_obs）、差分进化拟合得到的最优参数θ_opt、自定义模型函数y_model(x, θ)
• 计算最优参数对应的残差序列：res = y_obs - y_model(x, θ_opt)
• 进行多次重采样：每次从残差序列中有放回抽取和原样本量一致的残差，构造伪观测值y_boot = y_model(x, θ_opt) + 随机抽取的残差
• 对每一组伪观测值重新用差分进化拟合，得到一组参数估计结果
• 重复上述过程500~2000次，得到所有参数的经验分布，分布的标准差即为对应参数的不确定度，也可以直接取分位数得到参数的置信区间

简化代码示例：

import numpy as np
from scipy.optimize import differential_evolution

# 自定义模型示例（此处为线性模型，可替换为你的实际模型）
def model(x, a, b):
    return a * x + b

# 你的原始观测数据
x = np.linspace(0, 10, 100)
y_obs = 2 * x + 1 + np.random.normal(0, 0.5, size=100)

# 差分进化拟合得到最优参数
def loss_func(params):
    a, b = params
    return np.sum((y_obs - model(x, a, b)) ** 2)
bounds = [(-10, 10), (-10, 10)] # 替换为你的参数实际取值范围
theta_opt = differential_evolution(loss_func, bounds).x
res_opt = y_obs - model(x, *theta_opt)

# Bootstrap计算不确定度
n_boot = 1000
params_boot = np.zeros((n_boot, len(theta_opt)))
for i in range(n_boot):
    res_sample = np.random.choice(res_opt, size=len(res_opt), replace=True)
    y_boot = model(x, *theta_opt) + res_sample
    def loss_boot(params):
        a, b = params
        return np.sum((y_boot - model(x, a, b)) ** 2)
    params_boot[i] = differential_evolution(loss_boot, bounds).x

# 输出每个参数的标准不确定度
param_uncertainty = np.std(params_boot, axis=0)
print(f"参数1不确定度：{param_uncertainty[0]:.3f}, 参数2不确定度：{param_uncertainty[1]:.3f}")

2. 雅可比/海森矩阵近似法 #

该方法计算速度更快，但仅适用于模型在最优参数附近接近线性、残差服从正态同分布的场景：

• 拿到最优参数θ_opt后，用数值差分方法计算模型在最优参数处对所有待估参数的雅可比矩阵J
• 计算残差的均方误差：sigma_sq = np.sum(res_opt**2) / (样本量 - 待估参数个数)
• 按非线性最小二乘协方差公式计算：cov = sigma_sq * np.linalg.inv(J.T @ J)，协方差矩阵对角元开根号即为各参数的标准不确定度

提示：可以用scipy.optimize.approx_fprime快速完成数值差分计算雅可比矩阵，无需手动推导偏导。

3. MCMC后验采样法 #

如果采用贝叶斯估计框架，可以用该方法得到参数的完整后验分布：

• 把差分进化得到的最优参数作为MCMC采样的初始值，为所有参数设定合理的先验分布（无先验信息可设为均匀分布）
• 运行MCMC采样（可借助emcee等第三方库）得到参数的后验分布，后验分布的标准差即为参数不确定度

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