Scipy curve_fit与Excel岩崩幂律拟合结果差异大的原因咨询

问题描述 #

我采用幂律公式 y=aV^-b 拟合岩崩数据，使用Scipy Optimize的curve_fit函数得到的拟合结果与Excel差异显著，尤其在极端大值部分拟合效果很差。尽管R²值表现良好，但拟合曲线视觉上不符合预期。尝试截断小于100的数值后，问题仍存在。想咨询：

• 为何curve_fit拟合效果不佳？
• 是否应采用其他回归方法进行幂律拟合？

原始代码 #

x2 = Volume
y2 = cumulative_frequency

def rockfall(x,a1,b1):
    return a1*x**-b1

V = x2.values
quant = y2.values
c, cov = curve_fit(rockfall,V,quant)
print(c)

n = len(x2)
q = np.empty(n)
for i in range(n):
    q[i] = rockfall(x2[i],c[0],c[1])
    
from sklearn.metrics import r2_score
print('R^2: ',r2_score(quant,q))
print('Samples:',x2.size)
R2=r2_score(quant,q)

数据集 #

拟合效果差的原因 #

1. 最小二乘法的权重偏向小值点：curve_fit默认使用普通最小二乘法（OLS），最小化的是y的残差平方和。你的数据中，小体积样本占比极高（57个样本里仅8个体积大于1000），OLS会优先拟合数量多的小值点，大值点的残差影响被稀释，导致大值区域拟合偏差明显。
2. 与Excel拟合逻辑不同：Excel的幂律拟合是先对公式两边取对数，转换为线性模型ln(y) = ln(a) - b*ln(V)后做线性回归，本质是最小化ln(y)的残差平方和；而curve_fit直接拟合原幂律函数，最小化的是原始y的残差，两者的优化目标不一致，结果自然有差异。此外，大值点的y本身很小（如第一个点y=1），残差的权重天然更低，更容易被忽略。
3. 初始参数未优化：curve_fit默认初始参数为[1,1]，若数据与默认值差距较大，算法可能收敛到局部最优解，而非全局最优。

可行的改进方案 #

1. 改用对数线性回归（匹配Excel逻辑） #

先对数据取对数，转换为线性回归问题，再将参数转换回原幂律形式，结果会与Excel一致：

import numpy as np
from sklearn.linear_model import LinearRegression
import pandas as pd

# 构造数据集
data = pd.DataFrame({
    'V': [71197, 12594, 4780, 4578, 3590, 2624, 1699, 1025, 832.9, 654.74,
          572.55, 486.24, 391.46, 369.44, 361, 356.1, 314.91, 300, 294.7, 282.45,
          280.7, 276.6, 273, 258.96, 244.78, 223, 190, 189.7, 177.04, 176.16,
          175.8, 170.03, 168.23, 152.8, 141.6, 119.75, 102.76, 95.27, 90.16, 77.82,
          68.58, 59.92, 58.25, 49.44, 49.05, 42.9, 42.25, 39.55, 37.78, 36.31,
          30.84, 24.73, 23.2, 20.64, 18.67, 17.63, 11.13],
    'y': list(range(1, 58))
})

# 取对数转换
log_V = np.log(data['V'])
log_y = np.log(data['y'])

# 线性回归
model = LinearRegression()
model.fit(log_V.values.reshape(-1, 1), log_y)

# 转换为幂律参数
a = np.exp(model.intercept_)
b = -model.coef_[0]

print(f"幂律参数：a={a:.2f}, b={b:.2f}")

# 计算拟合结果与R²
y_pred = a * data['V'] ** (-b)
r2 = r2_score(data['y'], y_pred)
print(f"R²: {r2:.4f}")

2. 自定义加权损失函数 #

通过scipy.optimize.minimize自定义损失函数，给大体积点赋予更高权重，让拟合更关注大值区域：

from scipy.optimize import minimize

def weighted_loss(params, V, y):
    a, b = params
    y_pred = a * V ** (-b)
    # 以体积V作为权重，放大大值点的残差影响
    weights = V
    return np.sum(weights * (y - y_pred) ** 2)

# 用对数回归的结果作为初始参数，提升收敛效果
initial_guess = [a, b]
result = minimize(weighted_loss, initial_guess, args=(data['V'].values, data['y'].values))

a_weighted, b_weighted = result.x
print(f"加权拟合参数：a={a_weighted:.2f}, b={b_weighted:.2f}")

# 计算拟合结果与R²
y_pred_weighted = a_weighted * data['V'] ** (-b_weighted)
r2_weighted = r2_score(data['y'], y_pred_weighted)
print(f"加权拟合R²: {r2_weighted:.4f}")

3. 给curve_fit指定初始参数 #

如果坚持使用curve_fit，可以传入接近真实值的初始参数（比如Excel结果或对数回归结果），帮助算法找到最优解：

from scipy.optimize import curve_fit

# 用对数回归的a和b作为初始参数
c, cov = curve_fit(rockfall, data['V'].values, data['y'].values, p0=[a, b])
print(f"指定初始参数后的拟合结果：a={c[0]:.2f}, b={c[1]:.2f}")

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