python数据拟合

1.多项式拟合

1.1 多项式拟合描述

输入：多项式次数$n$；数据集合 { ( x i , y i ) } \{(x_i,y_i)\} {(xi​,yi​)}
输出：$f_n(x)$，使得 $\sum (f(x_i)-y_i{)}^2$ 最小

解决方案：polyfit ( x , y , deg , rcond = None , full = False , w = None )

1.2 多项式拟合实现

导入库：numoy、matplotlib

#导入库
import numpy as np
import matplotlib.pyplot as plt

生成测试数据：

#定义测试多项式函数
def func1(x):
    return 3*x*x*x-2*x*x+4

#生成测试数组
x1 = np.array([1,2,3,4,5,6,7,8])
y1 = func1(x1)

#添加噪声
n1 = np.random.normal(0.0,1.0,8)*0.1
y1 = y1*(1+n1)

数据拟合：

f1 = np.polyfit(x1, y1, 3)
p1 = np.poly1d(f1)
print('p1 is :\n',p1)

p1 is :
3 2
-1.897 x + 57.77 x - 204.8 x + 177.6

得到拟合函数：
$f_1=-1.897x^3+57.77x^2-204.8x+177.6$

绘制拟合曲线：

xx1 = np.arange(1,9,0.2)
yvals1 = p1(xx1) #拟合y值

#绘图
plot1 = plt.plot(x1, y1, 'o',label='original values')
plot2 = plt.plot(xx1, yvals1, 'r-',label='polyfit values')
plt.xlabel('x1')
plt.ylabel('y1')
plt.legend(loc=4) #指定legend的位置右下角
plt.title('polyfitting')
plt.show()

2.自定义函数拟合

2.1 自定义函数拟合描述

输入：自定义函数$f(x)$及参数${\sigma }_i$；数据集合 { ( x i , y i ) } \{(x_i,y_i)\} {(xi​,yi​)}
输出：使得 $\sum (f(x_i)-y_i{)}^2$ 最小参数取值

解决方案：scipy.optimize.curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, check_finite=True, bounds=(- inf, inf), method=None, jac=None, **kwargs)

2.1 自定义函数拟合的实现

导入库：numpy、matplotlib、scipy

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

设置拟合数据：

x2 = np.array([1,2,3,4,5,6,7,8,9,10])
y2 = np.array([0.16,0.63,1.60,3.00,8.00,33.0,73.0,125.0,211.0,310.0])
y2=y2*100

自定义拟合函数：

def func2(x, p, q,m):
    return m*(1-np.exp(-x*(p+q)))/(1+q/p*np.exp(-(p+q)*x))

非线性最小二乘法拟合：

popt, pcov = curve_fit(func2, x2, y2)

拟合y值

p = popt[0] 
q = popt[1]
m = popt[2]

yvals2 = func2(x2,p,q,m) #拟合y值
print('popt:', popt)
print('系数p:', p)
print('系数q:', q)
print('系数m:', m)
print('系数pcov:', pcov)
print('系数yvals2:', yvals2)

popt: [4.93963593e-04 7.86873973e-01 4.96871803e+04]
系数p: 0.0004939635925128215
系数q: 0.7868739729224568
系数m: 49687.18030400891
系数pcov: [[ 1.16165013e-08 -4.88182871e-06 3.60534694e-01]
[-4.88182871e-06 2.28100736e-03 -1.95238501e+02]
[ 3.60534694e-01 -1.95238501e+02 1.97747550e+07]]
系数yvals2: [ 37.30343602 119.0851606 297.86586198 686.25958059
1518.69723981 3252.26832861 6655.68174484 12625.30614774
21284.37699301 30920.18607139]

绘制拟合图像

xx2 = np.arange(1,21)
y2test = func2(xx2,p,q,m)

#绘图
plot1 = plt.plot(x2, y2, 's',label='original values')
plot2 = plt.plot(xx2, y2test, 'r',label='polyfit values')
plt.xlabel('x2')
plt.ylabel('y2')
plt.legend(loc=4) #指定legend的位置右下角
plt.title('curve_fit')
plt.show()

使用Python进行数据拟合

多项式拟合

任何一个函数都可以拆分成近似于这个函数的多项式表达。

多项式拟合需要用到的函数是numpy库当中的np.polyfit，它的使用方法为：

np.polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False)

