文档章节

scipy插值与拟合

o
 osc_fmg49rzg
发布于 2019/03/20 10:36
字数 2496
阅读 55
收藏 0

精选30+云产品,助力企业轻松上云!>>>

原文链接:https://zhuanlan.zhihu.com/p/28149195

 

1.最小二乘拟合

实例1

import numpy as np
import matplotlib.pyplot as plt from scipy.optimize import leastsq plt.figure(figsize=(9,9)) x=np.linspace(0,10,1000) X = np.array([8.19, 2.72, 6.39, 8.71, 4.7, 2.66, 3.78]) Y = np.array([7.01, 2.78, 6.47, 6.71, 4.1, 4.23, 4.05]) #计算以p为参数的直线和原始数据之间的误差 def f(p): k, b = p return(Y-(k*X+b)) #leastsq使得f的输出数组的平方和最小,参数初始值为[1,0] r = leastsq(f, [1,0]) k, b = r[0] print("k=",k,"b=",b) plt.scatter(X,Y, s=100, alpha=1.0, marker='o',label=u'数据点') y=k*x+b ax = plt.gca() #gca获取轴这个对象 ax.set_xlabel(..., fontsize=20) ax.set_ylabel(..., fontsize=20) #设置坐标轴标签字体大小 plt.plot(x, y, color='r',linewidth=5, linestyle=":",markersize=20, label=u'拟合曲线') plt.legend(loc=0, numpoints=1) leg = plt.gca().get_legend() ltext = leg.get_texts() plt.setp(ltext, fontsize='xx-large') plt.xlabel(u'安培/A') plt.ylabel(u'伏特/V') plt.xlim(0, x.max() * 1.1) plt.ylim(0, y.max() * 1.1) plt.xticks(fontsize=20) plt.yticks(fontsize=20) #刻度字体大小 plt.legend(loc='upper left') plt.show() 

 

 

实例2

#最小二乘拟合实例
import numpy as np from scipy.optimize import leastsq import pylab as pl def func(x, p): """  数据拟合所用的函数: A*cos(2*pi*k*x + theta)  """ A, k, theta = p return A*np.sin(k*x+theta) def residuals(p, y, x): """  实验数据x, y和拟合函数之间的差,p为拟合需要找到的系数  """ return y - func(x, p) x = np.linspace(0, 20, 100) A, k, theta = 10, 3, 6 # 真实数据的函数参数 y0 = func(x, [A, k, theta]) # 真实数据 y1 = y0 + 2 * np.random.randn(len(x)) # 加入噪声之后的实验数据 p0 = [10, 0.2, 0] # 第一次猜测的函数拟合参数 # 调用leastsq进行数据拟合 # residuals为计算误差的函数 # p0为拟合参数的初始值 # args为需要拟合的实验数据 plsq = leastsq(residuals, p0, args=(y1, x)) print (u"真实参数:", [A, k, theta] ) print (u"拟合参数", plsq[0]) # 实验数据拟合后的参数 pl.plot(x, y0, color='r',label=u"真实数据") pl.plot(x, y1, color='b',label=u"带噪声的实验数据") pl.plot(x, func(x, plsq[0]), color='g', label=u"拟合数据") pl.legend() pl.show() 

 

 

 

2. 插值

# -*- coding: utf-8 -*-
"""
Created on Thu Jul 27 16:42:30 2017

@author: Dell
""" import numpy as np import pylab as pl from scipy import interpolate import matplotlib.pyplot as plt x = np.linspace(0, 2*np.pi+np.pi/4, 10) y = np.sin(x) x_new = np.linspace(0, 2*np.pi+np.pi/4, 100) f_linear = interpolate.interp1d(x, y) tck = interpolate.splrep(x, y) y_bspline = interpolate.splev(x_new, tck) plt.xlabel(u'安培/A') plt.ylabel(u'伏特/V') plt.plot(x, y, "o", label=u"原始数据") plt.plot(x_new, f_linear(x_new), label=u"线性插值") plt.plot(x_new, y_bspline, label=u"B-spline插值") pl.legend() pl.show() 

 

 

 

 

实例分析

# -*- coding: utf-8 -*-
"""
Created on Thu Jul 27 16:53:21 2017

@author: Dell
""" import numpy as np from scipy import interpolate import pylab as pl #创建数据点集并绘制 pl.figure(figsize=(12,9)) x = np.linspace(0, 10, 11) y = np.sin(x) ax=pl.plot() pl.plot(x,y,'ro') #建立插值数据点 xnew = np.linspace(0, 10, 101) for kind in ['nearest', 'zero','linear','quadratic']: #根据kind创建插值对象interp1d f = interpolate.interp1d(x, y, kind = kind) ynew = f(xnew)#计算插值结果 pl.plot(xnew, ynew, label = str(kind)) pl.xticks(fontsize=20) pl.yticks(fontsize=20) pl.legend(loc = 'lower right') pl.show() 

 

 

 

 

