本文目录导读:

我来详细介绍如何使用Scipy进行曲线拟合,包含多个实际案例。
基础安装和导入
import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit from scipy import stats plt.rcParams['font.sans-serif'] = ['SimHei'] # 支持中文显示 plt.rcParams['axes.unicode_minus'] = False
线性拟合案例
# 生成示例数据
np.random.seed(42)
x_data = np.linspace(0, 10, 50)
# y = 2x + 1 + 噪声
y_data = 2 * x_data + 1 + np.random.normal(0, 0.5, len(x_data))
# 定义线性函数
def linear_func(x, a, b):
return a * x + b
# 进行曲线拟合
params, params_covariance = curve_fit(linear_func, x_data, y_data)
a_fit, b_fit = params
# 计算拟合优度
residuals = y_data - linear_func(x_data, *params)
ss_res = np.sum(residuals**2)
ss_tot = np.sum((y_data - np.mean(y_data))**2)
r_squared = 1 - (ss_res / ss_tot)
print(f"拟合参数: a = {a_fit:.3f}, b = {b_fit:.3f}")
print(f"R² = {r_squared:.3f}")
# 可视化
plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, label='原始数据', alpha=0.6)
plt.plot(x_data, linear_func(x_data, *params), 'r-', label=f'线性拟合: y={a_fit:.2f}x+{b_fit:.2f}')
plt.xlabel('x')
plt.ylabel('y')'线性曲线拟合')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
多项式拟合案例
# 生成二次函数数据
x_data = np.linspace(-5, 5, 100)
# y = x² - 2x + 3 + 噪声
y_data = x_data**2 - 2*x_data + 3 + np.random.normal(0, 2, len(x_data))
# 定义二次函数
def quadratic_func(x, a, b, c):
return a * x**2 + b * x + c
# 拟合
params, _ = curve_fit(quadratic_func, x_data, y_data)
a_fit, b_fit, c_fit = params
print(f"二次拟合: y = {a_fit:.2f}x² + {b_fit:.2f}x + {c_fit:.2f}")
# 可视化
x_smooth = np.linspace(-5, 5, 200)
y_smooth = quadratic_func(x_smooth, *params)
plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, label='原始数据', alpha=0.5, s=30)
plt.plot(x_smooth, y_smooth, 'r-', label='二次拟合曲线', linewidth=2)
plt.xlabel('x')
plt.ylabel('y')'二次多项式拟合')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
指数衰减拟合案例
# 生成指数衰减数据
t_data = np.linspace(0, 10, 50)
# y = 5 * exp(-0.5x) + 噪声
y_data = 5 * np.exp(-0.5 * t_data) + np.random.normal(0, 0.2, len(t_data))
# 定义指数函数
def exp_decay(t, A, k, c):
return A * np.exp(-k * t) + c
# 拟合(提供初始猜测值)
initial_guess = [1, 0.1, 1]
params, params_covariance = curve_fit(exp_decay, t_data, y_data, p0=initial_guess)
A_fit, k_fit, c_fit = params
print(f"指数拟合: y = {A_fit:.2f} * exp(-{k_fit:.3f}x) + {c_fit:.2f}")
# 可视化
t_smooth = np.linspace(0, 10, 200)
y_smooth = exp_decay(t_smooth, *params)
plt.figure(figsize=(10, 6))
plt.scatter(t_data, y_data, label='原始数据', alpha=0.6)
plt.plot(t_smooth, y_smooth, 'r-', label='指数衰减拟合', linewidth=2)
plt.xlabel('时间')
plt.ylabel('幅值')'指数衰减拟合')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
正弦波拟合案例
# 生成正弦波数据
np.random.seed(123)
x_data = np.linspace(0, 10, 100)
# y = 3 * sin(2x + π/4) + 1 + 噪声
y_data = 3 * np.sin(2 * x_data + np.pi/4) + 1 + np.random.normal(0, 0.2, len(x_data))
# 定义正弦函数
def sine_func(x, amplitude, frequency, phase, offset):
return amplitude * np.sin(frequency * x + phase) + offset
# 拟合(需要好的初始猜测)
initial_guess = [2, 2, 0, 0] # [振幅, 频率, 相位, 偏移]
params, _ = curve_fit(sine_func, x_data, y_data, p0=initial_guess)
A_fit, f_fit, p_fit, o_fit = params
print(f"正弦拟合: y = {A_fit:.2f} * sin({f_fit:.2f}x + {p_fit:.2f}) + {o_fit:.2f}")
# 可视化
x_smooth = np.linspace(0, 10, 200)
y_smooth = sine_func(x_smooth, *params)
plt.figure(figsize=(12, 6))
plt.scatter(x_data, y_data, label='原始数据', alpha=0.5, s=30)
plt.plot(x_smooth, y_smooth, 'r-', label='正弦拟合', linewidth=2)
plt.xlabel('x')
plt.ylabel('y')'正弦波拟合')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
高斯函数拟合案例
# 生成高斯分布数据
x_data = np.linspace(-5, 5, 100)
# 真实参数:均值=0,标准差=1,幅值=2
y_true = 2 * np.exp(-(x_data**2) / (2 * 1**2))
y_data = y_true + np.random.normal(0, 0.1, len(x_data))
# 定义高斯函数
def gaussian(x, amplitude, mean, std_dev):
return amplitude * np.exp(-(x - mean)**2 / (2 * std_dev**2))
# 拟合
initial_guess = [1, 0, 1] # [幅值, 均值, 标准差]
params, _ = curve_fit(gaussian, x_data, y_data, p0=initial_guess)
A_fit, mu_fit, sigma_fit = params
print(f"高斯拟合: amplitude={A_fit:.2f}, mean={mu_fit:.2f}, std={sigma_fit:.2f}")
# 可视化
x_smooth = np.linspace(-5, 5, 200)
y_smooth = gaussian(x_smooth, *params)
plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, label='原始数据', alpha=0.6, s=30)
plt.plot(x_smooth, y_smooth, 'r-', label='高斯拟合', linewidth=2)
plt.xlabel('x')
plt.ylabel('y')'高斯函数拟合')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
自定义复杂函数拟合
# 生成洛伦兹函数数据
def lorentzian(x, amplitude, center, gamma):
return amplitude * gamma**2 / ((x - center)**2 + gamma**2)
x_data = np.linspace(-10, 10, 100)
y_data = lorentzian(x_data, 10, 0, 2) + np.random.normal(0, 0.3, len(x_data))
# 拟合
initial_guess = [8, 0, 2]
params, params_covariance = curve_fit(lorentzian, x_data, y_data, p0=initial_guess)
# 计算参数误差(标准误)
perr = np.sqrt(np.diag(params_covariance))
print("洛伦兹拟合结果:")
print(f"幅值: {params[0]:.2f} ± {perr[0]:.2f}")
print(f"中心: {params[1]:.2f} ± {perr[1]:.2f}")
print(f"半高宽: {params[2]:.2f} ± {perr[2]:.2f}")
# 可视化
x_smooth = np.linspace(-10, 10, 200)
y_smooth = lorentzian(x_smooth, *params)
plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, label='原始数据', alpha=0.6, s=30)
plt.plot(x_smooth, y_smooth, 'r-', label='洛伦兹拟合', linewidth=2)
plt.fill_between(x_smooth, y_smooth - np.std(y_data - lorentzian(x_data, *params)),
y_smooth + np.std(y_data - lorentzian(x_data, *params)), alpha=0.2)
plt.xlabel('x')
plt.ylabel('y')'洛伦兹函数拟合(含误差带)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
实用工具函数
def calculate_fit_statistics(y_true, y_pred, n_params):
"""
计算拟合统计量
"""
residuals = y_true - y_pred
n = len(y_true)
# 残差平方和
rss = np.sum(residuals**2)
# 总平方和
tss = np.sum((y_true - np.mean(y_true))**2)
# R²
r_squared = 1 - rss / tss
# 调整R²
adj_r_squared = 1 - (rss / (n - n_params)) / (tss / (n - 1))
# 均方根误差 (RMSE)
rmse = np.sqrt(rss / n)
# 平均绝对误差 (MAE)
mae = np.mean(np.abs(residuals))
return {
'R²': r_squared,
'调整R²': adj_r_squared,
'RMSE': rmse,
'MAE': mae,
'RSS': rss
}
def batch_fitting_example():
"""
批量拟合示例:显示多种拟合方法的比较
"""
np.random.seed(42)
x = np.linspace(0, 10, 50)
y = 2 * np.exp(-0.3 * x) + 1 + np.random.normal(0, 0.2, len(x))
# 定义多种拟合函数
fit_functions = {
'Linear': lambda x, a, b: a * x + b,
'Exp': lambda x, a, b, c: a * np.exp(-b * x) + c,
'Quadratic': lambda x, a, b, c: a * x**2 + b * x + c,
'Power': lambda x, a, b: a * x**b
}
results = {}
x_smooth = np.linspace(0, 10, 200)
plt.figure(figsize=(12, 8))
plt.scatter(x, y, label='原始数据', alpha=0.6)
for name, func in fit_functions.items():
try:
# 初始猜测
if name == 'Exp':
p0 = [1, 0.1, 0]
else:
p0 = [1, 1, 1] if 'c' in func.__code__.co_varnames[:3] else [1, 1]
params, _ = curve_fit(func, x, y, p0=p0, maxfev=5000)
y_pred = func(x, *params)
stats = calculate_fit_statistics(y, y_pred, len(params))
results[name] = {
'params': params,
'stats': stats
}
# 绘制拟合曲线
if x_smooth is not None:
y_smooth = func(x_smooth, *params)
plt.plot(x_smooth, y_smooth, '-', label=f'{name} (R²={stats["R²"]:.3f})',
linewidth=2)
print(f"\n{name}拟合:")
print(f" 参数: {params}")
print(f" R²: {stats['R²']:.4f}, RMSE: {stats['RMSE']:.4f}")
except Exception as e:
print(f"\n{name}拟合失败: {e}")
plt.xlabel('x')
plt.ylabel('y')
plt.title('不同拟合方法比较')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
return results
# 运行批量拟合
# results = batch_fitting_example()
关键技术点
初始猜测的重要性
- 对于复杂函数,提供好的初始猜测可以避免收敛到局部最优
- 使用
p0参数设置初始猜测
拟合评估
# 获取拟合协方差矩阵(用于计算参数误差) params, params_covariance = curve_fit(func, x, y) perr = np.sqrt(np.diag(params_covariance)) # 参数标准误 # 计算置信区间 confidence = 0.95 n = len(y) p = len(params) dof = max(0, n - p) # 自由度 t_value = stats.t.ppf((1 + confidence) / 2, dof)
权重拟合
# 当数据点有不同的不确定性时 weights = 1 / (sigma_y**2) # 权重 params, _ = curve_fit(func, x, y, sigma=sigma_y, absolute_sigma=True)
边界约束
# 设置参数边界 bounds = ([0, -np.inf, 0], [np.inf, np.inf, np.inf]) # (下界, 上界) params, _ = curve_fit(func, x, y, bounds=bounds)
这些案例涵盖了大多数常见的数据拟合需求,可以根据实际应用场景选择合适的拟合函数和方法。