Python案例如何用Scipy做曲线拟合

wen python案例 2

本文目录导读:

Python案例如何用Scipy做曲线拟合

  1. 基础安装和导入
  2. 线性拟合案例
  3. 多项式拟合案例
  4. 指数衰减拟合案例
  5. 正弦波拟合案例
  6. 高斯函数拟合案例
  7. 自定义复杂函数拟合
  8. 实用工具函数
  9. 关键技术点

我来详细介绍如何使用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)

这些案例涵盖了大多数常见的数据拟合需求,可以根据实际应用场景选择合适的拟合函数和方法。

抱歉,评论功能暂时关闭!