B样条曲线插值
一维数据的插值运算可以通过 interp1d()实现。
其调用形式为:
Interp1d可以计算x的取值范围之内任意点的函数值,并返回新的数组。
interp1d(x, y, kind=‘linear’, …)
参数 x和y是一系列已知的数据点
参数kind是插值类型,可以是字符串或整数

B样条曲线插值
Kind给出了B样条曲线的阶数:
 ‘
zero‘ ‘nearest’ :0阶梯插值,相当于0阶B样条曲线
 ‘slinear’‘linear’ :线性插值,相当于1阶B样条曲线
 ‘quadratic’‘cubic’:2阶和3阶B样条曲线,更高阶的曲线可以直接使用整数值来指定

(1)#创建数据点集:

import numpy as np

x = np.linspace(0, 10, 11)

y = np.sin(x)

 

(2)#绘制数据点集:

import pylab as pl

pl.plot(x,y,'ro')

 

 

 

 

创建interp1d对象f、计算插值结果:
xnew = np.linspace(0, 10, 11)
from scipy import interpolate
f = interpolate.interp1d(x, y, kind = kind)
ynew = f(xnew)

根据kind类型创建interp1d对象f、计算并绘制插值结果:
xnew = np.linspace(0, 10, 11)
for kind in ['nearest', 'zero','linear','quadratic']:
#根据kind创建插值对象interp1d
f = interpolate.interp1d(x, y, kind = kind)
ynew = f(xnew)#计算插值结果
pl.plot(xnew, ynew, label = str(kind))#绘制插值结果

 

如果我们将代码稍作修改增加一个5阶插值

# -*- coding: utf-8 -*-
"""
Created on Thu Jul 27 16:53:21 2017

@author: Dell
""" import numpy as np from scipy import interpolate import pylab as pl #创建数据点集并绘制 pl.figure(figsize=(12,9)) x = np.linspace(0, 10, 11) y = np.sin(x) ax=pl.plot() pl.plot(x,y,'ro') #建立插值数据点 xnew = np.linspace(0, 10, 101) for kind in ['nearest', 'zero','linear','quadratic',5]: #根据kind创建插值对象interp1d f = interpolate.interp1d(x, y, kind = kind) ynew = f(xnew)#计算插值结果 pl.plot(xnew, ynew, label = str(kind)) pl.xticks(fontsize=20) pl.yticks(fontsize=20) pl.legend(loc = 'lower right') pl.show() 运行得到 

 

发现5阶已经很接近正弦曲线,但是如果x值选取范围较大,则会出现跳跃。

 

关于拟合与插值的数学基础可参见霍开拓:拟合与插值的区别?

 

左边插值,右边拟合

仔细看有啥不一样

插值曲线要过数据点,拟合曲线整体效果更好。

插值,对准了才可以插吗,那就一定得过数据点。拟合,就是要得到最接近的结果,是要看总体效果。

 

既然理想(思路)不一样,那么三观和行为(特点和策略)也就不一样啦。

插值是指已知某函数的在若干离散点上的函数值或者导数信息,通过求解该函数中待定形式的插值函数以及待定系数,使得该函数在给定离散点上满足约束。

所谓拟合是指已知某函数的若干离散函数值{f1,f2,…,fn},通过调整该函数中若干待定系数f(λ1, λ2,…,λn), 使得该函数与已知点集的差别(最小二乘意义)最小。如果待定函数是线性,就叫线性拟合或者线性回归(主要在统计中),否则叫作非线性拟合或者非线性回归。表达式也可以是分段函数,这种情况下叫作样条拟合。
从几何意义上将,拟合是给定了空间中的一些点,找到一个已知形式未知参数的连续曲面来最大限度地逼近这些点;而插值是找到一个( 或几个分片光滑的)连续曲面来穿过这些点。

 

插值法

以下引自某科

Lagrange插值
Lagrange插值是n次多项式插值,其成功地用构造插值基函数的 方法解决了求n次多项式插值函数问题。
★基本思想 将待求的n次多项式插值函数pn(x)改写成另一种表示方式,再利用插值条件⑴确定其中的待定函数,从而求出插值多项式。
Newton插值
Newton插值也是n次多项式插值,它提出另一种构造插值多项式的方法,与Lagrange插值相比,具有承袭性和易于变动节点的特点。
★基本思想 将待求的n次插值多项式Pn(x)改写为具有承袭性的形式,然后利用插值条件⑴确定Pn(x)的待定系数,以求出所要的插值函数。
Hermite插值
Hermite插值是利用未知函数f(x)在插值节点上的函数值及 导数值来构造插值多项式的,其提法为:给定n+1个互异的节点x0,x1,……,xn上的函数值和导数值求一个2n+1次多项式H2n+1(x)满足插值条件H2n+1(xk)=ykH'2n+1(xk)=y'k k=0,1,2,……,n ⒀如上求出的H2n+1(x)称为2n+1次Hermite插值函数,它与被插函数一般有更好的密合度.
★基本思想利用Lagrange插值函数的构造方法,先设定函数形式,再利用插值条件⒀求出插值函数.