使用最小二乘法原理，根据已知的x与y对应值，拟合一个下面形式的多项式。
$$P(x)=P{}_0{x}^{deg}+P{}_1{x}^{deg-1}\cdot \cdot \cdot +P_{deg}$$
返回一系列的系数$P$.

参数说明：
一般情况下，我们只需要用到前三个参数。

返回值：

我们还可以使用polyval()来计算我们需要预测多项式的值.

import matplotlib.pyplot as plt
import numpy as np

if __name__ == "__main__":
  x = np.arange(1, 31, 1)
  y = np.array([20, 23, 26, 29, 32, 35, 38, 45, 53, 62, 73, 86, 101, 118, 138, 161, 188, 220, 257, 300, 350, 409, 478, 558, 651, 760, 887, 1035, 1208, 1410])

  z1 = np.polyfit(x, y, 3)              # 曲线拟合，返回值为多项式的各项系数
  p1 = np.poly1d(z1)                    # 返回值为多项式的表达式，也就是函数式子
  print(p1)
  y_pred = p1(x)                        # 根据函数的多项式表达式，求解 y
  #print(np.polyval(p1, 29))             #根据多项式求解特定 x 对应的 y 值
  #print(np.polyval(z1, 29))             #根据多项式求解特定 x 对应的 y 值
  plot1 = plt.plot(x, y, '*', label='original values')
  plot2 = plt.plot(x, y_pred, 'r', label='fit values')
  plt.title('')
  plt.xlabel('')
  plt.ylabel('')
  plt.legend(loc=3, borderaxespad=0., bbox_to_anchor=(0, 0))
  plt.show()

输出结果:

        3         2
0.1215 x - 3.045 x + 28.62 x - 34.47

则其拟合的函数为$$y=0.1215x{}^3-3.045x{}^2+28.62x-34.47$$
将拟合的曲线与原数据各点进行对比发现拟合效果良好。

非多项式拟合

如果需要进行多项式拟合，前提是必须大体上知道散点的大致曲线形式，也就是大致的函数的形式。
比如，例子中的散点看起来像是指数的函数分布，因此可以给出假设的函数
$$y=ba{}^x+c$$
所以，只要给出具体的函数形式(可以是任意的，只要能写的出来皆可)，用最小二乘的方式去逼近和拟合，即求出函数的各项系数。

此时用到的是scipy.optimize包下的curve_fit函数了：

代码示例：

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

def func(x, a, b, c):
  return b * np.power(a, x) + c

if __name__ == "__main__":
  x = np.arange(1, 31, 1)
  y = np.array([20, 23, 26, 29, 32, 35, 38, 45, 53, 62, 73, 86, 101, 118, 138, 161, 188, 220, 257, 300, 350, 409, 478, 558, 651, 760, 887, 1035, 1208, 1410])

  popt, pcov = curve_fit(func, x, y)                # 曲线拟合，popt为函数的参数list
  y_pred = [func(i, popt[0], popt[1], popt[2]) for i in x]    # 直接用函数和函数参数list来进行y值的计算
  print(popt)

  plot1 = plt.plot(x, y, '*', label='original values')
  plot2 = plt.plot(x, y_pred, 'r', label='fit values')
  plt.title('')
  plt.xlabel('')
  plt.ylabel('')
  plt.legend(loc=3, borderaxespad=0., bbox_to_anchor=(0, 0))
  plt.show()

输出结果为：

[ 1.16791847 13.39168878  1.24633841]

