python数据拟合fit

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

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