貌似插值节点取的越多,差值曲线或曲面越接近原始曲线/曲面,因为采样多嘛。但事实总是不像广大人民群众想的那样,随着插值节点的增多,多项式次数也在增高,插值曲线在一些区域出现跳跃,并且越来越偏离原始曲线。这个现象被 Tolmé Runge 发现并解释,然后就以他的名字命名这种现象。It was discovered by Carl David Tolmé Runge (1901) when exploring the behavior of errors when using polynomial interpolation to approximate certain functions.

为了解决这个问题,人们发明了分段插值法。分段插值一般不会使用四次以上的多项式,而二次多项式会出现尖点,也是有问题的。所以就剩下线性和三次插值,最后使用最多的还是线性分段插值,这个好处是显而易见的。

 

拟合

最小二乘

如何找到最接近原始曲线或者数据点的拟合曲线,这不是一件容易操作的事。要想整体最接近,直接的想法就是拟合曲线的每一点到原始曲线的对应点的最接近,简单点说就是两曲线上所有点的函数值之差的绝对值之和最小。看似解决问题,但绝对值在数学上向来是个不好交流的语言障碍患者,那然后又该怎么办。数学家说了既然办不了你绝对值之和,那就办了你家亲戚,就看你平方之和长得像。于是就找了这个长得像的来背黑锅,大家都表示很和谐。然后给这种操作冠之名曰'最小二乘法'。

官方一点的表述 , 选择参数c使得拟合模型与实际观测值在曲线拟合各点的残差(或离差)ek=yk-f(xk,c)的加权平方和达到最小,此时所求曲线称作在加权最小二乘意义下对数据的拟合曲线,这种方法叫做最小二乘法。

o
粉丝 0
博文 500
码字总数 0
作品 0
私信 提问
加载中
请先登录后再评论。
python scipy样条插值函数大全(interpolate里interpld函数)

scipy样条插值 scipy样条插值 1、样条插值法是一种以可变样条来作出一条经过一系列点的光滑曲线的数学方法。插值样条是由一些多项式组成的,每一个多项式都是由相邻的两个数据点决定的,这样...

osc_evqmj7fy
2019/07/15
32
0
SciPy 优化

<div class="article-child "><h2>章节</h2><ul><li class="page_item page-item-3474"><a href="https://www.qikegu.com/docs/3474">SciPy 介绍</a></li><li class="page_item page-item-3......

osc_6pkt76kw
2019/11/08
13
0
杏彩平台源码下载带龙虎和玩法

杏彩平台源码下载 带龙虎和玩法 地址一:【hubawl.com】 地址二:【bbscherry.com】 scipy包含致力于科学计算中常见问题的各个工具箱。它的不同子模块相应于不同的应用。像插值,积分,优化,...

asdasd12
2018/06/29
0
0
SciPy 安装

<div class="article-child "><h2>章节</h2><ul><li class="page_item page-item-3474"><a href="https://www.qikegu.com/docs/3474">SciPy 介绍</a></li><li class="page_item page-item-3......

osc_7hoa7os1
2019/11/05
8
0
Scipy的应用

首先总体概括一下Scipy的用处 >>> #Scipy依赖于numpy >>> #Scipy提供了真正的矩阵 >>> #Scipy包含的功能:最优化,线性代数,积分,插值,拟合,特殊函数,快速傅里叶变换,信号处理,图形处...

osc_g419uyhd
2018/07/26
7
0

没有更多内容

加载失败,请刷新页面

加载更多

MongoDB入门系列——3.可视化工具篇

点击上方,轻松关注!! 前面我们已经介绍了MongoDB怎么安装,接下来要安装他的可视化工具——Studio 3T。 先到这下载一个压缩包,百度网盘,https://pan.baidu.com/s/1M8mlWo334KE8I1_UA2Da...

学习Java的小姐姐
2018/11/08
0
0
分层图的绘制 python(来自国外课程)

Exercise 10: Hierarchical clustering of the grain data In the video, you learnt that the SciPy linkage() function performs hierarchical clustering on an array of samples. Use th......

齐勇cn
27分钟前
13
0
微信小程序超简单的双向绑定(类似vue的v-model)

<input model:value="{{value}}" />

祖达
27分钟前
9
0
为什么AngularJS在select中包含一个空选项? - Why does AngularJS include an empty option in select?

问题: I've been working with AngularJS for the last few weeks, and the one thing which is really bothering me is that even after trying all permutations or the configuration de......

技术盛宴
30分钟前
13
0
centos宝塔面板安装及常见错误处理(超级详细)

原文连接:https://www.wjcms.net/archives/centos%E5%AE%9D%E5%A1%94%E9%9D%A2%E6%9D%BF%E5%AE%89%E8%A3%85%E5%8F%8A%E5%B8%B8%E8%A7%81%E9%94%99%E8%AF%AF%E5%A4%84%E7%90%86%E8%B6%85%E7%......

神兵小将
51分钟前
17
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部