# Python曲线拟合 - 完整代码示例# 曲线拟合通过优化参数使模型函数最佳匹配观测数据import numpy as npimport matplotlib.pyplot as pltfrom scipy import optimize, stats# 1. 多项式拟合numpy.polyfit# 使用最小二乘法拟合多项式到数据点np.random.seed(42)x_data np.linspace(-3, 3, 30)# 真实模型 y 2x² - 3x 1 噪声y_true 2 * x_data**2 - 3 * x_data 1y_observed y_true 2.0 * np.random.randn(30) # 添加噪声# 拟合不同阶数的多项式deg2_coeffs np.polyfit(x_data, y_observed, 2) # 二次多项式deg5_coeffs np.polyfit(x_data, y_observed, 5) # 五次多项式可能过拟合# 使用 poly1d 从系数创建多项式函数p2 np.poly1d(deg2_coeffs)p5 np.poly1d(deg5_coeffs)print(多项式拟合结果:)print(f真实系数: [2, -3, 1])print(f二次拟合系数: {deg2_coeffs})print(f五次拟合系数: {deg5_coeffs})# 2. scipy.optimize.curve_fit 非线性曲线拟合# 定义要拟合的模型函数高斯函数def gaussian(x, A, mu, sigma, offset):高斯函数模型A * exp(-(x-mu)²/(2*sigma²)) offsetreturn A * np.exp(-(x - mu)**2 / (2 * sigma**2)) offset# 生成高斯型数据x_gauss np.linspace(-5, 5, 50)y_gauss_true 3.0 * np.exp(-(x_gauss - 0.5)**2 / (2 * 0.8**2)) 0.5y_gauss_obs y_gauss_true 0.3 * np.random.randn(50)# 使用 curve_fit 拟合需要提供初始猜测initial_guess [2.0, 0.0, 1.0, 0.0] # [A, mu, sigma, offset]popt, pcov optimize.curve_fit(gaussian, x_gauss, y_gauss_obs, p0initial_guess)# popt 是最优参数pcov 是协方差矩阵用于计算参数误差perr np.sqrt(np.diag(pcov)) # 参数的标准误差print(f\n高斯拟合结果:)print(f参数 真实值 拟合值 标准误差)params_names [A, mu, sigma, offset]params_true [3.0, 0.5, 0.8, 0.5]for i, name in enumerate(params_names):print(f{name:8s} {params_true[i]:8.3f} {popt[i]:8.3f} {perr[i]:8.4f})# 3. 拟合优度评估y_pred gaussian(x_gauss, *popt)residuals y_gauss_obs - y_predss_res np.sum(residuals**2) # 残差平方和ss_tot np.sum((y_gauss_obs - np.mean(y_gauss_obs))**2) # 总平方和r_squared 1 - ss_res / ss_tot # R² 决定系数# 计算调整 R² (考虑参数数量)n len(y_gauss_obs)k len(popt)adj_r_squared 1 - (1 - r_squared) * (n - 1) / (n - k - 1)print(f\n拟合优度:)print(fR² {r_squared:.4f})print(f调整 R² {adj_r_squared:.4f})print(f残差标准差 {np.std(residuals):.4f})# 4. 置信带和预测带# 使用协方差矩阵计算拟合值的置信区间alpha 0.05 # 显著性水平 5%t_value stats.t.ppf(1 - alpha/2, n - k - 1) # t 分布临界值# 构建设计矩阵模型在 x 处的雅可比矩阵J np.zeros((len(x_gauss), len(popt)))J[:, 0] np.exp(-(x_gauss - popt[1])**2 / (2 * popt[2]**2)) # d/dAJ[:, 1] popt[0] * np.exp(-(x_gauss - popt[1])**2 / (2 * popt[2]**2)) * (x_gauss - popt[1]) / popt[2]**2 # d/dmuJ[:, 2] popt[0] * np.exp(-(x_gauss - popt[1])**2 / (2 * popt[2]**2)) * (x_gauss - popt[1])**2 / popt[2]**3 # d/dsigmaJ[:, 3] np.ones(len(x_gauss)) # d/doffset# 拟合值的方差y_pred_var np.sum(J pcov * J, axis1)y_pred_se np.sqrt(y_pred_var) # 标准误差ci_half t_value * y_pred_se # 置信区间半宽# 预测区间包含观测误差pi_half t_value * np.sqrt(y_pred_var np.var(residuals))print(f\n在 x0.5 处的区间估计:)print(f拟合值: {gaussian(0.5, *popt):.4f})print(f95% 置信区间: [{gaussian(0.5, *popt) - ci_half[np.argmin(abs(x_gauss - 0.5))]:.4f}, f{gaussian(0.5, *popt) ci_half[np.argmin(abs(x_gauss - 0.5))]:.4f}])# 5. least_squares 更灵活的非线性最小二乘# 当 curve_fit 不够灵活时直接使用 least_squaresdef residual_func(params, x, y):残差函数返回每个数据点的残差A, mu, sigma, offset paramsy_model gaussian(x, A, mu, sigma, offset)return y - y_model # 返回残差向量# 使用边界约束进行拟合防止 sigma 为负bounds ([-10, -5, 0.01, -5], [10, 5, 5, 5]) # (下界, 上界)result_ls optimize.least_squares(residual_func, initial_guess,args(x_gauss, y_gauss_obs),boundsbounds,methodtrf # Trust Region Reflective)print(f\nleast_squares 拟合结果:)print(fA {result_ls.x[0]:.3f}, mu {result_ls.x[1]:.3f})print(fsigma {result_ls.x[2]:.3f}, offset {result_ls.x[3]:.3f})print(f优化成功: {result_ls.success})print(f迭代次数: {result_ls.nfev})# 6. 加权拟合当数据点具有不同权重时# 假设越靠近中心的数据越可靠权重越高weights np.exp(-x_gauss**2 / 10) # 权重函数def weighted_residual(params, x, y, w):加权残差函数return w * residual_func(params, x, y)result_weighted optimize.least_squares(lambda p, x, y: weighted_residual(p, x, y, weights),initial_guess, args(x_gauss, y_gauss_obs))print(f\n加权拟合 A {result_weighted.x[0]:.3f})# 7. 模型选择比较不同模型的 AIC/BIC# AIC n * ln(RSS/n) 2k, BIC n * ln(RSS/n) k*ln(n)for deg in [1, 2, 3, 4, 5]:coeffs np.polyfit(x_data, y_observed, deg)y_pred np.polyval(coeffs, x_data)rss np.sum((y_observed - y_pred)**2)k deg 1 # 参数个数aic n * np.log(rss / n) 2 * kbic n * np.log(rss / n) k * np.log(n)print(f阶数 {deg}: AIC{aic:.1f}, BIC{bic:.1f}, RSS{rss:.2f})print(\n曲线拟合总结选择合适模型复杂度避免过拟合)print(curve_fit 适合标准模型least_squares 提供更多控制)
Python曲线拟合
# Python曲线拟合 - 完整代码示例# 曲线拟合通过优化参数使模型函数最佳匹配观测数据import numpy as npimport matplotlib.pyplot as pltfrom scipy import optimize, stats# 1. 多项式拟合numpy.polyfit# 使用最小二乘法拟合多项式到数据点np.random.seed(42)x_data np.linspace(-3, 3, 30)# 真实模型 y 2x² - 3x 1 噪声y_true 2 * x_data**2 - 3 * x_data 1y_observed y_true 2.0 * np.random.randn(30) # 添加噪声# 拟合不同阶数的多项式deg2_coeffs np.polyfit(x_data, y_observed, 2) # 二次多项式deg5_coeffs np.polyfit(x_data, y_observed, 5) # 五次多项式可能过拟合# 使用 poly1d 从系数创建多项式函数p2 np.poly1d(deg2_coeffs)p5 np.poly1d(deg5_coeffs)print(多项式拟合结果:)print(f真实系数: [2, -3, 1])print(f二次拟合系数: {deg2_coeffs})print(f五次拟合系数: {deg5_coeffs})# 2. scipy.optimize.curve_fit 非线性曲线拟合# 定义要拟合的模型函数高斯函数def gaussian(x, A, mu, sigma, offset):高斯函数模型A * exp(-(x-mu)²/(2*sigma²)) offsetreturn A * np.exp(-(x - mu)**2 / (2 * sigma**2)) offset# 生成高斯型数据x_gauss np.linspace(-5, 5, 50)y_gauss_true 3.0 * np.exp(-(x_gauss - 0.5)**2 / (2 * 0.8**2)) 0.5y_gauss_obs y_gauss_true 0.3 * np.random.randn(50)# 使用 curve_fit 拟合需要提供初始猜测initial_guess [2.0, 0.0, 1.0, 0.0] # [A, mu, sigma, offset]popt, pcov optimize.curve_fit(gaussian, x_gauss, y_gauss_obs, p0initial_guess)# popt 是最优参数pcov 是协方差矩阵用于计算参数误差perr np.sqrt(np.diag(pcov)) # 参数的标准误差print(f\n高斯拟合结果:)print(f参数 真实值 拟合值 标准误差)params_names [A, mu, sigma, offset]params_true [3.0, 0.5, 0.8, 0.5]for i, name in enumerate(params_names):print(f{name:8s} {params_true[i]:8.3f} {popt[i]:8.3f} {perr[i]:8.4f})# 3. 拟合优度评估y_pred gaussian(x_gauss, *popt)residuals y_gauss_obs - y_predss_res np.sum(residuals**2) # 残差平方和ss_tot np.sum((y_gauss_obs - np.mean(y_gauss_obs))**2) # 总平方和r_squared 1 - ss_res / ss_tot # R² 决定系数# 计算调整 R² (考虑参数数量)n len(y_gauss_obs)k len(popt)adj_r_squared 1 - (1 - r_squared) * (n - 1) / (n - k - 1)print(f\n拟合优度:)print(fR² {r_squared:.4f})print(f调整 R² {adj_r_squared:.4f})print(f残差标准差 {np.std(residuals):.4f})# 4. 置信带和预测带# 使用协方差矩阵计算拟合值的置信区间alpha 0.05 # 显著性水平 5%t_value stats.t.ppf(1 - alpha/2, n - k - 1) # t 分布临界值# 构建设计矩阵模型在 x 处的雅可比矩阵J np.zeros((len(x_gauss), len(popt)))J[:, 0] np.exp(-(x_gauss - popt[1])**2 / (2 * popt[2]**2)) # d/dAJ[:, 1] popt[0] * np.exp(-(x_gauss - popt[1])**2 / (2 * popt[2]**2)) * (x_gauss - popt[1]) / popt[2]**2 # d/dmuJ[:, 2] popt[0] * np.exp(-(x_gauss - popt[1])**2 / (2 * popt[2]**2)) * (x_gauss - popt[1])**2 / popt[2]**3 # d/dsigmaJ[:, 3] np.ones(len(x_gauss)) # d/doffset# 拟合值的方差y_pred_var np.sum(J pcov * J, axis1)y_pred_se np.sqrt(y_pred_var) # 标准误差ci_half t_value * y_pred_se # 置信区间半宽# 预测区间包含观测误差pi_half t_value * np.sqrt(y_pred_var np.var(residuals))print(f\n在 x0.5 处的区间估计:)print(f拟合值: {gaussian(0.5, *popt):.4f})print(f95% 置信区间: [{gaussian(0.5, *popt) - ci_half[np.argmin(abs(x_gauss - 0.5))]:.4f}, f{gaussian(0.5, *popt) ci_half[np.argmin(abs(x_gauss - 0.5))]:.4f}])# 5. least_squares 更灵活的非线性最小二乘# 当 curve_fit 不够灵活时直接使用 least_squaresdef residual_func(params, x, y):残差函数返回每个数据点的残差A, mu, sigma, offset paramsy_model gaussian(x, A, mu, sigma, offset)return y - y_model # 返回残差向量# 使用边界约束进行拟合防止 sigma 为负bounds ([-10, -5, 0.01, -5], [10, 5, 5, 5]) # (下界, 上界)result_ls optimize.least_squares(residual_func, initial_guess,args(x_gauss, y_gauss_obs),boundsbounds,methodtrf # Trust Region Reflective)print(f\nleast_squares 拟合结果:)print(fA {result_ls.x[0]:.3f}, mu {result_ls.x[1]:.3f})print(fsigma {result_ls.x[2]:.3f}, offset {result_ls.x[3]:.3f})print(f优化成功: {result_ls.success})print(f迭代次数: {result_ls.nfev})# 6. 加权拟合当数据点具有不同权重时# 假设越靠近中心的数据越可靠权重越高weights np.exp(-x_gauss**2 / 10) # 权重函数def weighted_residual(params, x, y, w):加权残差函数return w * residual_func(params, x, y)result_weighted optimize.least_squares(lambda p, x, y: weighted_residual(p, x, y, weights),initial_guess, args(x_gauss, y_gauss_obs))print(f\n加权拟合 A {result_weighted.x[0]:.3f})# 7. 模型选择比较不同模型的 AIC/BIC# AIC n * ln(RSS/n) 2k, BIC n * ln(RSS/n) k*ln(n)for deg in [1, 2, 3, 4, 5]:coeffs np.polyfit(x_data, y_observed, deg)y_pred np.polyval(coeffs, x_data)rss np.sum((y_observed - y_pred)**2)k deg 1 # 参数个数aic n * np.log(rss / n) 2 * kbic n * np.log(rss / n) k * np.log(n)print(f阶数 {deg}: AIC{aic:.1f}, BIC{bic:.1f}, RSS{rss:.2f})print(\n曲线拟合总结选择合适模型复杂度避免过拟合)print(curve_fit 适合标准模型least_squares 提供更多控制)
