Python:使用piecewise与curve_fit进行三段拟合
x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15,16,17,18,19,20,21], dtype=float) y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03,145,147,149,151,153,155]) plt.scatter(x,y,s=30,c='b')
得到如下散点图:
定义分段函数
#6个未知参数 x0x1,y0,y1分别是2个分割间断点的横纵坐标 k0,k1是第一和第三段直线的斜率 def piecewise(x,x0,x1,y0,y1,k0,k1): return np.piecewise(x , [x <= x0, np.logical_and(x0<x, x<= x1),x>x1] , [lambda x:k0*(x-x0) + y0,#根据点斜式构建函数 lambda x:(x-x0)*(y1-y0)/(x1-x0)+y0,#根据两点式构建函数 lambda x:k1*(x-x1) + y1])
根据分段函数进行拟合,通过迭代寻找最优的p,即为p_best
注:p(p_best)中包含的是拟合之后求得的所有未知参数
perr_min = np.inf p_best = None for n in range(100): k = np.random.rand(6)*20 p , e = optimize.curve_fit(piecewise, x, y,p0=k) perr = np.sum(np.abs(y-piecewise(x, *p))) if(perr < perr_min): perr_min = perr p_best = p
根据p_best调用curve_fit函数绘制拟合图像
xd = np.linspace(0, 21, 100) plt.figure() plt.plot(xd, piecewise(xd, *p_best)) xx=(p_best[0],p_best[1]) yy=(p_best[2],p_best[3]) plt.scatter(xx,yy,s=30,c='black') plt.show()
结果如下:
完整代码:
from scipy import optimize import matplotlib.pyplot as plt import numpy as np #6个未知参数 x0x1,y0,y1分别是2个分割间断点的横纵坐标 k0,k1是第一和第三段直线的斜率 def piecewise(x,x0,x1,y0,y1,k0,k1): return np.piecewise(x , [x <= x0, np.logical_and(x0<x, x<= x1),x>x1] , [lambda x:k0*(x-x0) + y0,#根据点斜式构建函数 lambda x:(x-x0)*(y1-y0)/(x1-x0)+y0,#根据两点式构建函数 lambda x:k1*(x-x1) + y1]) x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15,16,17,18,19,20,21], dtype=float) y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03,145,147,149,151,153,155]) plt.scatter(x,y,s=30,c='b') perr_min = np.inf p_best = None for n in range(100): k = np.random.rand(6)*20 p , e = optimize.curve_fit(piecewise, x, y,p0=k) perr = np.sum(np.abs(y-piecewise(x, *p))) if(perr < perr_min): perr_min = perr p_best = p xd = np.linspace(0, 21, 100) plt.figure() plt.plot(xd, piecewise(xd, *p_best)) xx=(p_best[0],p_best[1]) yy=(p_best[2],p_best[3]) plt.scatter(xx,yy,s=30,c='black') plt.show()
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 浏览器原生「磁吸」效果!Anchor Positioning 锚点定位神器解析
· 没有源码,如何修改代码逻辑?
· 一个奇形怪状的面试题:Bean中的CHM要不要加volatile?
· [.NET]调用本地 Deepseek 模型
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· PowerShell开发游戏 · 打蜜蜂
· 在鹅厂做java开发是什么体验
· 百万级群聊的设计实践
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
· 永远不要相信用户的输入:从 SQL 注入攻防看输入验证的重要性