线性拟合polyfit_Python曲线拟合详解
转⼀个超级详细的Python曲线拟合详解⽂章(怕以后不到了),本栏⽬初学者不⽤细看,当⼿册查就好了。原⽂在这⾥:04.04 curve fitting,侵删。
导⼊基础包:
In [1]:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
多项式拟合
导⼊线多项式拟合⼯具:
In [2]:
from numpy import polyfit, poly1d
产⽣数据:
In [3]:
x = np.linspace(-5, 5, 100)
y = 4 * x + 1.5
noise_y = y + np.random.randn(y.shape[-1]) * 2.5
画出数据:
In [4]:
%matplotlib inline
p = plt.plot(x, noise_y, 'rx')
p = plt.plot(x, y, 'b:')
进⾏线性拟合,polyfit 是多项式拟合函数,线性拟合即⼀阶多项式:
In [5]:
coeff = polyfit(x, noise_y, 1)
print coeff
[ 3.93921315 1.59379469]
⼀阶多项式 y=a1x+a0 拟合,返回两个系数 [a1,a0]。
画出拟合曲线:
In [6]:
p = plt.plot(x, noise_y, 'rx')
p = plt.plot(x, coeff[0] * x + coeff[1], 'k-')
p = plt.plot(x, y, 'b--')
还可以⽤ poly1d ⽣成⼀个以传⼊的 coeff 为参数的多项式函数:In [7]:
f = poly1d(coeff)
p = plt.plot(x, noise_y, 'rx')
p = plt.plot(x, f(x))
In [8]:
f
Out[8]:
poly1d([ 3.93921315, 1.59379469])
显⽰ f:
In [9]:
print f
3.939 x + 1.594
还可以对它进⾏数学操作⽣成新的多项式:
In [10]:
print f + 2 * f ** 2
2
31.03 x + 29.05 x + 6.674
多项式拟合正弦函数
正弦函数:
In [11]:
x = np.linspace(-np.pi,np.pi,100)
y = np.sin(x)
⽤⼀阶到九阶多项式拟合,类似泰勒展开:
In [12]:
y1 = poly1d(polyfit(x,y,1))
y3 = poly1d(polyfit(x,y,3))
y5 = poly1d(polyfit(x,y,5))
y7 = poly1d(polyfit(x,y,7))
y9 = poly1d(polyfit(x,y,9))
In [13]:
linspace函数pythonx = np.linspace(-3 * np.pi,3 * np.pi,100)
p = plt.plot(x, np.sin(x), 'k')
p = plt.plot(x, y1(x))
p = plt.plot(x, y3(x))
p = plt.plot(x, y5(x))
p = plt.plot(x, y7(x))
p = plt.plot(x, y9(x))
a = plt.axis([-3 * np.pi, 3 * np.pi, -1.25, 1.25])
⿊⾊为原始的图形,可以看到,随着多项式拟合的阶数的增加,曲线与拟合数据的吻合程度在逐渐增⼤。最⼩⼆乘拟合
导⼊相关的模块:
In [14]:
from scipy.linalg import lstsq
from scipy.stats import linregress
In [15]:
x = np.linspace(0,5,100)
y = 0.5 * x + np.random.randn(x.shape[-1]) * 0.35
plt.plot(x,y,'x')
Out[15]:
[<matplotlib.lines.Line2D at 0xbc98518>]
⼀般来书,当我们使⽤⼀个 N-1 阶的多项式拟合这 M 个点时,有这样的关系存在:即
Scipy.linalg.lstsq 最⼩⼆乘解
要得到 C ,可以使⽤ scipy.linalg.lstsq 求最⼩⼆乘解。
这⾥,我们使⽤ 1 阶多项式即 N = 2,先将 x 扩展成 X:
In [16]:
X = np.hstack((x[:,np.newaxis], np.ones((x.shape[-1],1))))
X[1:5]
Out[16]:
array([[ 0.05050505, 1. ],
[ 0.1010101 , 1. ],
[ 0.15151515, 1. ],
[ 0.2020202 , 1. ]])
求解:
In [17]:
C, resid, rank, s = lstsq(X, y)
C, resid, rank, s
Out[17]:
(array([ 0.50432002, 0.0415695 ]),
12.182942535066523,
2,
array([ 30.23732043, 4.82146667]))
画图:
In [18]:
p = plt.plot(x, y, 'rx')
p = plt.plot(x, C[0] * x + C[1], 'k--')
print "sum squared residual = {:.3f}".format(resid)
print "rank of the X matrix = {}".format(rank)
print "singular values of X = {}".format(s)
sum squared residual = 12.183
rank of the X matrix = 2
singular values of X = [ 30.23732043 4.82146667]
Scipy.stats.linregress 线性回归
对于上⾯的问题,还可以使⽤线性回归进⾏求解:
In [19]:
slope, intercept, r_value, p_value, stderr = linregress(x, y)
slope, intercept
Out[19]:
(0.50432001884393252, 0.041569499438028901)
In [20]:
p = plt.plot(x, y, 'rx')
p = plt.plot(x, slope * x + intercept, 'k--')
print "R-value = {:.3f}".format(r_value)
print "p-value (probability there is no correlation) = {:.3e}".format(p_value) print "Root mean squared error of the fit = {:.3f}".format(np.sqrt(stderr))
R-value = 0.903
p-value (probability there is no correlation) = 8.225e-38
Root mean squared error of the fit = 0.156
版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。
发表评论