python数据拟合fit

文章目录

第0部分：多项式拟合数学基础
- 举例
第一部分：多项式拟合
第二部分
- 最小二乘法拟合（参考python科学计算）
- - 使用幂律谱
  - 使用e指数
  - 三种方法总结
第三部分：使用窗口平滑化处理(scipy.signal.convolve)

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

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

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

上面的图片告诉我们一个很重要的信息，即有自变量xix_ixi组成的矩阵是一样的，那么我们在对(xi,yi+zi)(x_i,y_i+z_i)(xi,yi+zi)的拟合系数等于单独对(xi,yi)(x_i,y_i)(xi,yi)的拟合系数与(xi,zi)(x_i,z_i)(xi,zi)的拟合系数的和。这表明当我们对两条曲线的和做拟合时得到的拟合曲线是pn(x)=∑k=0n(ak+bk)xkp_n(x)=\sum_{k=0}^n(a_k+b_k)x^kpn(x)=∑k=0n(ak+bk)xk,那么两条曲线与拟合曲线相减，剩余为wi=(y+z)i−pni(x)=(y+z)i−∑k=0n(ak+bk)xkw_i=(y+z)_i-p_n^i(x)=(y+z)_i-\sum_{k=0}^n(a_k+b_k)x^kwi=(y+z)i−pni(x)=(y+z)i−∑k=0n(ak+bk)xk,如果要计算其中一条曲线被减了多少，或者要知道剩余中每条曲线分别还含有多少比例，那么可以放心大胆的再次使用这个nnn阶多项式单独对每条曲线进行fit再相减，即yyy的剩余为yi−∑k=0nakxky_i-\sum_{k=0}^na_kx^kyi−∑k=0nakxk,zzz的剩余为zi−∑k=0nbkxkz_i-\sum_{k=0}^nb_kx^kzi−∑k=0nbkxk.这里要注意必须使用相同阶数，比如对(y+z)(y+z)(y+z)使用10阶多项式拟合，那么对单个成分也必须是10阶。

举例

第一部分：多项式拟合

1）一次式

import numpy as np
import matplotlib.pyplot as pltx,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×xby=a\times x^by=a×xb
如果用矢量p表示函数中需要确定的参数，则目标是找到一组p使得函数S的值最小：
S(p)=∑i=1m[yi−f(xi,p)]2S(\boldsymbol p)=\sum_{i=1}^{m}[y_i-f(x_i,\boldsymbol p)]^2S(p)=i=1∑m[yi−f(xi,p)]2

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsqdef residuals(p):a,b = p    # p表示参数，所以这里有两个参数a,b（不可少）return Ad-(a*fsky**b)   #输出真实值与拟合值的差，这里使用的是幂律谱y_fit = a*x**bfsky = 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**bplt.scatter(fsky,Ad)
plt.plot(x,Ad_fit)
plt.show()

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

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

使用e指数

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

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsqdef residuals(p):a,b = preturn 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 leastsqdef polynomial_fit(x,y,N):return np.polyfit(x,y,N)def exponential(p):a,b = preturn Ad-(a*np.exp(b*fsky))def power_law(p):a,b = preturn 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 leastsqdef residuals(p):a,b = p    # p表示参数，所以这里有两个参数a,b（不可少）return(Ad-(a*(fsky-2001)**b))   #输出真实值与拟合值的差，这里使用的是幂律谱y_fit = a*x**bfsky = 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)**bdef residualss(p):c,d = preturn(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**bfsky = 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 = preturn(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)

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