Python-曲线拟合（三）

非线性函数拟合

1. 非线性函数不转化成线性函数
2. 直接使用$scipy.optimize$ 下面的$curve$_$fit()$拟合

scipy.optimize.curve_fit(f, xdata, ydata, p0=None, sigma=None)

输入参数：
$f$: 拟合函数
$xdata$,$ydata$: 被拟合数据点
$p0$: 拟合参数初始值
$sigma$: 数据点$ydata$的误差
返回值 :
$popt$: 拟合参数
$pcov$: 拟合参数的协方差
矩阵:
参数标准偏差 :

perr = np.diag(pcov)

实例：
发射源的发射强度具有指数形式$I=I_0{e}^{-\alpha t}$.现有观测数据如下，试确定$I_0$和$\alpha$.

$\begin{array}{clllllll}{\mathbf{t}}_{\mathbf{i}}& 0.2& 0.3& 0.4& 0.5& 0.6& 0.7& 0.8\\ I_i& 3.16& 2.38& 1.75& 1.34& 1& 0.74& 0.56\end{array}$
Python代码实现：

import numpy as np
from scipy.optimize import curve_fit
x = np.arange(0.2, 0.9, 0.1)
y = np.array([3.16,2.38,1.75,1.34,1,0.74,
0.56])
def func(x, a, b): #定义拟合模型
    return a * np.exp(b * x)
pars, pcov = curve_fit(func, x, y)
print(pars) #拟合参数
print(np.sqrt(np.diag(pcov))) #拟合误差

运行结果：

[ 5.63606121 -2.89059333]
[0.03875286 0.01987778]

置信区间

拟合参数的置信区间：
$\begin{array}{l}\widehat{{\beta }_0}={\beta }_0\pm t_{\frac{\alpha }{2}}(n-p){\sigma }_{{\beta }_0}\\ \widehat{{\beta }_1}={\beta }_1\pm t_{\frac{\alpha }{2}}(n-p){\sigma }_{{\beta }_1}\end{array}$
𝑛: 数据点的个数
𝑝: 拟合参数的个数
${\sigma }_{{\beta }_0}$, ${\sigma }_{{\beta }_1}$, 拟合参数的标准偏差
拟合参数的置信区间（Python代码）

import numpy as np
from scipy.stats.distributions import t
from scipy.optimize import curve_fit
x = np.arange(0.2, 0.9, 0.1)
y = np.array([3.16,2.38,1.75,1.34,1,0.74,
0.56])
def func(x, a, b): #定义拟合模型
    return a * np.exp(b * x)
pars, pcov = curve_fit(func, x, y)
print(pars) #拟合参数
print(np.sqrt(np.diag(pcov))) #拟合误差
n = len(y) # number of data points
p = len(pars) # number of parameters
alpha = 0.05 # 95% confidence interval = 100*(1-alpha)
dof = n - p #自由度
# student-t value for the dof and confidence level
tval = t.ppf(1.0-alpha/2., dof) #t_0.975(n-p)
sigma = np.sqrt(np.diag(pcov)) #标准偏差
for i in range(p):
    print('p{0}: {1:5.3f} [{2:5.3f} {3:5.3f}]'\
          .format(i,pars[i], #参数
                  pars[i]-sigma[i]*tval, #置信下限
                  pars[i]+sigma[i]*tval)) #置信上限

运行结果：

[ 5.63606121 -2.89059333]
[0.03875286 0.01987778]
p0: 5.636 [5.536 5.736]
p1: -2.891 [-2.942 -2.839]