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()

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)的加权平方和达到最小,此时所求曲线称作在加权最小二乘意义下对数据的拟合曲线,这种方法叫做最小二乘法。

数值运算pythonmopn_Python SciPy库——拟合与插值相关推荐

  1. java的拟合与插值库,插值和拟合区别

    插值和拟合都是函数逼近或者数值逼近的重要组成部分.他们的共同点都是通过已知一些离散点集M上的约束,求取一个定义在连续集合S(M包含于S)的未知连续函数,从而达到获取整体规律目的,即通过"窥几 ...

  2. 使用计算机进行数值运算 可根据需要达到,计算机统考题库.pdf

    一.关于计算机的诞生与发展 1. 一般认为,世界上第一台电子数字计算机诞生于A.1946年 2. 当前的计算机一般被认为是第四代计算机,它所采用的逻辑元 件是C.大规模集成电路 3. 下列关于世界上第 ...

  3. Python库之Scipy库的简介、安装详细

    目录 Scipy库的简介 Scipy库的安装 Scipy库的简介 Scipy高级科学计算库:和Numpy联系很密切,Scipy一般都是操控Numpy数组来进行科学计算.统计分析,所以可以说是基于Num ...

  4. Python数据分析入门--SciPy库学习笔记

    文章目录 前言 Scipy库简单入门 1.cluster模块 2. constants模块 3. fftpack模块 4. integrate 模块 5. interpolate 模块 6. lina ...

  5. matlab实验二数值运算报告,MATLAB数值运算实验报告.docx

    MATLAB数值运算实验报告 实验报告系 (部): 信息工程 班 级: 姓 名: 学 号: 课 程: MATLAB 实验名称: Matlab数值运算目录一 . 实验目的2二 . 实验内容2三 . 实验 ...

  6. python scipy库总结

    本节导图:https://www.processon.com/view/link/5fd9f5dc5653bb06f344d655 大纲 文章目录 大纲 1. 数值计算库SciPy 2. 拟合与优化 ...

  7. python解非线性规划问题讲析_python中线性规划中的单纯形法、scipy库与非线性规划求解问题...

    单纯形法.scipy库与非线性规划求解问题 单纯形法的基本定义 大M法求解线性规划的原理 excel求解 Python调用optimize包和scipy求解线性规划 Python编程实现单纯形法 对比 ...

  8. 数据分析与AI(七)傅里叶对登月图片降噪/scipy库对图片进行处理/

    登月图片消噪 scipy.fftpack模块用来计算快速傅里叶变换 速度比传统傅里叶变换更快,是对之前算法的改进 黑白图片是二维数据,注意使用fftpack的二维转变方法 import numpy a ...

  9. Python中的数值运算与逻辑运算

    Python能够实现数值运算和逻辑运算. 1.数值运算 打开Python命令行,输入以下命令: >>>1 + 2 >>>3 >>>3.5 - 4. ...

最新文章

  1. CPU 内部结构解析
  2. 课程 | 中科院教授带你快速入门机器学习
  3. python实现gauss-seidel迭代公式_python实现高斯(Gauss)迭代法的例子
  4. 批阅论文和作业Python程序助手
  5. 【转】centos 6.2 安装mysql-5.5.17
  6. Microsfot.Web.UI.WebControls.TreeView JavaScript控制方法研究(转)
  7. Hibernate 实体映射类的状态值自动转换
  8. php 匹配标记,php – 正则表达式匹配没有标记的链接
  9. html 正则表达式 中文,正则表达式的中文搜索
  10. 前端开发工程师养成记
  11. 功率谱和频谱的区别、联系
  12. C# 线程问题之死锁
  13. 图文|Android 使用Thread 和多线程使用互斥锁
  14. [洛谷2397]yyy loves Maths VI
  15. C#部分---arraylist集合、arraylist集合中的object数据转换成int类string类等;间隔时间的表示方法;...
  16. paip.jquery ajax 请求JSON数据填充SELECT全过程纪录
  17. 车牌号识别 OpenCV
  18. 编程猫编程平台的使用介绍
  19. [翻译]X窗口管理器的原理剖析(一)
  20. H81主板 安装XP 网卡驱动

热门文章

  1. c语言api文档_初学 C 语言没有项目练手?这 20 个小项目拿走不谢
  2. java object 源码_java中Object类 源代码详解
  3. 7添加静态路由 hat red_不同VPC路由器通过静态路由、动态路由(OSPF)实现网络互通实战...
  4. uniapp手写_【转】uni-app框架纯手写微信小程序开发左侧滑动菜单
  5. ubuntu 18.04安装与配置 Redis
  6. python统计缺失值
  7. php 购物车案例教程,php初步实现购物车功能的实例分析
  8. 用虚拟机配置Linux实验环境
  9. CSS——可视化格式模型
  10. BZOJ2654/COGS1764 [2012国家集训队]tree(陈立杰) [生成树,二分]