文章目录
第0部分:多项式拟合数学基础
参考文献
多项式拟合采用的是最小二乘拟合
这里最重要的就是平方误差条件和公式(4)。
公式4表明,
1) 我们在计算系数a的时候可以直接通过矩阵来计算。
2) 对两条曲线和的拟合等于单独对两条曲线拟合的和。
比如,现在有两组数据
(
x
i
,
y
i
)
和
(
x
i
,
z
i
)
(x_i,y_i)和(x_i,z_i)
(xi,yi)和(xi,zi),多项式对单组数据拟合和对
(
x
i
,
y
i
+
z
i
)
(x_i,y_i+z_i)
(xi,yi+zi)的拟合是可以直接相加的,我们可以由(4)式子有
上面的图片告诉我们一个很重要的信息,即有自变量
x
i
x_i
xi组成的矩阵是一样的,那么我们在对
(
x
i
,
y
i
+
z
i
)
(x_i,y_i+z_i)
(xi,yi+zi)的拟合系数等于单独对
(
x
i
,
y
i
)
(x_i,y_i)
(xi,yi)的拟合系数与
(
x
i
,
z
i
)
(x_i,z_i)
(xi,zi)的拟合系数的和。这表明当我们对两条曲线的和做拟合时得到的拟合曲线是
p
n
(
x
)
=
∑
k
=
0
n
(
a
k
+
b
k
)
x
k
p_n(x)=\sum_{k=0}^n(a_k+b_k)x^k
pn(x)=∑k=0n(ak+bk)xk,那么两条曲线与拟合曲线相减,剩余为
w
i
=
(
y
+
z
)
i
−
p
n
i
(
x
)
=
(
y
+
z
)
i
−
∑
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
wi=(y+z)i−pni(x)=(y+z)i−∑k=0n(ak+bk)xk,如果要计算其中一条曲线被减了多少,或者要知道剩余中每条曲线分别还含有多少比例,那么可以放心大胆的再次使用这个
n
n
n阶多项式单独对每条曲线进行fit再相减,即
y
y
y的剩余为
y
i
−
∑
k
=
0
n
a
k
x
k
y_i-\sum_{k=0}^na_kx^k
yi−∑k=0nakxk,
z
z
z的剩余为
z
i
−
∑
k
=
0
n
b
k
x
k
z_i-\sum_{k=0}^nb_kx^k
zi−∑k=0nbkxk.这里要注意必须使用相同阶数,比如对
(
y
+
z
)
(y+z)
(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
×
x
b
y=a\times x^b
y=a×xb
如果用矢量p表示函数中需要确定的参数,则目标是找到一组p使得函数S的值最小:
S
(
p
)
=
∑
i
=
1
m
[
y
i
−
f
(
x
i
,
p
)
]
2
S(\boldsymbol p)=\sum_{i=1}^{m}[y_i-f(x_i,\boldsymbol p)]^2
S(p)=i=1∑m[yi−f(xi,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
ℓ
C_\ell
Cℓ,它们在log坐标下的关系如下:
现在要用幂律谱拟合这条曲线,由于是在ipython环境,所以就直接这么写了
使用e指数
e指数有如下形式
y
=
a
×
e
b
×
x
y=a\times e^{b\times x}
y=a×eb×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)
下图蓝色是平滑前,黄色是平滑后,可以降低噪声水平