Python中如何对数据进行函数拟合并优化简洁模型拟合效果

你当前用的1/(a*x² + b*x + c)模型拟合效果差，核心有两个原因：

• 初始参数设置完全不合理：你的X取值范围是4~112，x²项量级可达1e4，初始值设为[1,1,1]会让二次项在拟合初始阶段权重过大，curve_fit很容易陷入局部最优。
• 模型形式和数据趋势不匹配：从你的数据看，y随x增大单调递减，但下降速率持续放缓，x到112时y仍接近1，没有快速趋近于0的趋势；而二次倒数模型当x→∞时y会快速趋近于0，和数据长期趋势不符。

方案1：优化现有二次倒数模型（不换模型形式） #

不需要加高阶项，只要给curve_fit传递合理的初始参数即可，不需要瞎猜初始值：先对y取倒数得到z=1/y，对z和X做二次多项式线性拟合，得到的系数直接作为初始参数，线性拟合得到的参数已经非常接近全局最优，不会出现局部收敛问题。
代码示例：

import numpy as np
from scipy.optimize import curve_fit

# 你的原始数据
y = np.array([1.45952016, 1.36947283, 1.31433227, 1.24076599, 1.20577963,
       1.14454815, 1.13068077, 1.09638278, 1.08121406, 1.04417094,
       1.02251471, 1.01268524, 0.98535659, 0.97400591])
X = np.array([4.571428571362048, 8.771428571548313, 12.404761904850602, 17.904761904850602,
            22.904761904850602, 31.238095237873495, 37.95833333302289, 
            44.67857142863795, 51.39880952378735, 64.83928571408615, 
            71.5595238097012, 85., 98.55357142863795, 112.1071428572759])

def func_quad_inv(x, a, b, c):
    return 1/(a*(x**2) + b*x + c)

# 线性预拟合得到初始参数
z = 1 / y
p0 = np.polyfit(X, z, 2)  # 输出顺序是a,b,c，和func参数顺序一致
c_opt, cov = curve_fit(func_quad_inv, X, y, p0=p0)

# 生成预测值可以不用循环，直接向量化计算
test_ar = np.arange(min(X), max(X), 0.25)
pred = func_quad_inv(test_ar, *c_opt)

这个优化后二次倒数模型的R²可以到0.99以上，但要注意这个模型外推x>150之后y会快速下降，不符合数据的平缓趋势。如果你想要更简洁的形式，甚至可以直接去掉二次项，用两参数模型y=1/(a*x + b)，同样用z=1/y做线性拟合得到初始值，在当前数据集上R²也能达到0.99以上，完全满足精度要求。

方案2：换用更适配趋势的三参数简洁模型 #

如果要兼顾简洁性（和你当前模型参数数量一致，都是3个参数）和趋势合理性，优先选以下两种形式：

（1）带渐近值的幂函数模型 #

模型形式：y = a * x**(-b) + c

• 物理意义明确：c是x趋近无穷大时y的渐近值，b控制衰减速率，a是缩放系数；
• 完美适配y递减越来越慢、最终趋于平稳的趋势，外推合理性远好于二次倒数模型；
• 初始值可以简单设为p0=[1, 0.1, 0.9]（c接近你数据末尾的y值即可）。
  代码示例：

def func_power(x, a, b, c):
    return a * np.power(x, -b) + c

p0_power = [2, 0.2, 0.9]
c_power, cov_power = curve_fit(func_power, X, y, p0=p0_power, bounds=([0,0,0], [10,1,2]))
pred_power = func_power(test_ar, *c_power)

这个模型在你数据集上的R²可以到0.995以上，比优化后的二次倒数模型精度更高，外推也更合理。

（2）可变幂次的倒数模型 #

模型形式：y = 1/(a * x**k + b)

• 相当于把你现有模型里固定的x²的幂次改成可学习的k，也是3个参数，形式和你原有模型接近，灵活度更高；
• 拟合得到的k一般会远小于2（大概在0.3~0.5区间），不会出现x大的时候y快速掉0的问题。
  代码示例：

def func_inv_power(x, a, k, b):
    return 1/(a * np.power(x, k) + b)

p0_inv = [0.01, 0.5, 0.6]
c_inv, cov_inv = curve_fit(func_inv_power, X, y, p0=p0_inv)
pred_inv = func_inv_power(test_ar, *c_inv)

拟合效果验证 #

你可以用以下代码计算不同模型的R²，对比精度：

def r2(y_true, y_pred):
    ss_res = np.sum((y_true - y_pred)**2)
    ss_tot = np.sum((y_true - np.mean(y_true))**2)
    return 1 - ss_res/ss_tot

print("优化后二次倒数模型R²：", r2(y, func_quad_inv(X, *c_opt)))
print("幂函数模型R²：", r2(y, func_power(X, *c_power)))
print("可变幂次倒数模型R²：", r2(y, func_inv_power(X, *c_inv)))

补充注意：如果你的数据有明确的物理意义，优先选符合物理规律的模型形式，不要只看拟合R²。

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