利用Python进行曲线拟合

曲线拟合

参考博客：http://lijin-thu.github.io/04.%20scipy/04.04%20curve%20fitting.html

最小二乘多项式拟合

    寻找一个曲线满足$y=f(x)$，线性最小二乘法拟合是解决曲线拟合最常用的方法，基本思路是令

$$f(x)=a_1{r}_1(x)+a_2{r}_2(x)+\cdots +a_m{r}_m(x)$$

式中：$r_k(x)$为实现选定的一组线性无关函数；$a_k$为待定系数($k=1,2,3,\cdots ,m;m<n$)

    拟合的准则是让$y_i(i=1,2,3,\cdots ,m)$与$f(x_i)$的距离的平方和最小，我们将这个条件成为最小二乘准则。这个就不详细表述相关的原理，值得注意的是用最小二乘多项式法进行拟合时，我们要注意的就是$r(x)$的取值常用的曲线有一下几类：

• $y=a_{0}x+a_1$
• 多项式 $y=a_1{x}^m+a_2{x}^{m-1}+\cdots +a_{m}x+a_{m+1}$
• 一支双曲线 $y=\frac{a_1}{x}+a_2$
• 指数函数 $y=a_1{e}^{a_{2}x}$

Python中也有相应的拟合方法，首先介绍第一个线性多项式拟合工具polyfit，首先导入需要用到的包：

from numpy import polyfit,poly1d
import numpy as np
import matplotlib.pyplot as plt

介绍一下polyfit这个方法，现在的需求只需要前三个参数，日常传送门

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

第一个X参数和第二个Y参数表示的都是数据点的(x,y)坐标值，第三个deg表示的是要拟合的函数的阶数，如填写1则为线性，2及其以上都为多项式拟合。
例子：

#导入所需要用到的包
from numpy import polyfit,poly1d
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0,10,100)#生成100个在区间[0,10]等间距的点
y = x**2+x#计算y值
y_noise =y + np.random.randn(y.shape[0])*5 #产生噪音y,randn是生成一个形状如y的随机数矩阵，其中矩阵中元素的值在[-1,1]之间，如果生成非负随机数则使用rand(numpt,shape)方法
plt.plot(x,y,'--',x,y_noise,'r+')
plt.legend(['Origin','Noise'])
plt.show()

效果图：

下面我们要对噪声点进行各种次数的多项式拟合：

x = np.linspace(0,10,100)
y = x**2+x
y_noise =y + np.random.randn(y.shape[0])*5
fit_result_1 = np.polyfit(x,y_noise,1)#进行一次拟合，这个只能返回一个拟合后的系数矩阵，如：[11.256762634722696, -18.173869106185016]，并不能生成图像函数
fit_ploy_1 = np.poly1d(fit_result_1)#将根据系数矩阵，转化成多项式函数
fit_result_2 = np.polyfit(x,y_noise,2)
fit_ploy_2 = np.poly1d(fit_result_2)
fit_result_5 = np.polyfit(x,y_noise,5)
fit_ploy_5 = np.poly1d(fit_result_5)
plt.plot(x,y_noise,'r+',x,fit_ploy_1(x),'--',x,fit_ploy_2(x))
plt.plot(x,fit_ploy_5(x))
plt.legend(['Noise','d-1','d-2','d-5'],loc='upper left')
plt.show()

效果图：

非线性最小二乘

上面的函数比较简单好用，但是这个也有局限性就是只能通过多项式函数或者线性函数来进行拟合，那么有些指数函数，或者一些定制化函数怎么进行拟合呢？scipy.optimize模块提供了非线性最小二乘估计方法curve_fit。
介绍一下这个方法：

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)

解释一下我觉得重要的参数的意思：
前三个非常重要，f为要拟合的函数，要求必须是独立出来的，内部第一个参数为x的取值，后面全是参数，如下所示：

ydata = f(xdata, *params) + eps.

xdata就是观测点选取得值；ydata就是观测出的函数值；p0为初次猜想的参数得大致取值，不选默认都是1；bounds是边界值，指的是参数的边界，是一个元组(最小值，最大值)，如果不给默认为正负无穷；
method指的是采用什么方法进行拟合，有三种可选择的方法这三种算法我也不是特别了解，具体细节网上搜吧，我这里只说个大概的情况信息传送门：

• ‘lm’：Levenberg-Marquardt algorithm as implemented in MINPACK. Doesn’t handle bounds and sparse Jacobians. Usually the most efficient method for small unconstrained problems.
• ‘trf’ : Trust Region Reflective algorithm, particularly suitable for large sparse problems with bounds. Generally robust method.
• ‘dogbox’ : dogleg algorithm with rectangular trust regions, typical use case is small problems with bounds. Not recommended for problems with rank-deficient Jacobian.

默认得方法就是‘lm’，实在翻译水平有限，怕理解错误。
例子：

def function(x, a , b, f, phi):
    return a * np.exp(-b * np.sin(f * x + phi))
real_param = [3,2,1.25,np.pi/4]
x = np.linspace(0, 2 * np.pi, 100)
y = function(x,*real_param)
y_noise = y+np.random.randn(y.shape[0])
p_est_d,err_est_d = curve_fit(function,x,y_noise,method= 'dogbox')
p_est_l,err_est_l = curve_fit(function,x,y_noise,method='lm')
p_est_t,err_est_t = curve_fit(function,x,y_noise,method = 'trf')
plt.plot(x,y,'--',x,y_noise,'r+')
plt.plot(x,function(x,*p_est_d),x,function(x,*p_est_l),x,function(x,*p_est_t),'o--')
plt.legend(['real','noise','dogbox','lm','trf'],loc='best')
print("dogbox:",err_est_d)
print("lm:",err_est_l)
print("trf:",err_est_t)
print ("   a\t  b\t f\tphi")
print(err_est_d.diagonal())
plt.show()

图像显示：

解释一下：这个并不是只有最后一个图像画上去了，而是三个拟合曲线基本完全一致，我放大到y值差别小于0.0001，才能看见微小的差别。下面给出了四个参数两两之间得协方差矩阵，仔细观察发现也没多大差别。

dogbox: 
[[ 0.17292217 -0.05669224  0.02034361 -0.06410764]
 [-0.05669224  0.0188186  -0.00653609  0.02059677]
 [ 0.02034361 -0.00653609  0.00281707 -0.00887735]
 [-0.06410764  0.02059677 -0.00887735  0.02814787]]
lm: 
[[ 0.17286838 -0.05667102  0.02033665 -0.06408569]
 [-0.05667102  0.01881047 -0.00653339  0.02058824]
 [ 0.02033665 -0.00653339  0.00281619 -0.00887458]
 [-0.06408569  0.02058824 -0.00887458  0.02813912]]
trf: 
[[ 0.17292369 -0.05669284  0.0203438  -0.06410826]
 [-0.05669284  0.01881883 -0.00653617  0.02059701]
 [ 0.0203438  -0.00653617  0.00281709 -0.00887743]
 [-0.06410826  0.02059701 -0.00887743  0.02814811]]
 #对角线表示得是各个参数方差。
   a      b      f      phi
[0.17292217 0.0188186  0.00281707 0.02814787]

其实对于线性拟合问题我们还有其他方法：

1. 比如使用Scipy.stats.linregress来做线性回归
2. 使用Scipy.linalg.lstsq来求出最小二乘解
   这些都是可以解决问题的。