那么就有
$$\{\begin{array}{rlr}a& =& 1.16791847\\ b& =& 13.39168878\\ c& =& 1.24633841\end{array}$$
于是我们得到了拟合的函数为

$$y=13.39168878\ast 1.16791847{}^x+1.24633841$$
将拟合的曲线与原数据各点进行对比发现拟合效果良好。

第0部分：多项式拟合数学基础

参考文献
多项式拟合采用的是最小二乘拟合

这里最重要的就是平方误差条件和公式(4)。
公式4表明，
1） 我们在计算系数a的时候可以直接通过矩阵来计算。
2） 对两条曲线和的拟合等于单独对两条曲线拟合的和。
比如，现在有两组数据$(x_i,y_i)和(x_i,z_i)$,多项式对单组数据拟合和对$(x_i,y_i+z_i)$的拟合是可以直接相加的，我们可以由（4）式子有

上面的图片告诉我们一个很重要的信息，即有自变量$x_i$组成的矩阵是一样的，那么我们在对$(x_i,y_i+z_i)$的拟合系数等于单独对$(x_i,y_i)$的拟合系数与$(x_i,z_i)$的拟合系数的和。这表明当我们对两条曲线的和做拟合时得到的拟合曲线是$p_n(x)={\sum }_{k=0}^n(a_k+b_k)x^k$,那么两条曲线与拟合曲线相减，剩余为$w_i=(y+z{)}_i-p_n^i(x)=(y+z{)}_i-{\sum }_{k=0}^n(a_k+b_k)x^k$,如果要计算其中一条曲线被减了多少，或者要知道剩余中每条曲线分别还含有多少比例，那么可以放心大胆的再次使用这个$n$阶多项式单独对每条曲线进行fit再相减，即$y$的剩余为$y_i-{\sum }_{k=0}^n{a}_k{x}^k$,$z$的剩余为$z_i-{\sum }_{k=0}^n{b}_k{x}^k$.这里要注意必须使用相同阶数，比如对$(y+z)$使用10阶多项式拟合，那么对单个成分也必须是10阶。

举例

第一部分：多项式拟合

• 1）一次式

import numpy as np
import matplotlib.pyplot as plt

x,y = [1,1.25,1.33,1.60,2.0,2.4,2.5,3.0,3.2],[1000,1100,1300,1500,1861.48\
      ,2233.6,2326.64,2791.94,2978.09]
x,y = np.array(x),np.array(y) #变成数组形式，方便数学计算
z1 = np.polyfit(x,y,1)   #使用一次式拟合,使用N次式则将1改为N即可
p1 = np.poly1d(z1)    #一维
print(p1)         #输出p1的式子
y_fit = p1(x)      #fit之后的结果

    
plt.scatter(x,y_fit,label='fit')    #画fit后的图
plt.scatter(x,y,label='data')    #画原始图

plt.legend()
plt.show()

输出的一次函数形式：

得到的图像：

第二部分

最小二乘法拟合（参考python科学计算）

使用幂律谱

幂律谱有如下形式：
$$y=a\times x^b$$
如果用矢量p表示函数中需要确定的参数，则目标是找到一组p使得函数S的值最小：
$$S(\mathbf{p})=\sum _{i=1}^m[y_i-f(x_i,\mathbf{p}){]}^2$$

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq

def residuals(p):
    a,b = p    # p表示参数，所以这里有两个参数a,b（不可少）
    return Ad-(a*fsky**b)   #输出真实值与拟合值的差，这里使用的是幂律谱y_fit = a*x**b

fsky = np.array([0.24,0.33,0.42,0.52,0.62,0.71])
Ad = np.array([0.48*34.3,0.45*47.3,0.50*74.7,0.53*120.1,0.53*190.7,0.53*315.4])

popt = leastsq(residuals,[1,0])  #这里给出residual得到的数组的平方和，即上面公式的S（p）
#[1,0]分别是a,b的初始值
a,b = popt[0]   # 两个参数
x = np.linspace(0,1,10)
Ad_fit = a*x**b

plt.scatter(fsky,Ad)
plt.plot(x,Ad_fit)
plt.show()

更新（2021.7.26）
更新对于log坐标下的幂律谱拟合
比如有个$\ell$和$C_{\ell }$，它们在log坐标下的关系如下：

现在要用幂律谱拟合这条曲线，由于是在ipython环境，所以就直接这么写了

使用e指数

e指数有如下形式
$$y=a\times e^{b\times x}$$

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq

def residuals(p):
    a,b = p
    return Ad-(a*np.exp(b*fsky))

fsky = np.array([0.24,0.33,0.42,0.52,0.62,0.71])
Ad = np.array([0.48*34.3,0.45*47.3,0.50*74.7,0.53*120.1,0.53*190.7,0.53*315.4])

popt = leastsq(residuals,[1,0])
a,b = popt[0]
x = np.linspace(0,1,10)
Ad_fit = a*np.exp(b*x)

plt.scatter(fsky,Ad)
plt.plot(x,Ad_fit)
plt.show()

三种方法总结

将上面提到的三种方法统一写在下面的这个代码中，对输入fsky, Ad做拟合

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq

def polynomial_fit(x,y,N):
    return np.polyfit(x,y,N)

def exponential(p):
    a,b = p
    return Ad-(a*np.exp(b*fsky))

def power_law(p):
    a,b = p
    return Ad-(a*fsky**b)

fsky = np.array([0.24,0.33,0.42,0.52,0.62,0.71])
Ad = np.array([0.48*34.3,0.45*47.3,0.50*74.7,0.53*120.1,0.53*190.7,0.53*315.4])
# x range
x = np.linspace(0,1,100)

# polynomial_fit
popt_poly = polynomial_fit(fsky,Ad,3)
p1 = np.poly1d(popt_poly)
Ad_poly = p1(x)
print(p1)

#exponential form
popt_exp = leastsq(exponential,[1,0])
a_exp,b_exp = popt_exp[0]
#power law
popt_power_law = leastsq(power_law,[1,0])
a_pow,b_pow = popt_power_law[0]

Ad_fit_exp = a_exp*np.exp(b_exp*x)
Ad_fit_pow = a_pow*x**b_pow

#plt.figure()
plt.scatter(fsky,Ad)
plt.plot(x,Ad_poly,label='polynomial fit')
plt.plot(x,Ad_fit_exp,label='exponential fit')
plt.plot(x,Ad_fit_pow,label='power_law fit')
plt.xlabel('fsky')
plt.ylabel('Ad')

plt.legend()
plt.show()

例题：拟合中国GDP最近10年的增速，计算多少年可以赶上美国

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq

def residuals(p):
    a,b = p    # p表示参数，所以这里有两个参数a,b（不可少）
    return(Ad-(a*(fsky-2001)**b))   #输出真实值与拟合值的差，这里使用的是幂律谱y_fit = a*x**b

fsky = np.array([2010,2011,2012,2013,2014,2015,2016,2017,2018,2019])  #这里是时间，我就直接使用sky表示了
Ad = np.array([10.9,9.5,7.9,7.8,7.3,6.9,6.7,6.8,6.6,6.3])#这里是中国GDP最近10年增速

popt = leastsq(residuals,[1,1])  #这里给出residual得到的数组的平方和，即S（p）
a,b = popt[0]   # 两个参数
x = np.arange(2010,2030)
Ad_fit = a*(x-2001)**b

def residualss(p):
    c,d = p
    return(Ad-(c*np.exp(d*(fsky-2001))))

popts = leastsq(residualss,[1,0])
c,d = popts[0]
Ad_fit1 = c*np.exp(d*(x-2001))

plt.plot(x,Ad_fit1,'k',label='exponential fit')
plt.scatter(fsky,Ad)
plt.plot(x,Ad_fit,'r',label='power_law fit')
plt.legend()
plt.show()

这种情况需要到2038年

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq

#power law 形式
def residuals(p):
    a,b = p    # p表示参数，所以这里有两个参数a,b（不可少）
    return(Ad-(a*(fsky-2001)**b))   #输出真实值与拟合值的差，这里使用的是幂律谱y_fit = a*x**b

fsky = np.array([2010,2011,2012,2013,2014,2015,2016,2017,2018,2019])
Ad = np.array([10.9,9.5,7.9,7.8,7.3,6.9,6.7,6.8,6.6,6.3])

popt = leastsq(residuals,[1,1])  #这里给出residual得到的数组的平方和，即S（p）
a,b = popt[0]   # 两个参数
x = np.arange(2010,2040)
Ad_fit = a*(x-2001)**b


#e指数形式
def residualss(p):
    c,d = p
    return(Ad-(c*np.exp(d*(fsky-2001))))

popts = leastsq(residualss,[1,0])
c,d = popts[0]
Ad_fit1 = c*np.exp(d*(x-2001))

plt.plot(x,Ad_fit,'r',label='power_law fit')
plt.plot(x,Ad_fit1,'k',label='exponential fit')
plt.scatter(fsky,Ad)

plt.legend()
plt.show()

x1 = np.arange(2018,2045)
w = a*(x1-2001)**b
#使用power law 形式，乘以0.001是从百分比化成小数
#计算需要用到append函数来保存
y1 = [13.6]
w = 13.6
for i in range(2019,2045):
    w = w*(1+0.01*a*(i-2001)**b)
    y1.append(w)
y1 = y1
y2 = 20.5*np.power(1+0.023,x1-2018)
plt.figure()
plt.plot(x1,y1,label='China')
plt.plot(x1,y2,label='United States')
plt.legend()
plt.show()

第三部分：使用窗口平滑化处理(scipy.signal.convolve)

如果信号噪声比较多，或者图片/数组中高频部分多，那么可以考虑使用scipy的卷积来平滑smooth，

#convolution (smoothing)
win = signal.hann(10)
y1 = signal.convolve(y,win,mode='same')/sum(win)

下图蓝色是平滑前，黄色是平滑后，可以降低噪声水平