如何在双对数坐标轴下实现指数函数的稳健自动化拟合

核心原理 #

当前使用的指数模型可通过对数变换转换为线性模型，无需手动设置初始参数：
对于目标拟合函数 y = exponent1 * a * exp(- exponent2 * b * x)，两边取自然对数可得：
ln(y) = ln(exponent1 * a) - b * exponent2 * x
对拟合区间内的x和y做线性回归即可直接得到a、b的最优初始值，传入curve_fit即可实现全自动拟合。

优化后可批量运行的代码 #

import numpy as np
import scipy.optimize as sc
import scipy.stats as stats
import matplotlib.pyplot as plt

# 固定参数配置
start_exp = 21  # 拟合区间起点
end_exp   = 40  # 拟合区间终点
exponent1 = 3.      
exponent2 = 2.    

# 原始数据
bin_centers = np.array([9.31939514e-01, 1.71773809e+00, 2.58227857e+00, 3.16611120e+00,
       3.88194374e+00, 4.75962032e+00, 5.83573259e+00, 7.15514527e+00,
       8.77286667e+00, 1.07563420e+01, 1.31882653e+01, 1.61700271e+01,
       1.98259415e+01, 2.43084291e+01, 2.98043714e+01, 3.65429025e+01,
       4.48049619e+01, 5.49350072e+01, 6.73553750e+01, 8.25838892e+01,
       1.01255449e+02, 1.24148500e+02, 1.52217488e+02, 1.86632650e+02,
       2.28828807e+02, 2.80565179e+02, 3.43998732e+02, 4.21774108e+02,
       5.17133876e+02, 6.34053729e+02, 7.77408231e+02, 9.53174043e+02,
       1.16867911e+03])

norm_counts = np.array([7.47683145e-02, 2.44802278e-02, 6.40790461e-02, 4.24335863e-02,
       2.91339699e-02, 2.34427016e-02, 1.63884469e-02, 2.19590994e-02,
       1.42326906e-02, 1.64066907e-02, 1.03597068e-02, 1.03504798e-02,
       9.22668847e-03, 6.22943476e-03, 6.37955519e-03, 5.30701726e-03,
       3.61688252e-03, 3.48879003e-03, 2.46230484e-03, 1.91634352e-03,
       1.45802119e-03, 9.72116428e-04, 8.17790259e-04, 4.57538385e-04,
       3.66533958e-04, 2.17783386e-04, 1.36803586e-04, 7.28849304e-05,
       4.84365929e-05, 2.93293649e-05, 2.24564621e-05, 1.19448771e-05,
       4.87111765e-06])

def expo(x, a, b):
    return exponent1 * a * np.exp(-b * x * exponent2)

# 提取拟合区间数据，自动过滤y为非正的异常点避免取对数报错
x_exp = bin_centers[start_exp:end_exp]
y_exp = norm_counts[start_exp:end_exp]
valid_mask = y_exp > 0
x_valid = x_exp[valid_mask]
y_valid = y_exp[valid_mask]

# 自动计算初始参数：线性回归ln(y)和x的关系
ln_y = np.log(y_valid)
slope, intercept, r_val, p_val, std_err = stats.linregress(x_valid, ln_y)
# 转换为拟合参数初始值
a0 = np.exp(intercept) / exponent1
b0 = -slope / exponent2

# 全自动拟合，无需手动输入p0
a_fit, b_fit = sc.curve_fit(expo, x_exp, y_exp, p0=(a0, b0))[0]

# 绘图
expon = expo(bin_centers[start_exp:end_exp], a_fit, b_fit)
plt.plot(bin_centers[start_exp:end_exp], expon,'-.',color='black',linewidth=0.8,label='${\propto} $ $e^{%.4f}$' %(-b_fit*exponent2))
plt.plot(bin_centers, norm_counts)
plt.yscale('log')
plt.xscale('log')
plt.legend()
plt.show()

高噪声场景下的稳健性优化 #

如果待拟合曲线存在大量异常点，可以将普通线性回归替换为稳健回归算法，进一步降低拟合失败概率：

• 替换stats.linregress为stats.theilslopes，该算法对异常值的耐受度远高于普通最小二乘，只需修改初值计算部分代码即可：

# 稳健线性回归计算初值
slope, intercept, _, _ = stats.theilslopes(ln_y, x_valid)
a0 = np.exp(intercept) / exponent1
b0 = -slope / exponent2

该方案完全无需人工干预，可直接嵌入批量处理脚本对大量曲线执行拟合操作，拟合成功率接近100%。

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