曲线拟合

曲线拟合

理论在WIKIPEDIA里介绍的很详细

Curve fitting - Wikipedia

数据拟合的最小二乘法

数据拟合的最小二乘问题是：根据给定的数据组$(x_i,y_i),i=1,2,\dots ,n,$选取近似函数形式，即给定函数类$\phi (x)\in H$,使得
$$\sum _{i=1}^n{\delta }_i^2=\sum _{i=1}^n[y_i-\phi (x_i){]}^2$$
最小，即
$$\sum _{i=1}^n[y_i-\phi (x_i){]}^2=\underset{\phi \in H}{\min }\sum _{i=1}^n[y_i-\varphi (x_i){]}^2$$
函数$\phi (x)$称为这组数据的最小二乘函数，通常取$H$为一些比较简单函数的集合，如低次多项式，指数函数等。通常是先把已知数据先绘图出来，观察其数据分布特征，从直观上选择用什么拟合函数合适，当然，很多时候也会有分成几段使用不同拟合函数来交叉拟合的情况。

我是从二次多项式（$\phi (x)=a_0+a_{1}x+a_2{x}^2$）拟合举例说明这个曲线拟合到底是个什么流程。个人感觉理解效果会好点。

观察数据特征

import numpy as np
import matplotlib.pyplot as plt 
x = np.array([1,3,4,5,6,7,8,9,10]).reshape(-1,1)
y = np.array([10,5,4,2,1,1,2,3,4]).reshape(-1,1)
plt.plot(x,y)
plt.plot(x,y,'ro', label = 'Original data')
plt.legend()

选取拟合函数

通过观察上图结果，我们其实肉眼就可以看出来，有点像我们学过的二次函数。所以这里我们选择二次多项式来拟合。
$$\phi (x)=a_0+a_{1}x+a_2{x}^2$$
其中，$\phi (x_1)=1,\phi (x_2)=x,\phi (x_3)=x^2$

def fun(x,a0,a1,a2):
    y = a0 + a1 * x + a2 * x**2
    return y

拟合

正则方程组为：

$${\Phi }^T\Phi a={\Phi }^{T}y$$
其中，
$${\Phi }^T\Phi =\left(\begin{array}{cccc}({\phi }_0,{\phi }_0)& ({\phi }_0,{\phi }_1)& \dots & ({\phi }_0,{\phi }_n)\\ ({\phi }_1,{\phi }_0)& ({\phi }_1,{\phi }_1)& \dots & ({\phi }_1,{\phi }_n)\\ ⋮& \ddots & ⋮& \\ ({\phi }_n,{\phi }_0)& ({\phi }_n,{\phi }_1)& \dots & ({\phi }_n,{\phi }_n)\end{array}\right)$$
$$a=[a_1,a_2,\dots ,a_n{]}^T$$

根据最小二乘估计
$$a=({\Phi }^T\Phi {)}^{-1}{\Phi }^{T}y$$

def get_params(x,y):
    phi_matrix = np.hstack((np.ones(x.shape[0]).reshape(-1,1),x,x ** 2))
    A = phi_matrix.T @ phi_matrix
    b = phi_matrix.T @ y
    params = np.linalg.pinv(A) @ b
    return params

params = get_params(x,y)
print(params)

[[13.45966387]
 [-3.6053094 ]
 [ 0.26757066]]

根据上述求解参数，拟合函数为
$$\varphi (x)=13.45966387-3.6053094x+0.26757066x^2$$

带入拟合函数，求一下根据我们拟合函数得到的y值，记为y_fitted

y_fitted = fun(x,params[0],params[1],params[2])
print(y_fitted)

[[10.12192513]
 [ 5.05187166]
 [ 3.31955691]
 [ 2.1223835 ]
 [ 1.46035141]
 [ 1.33346066]
 [ 1.74171123]
 [ 2.68510313]
 [ 4.16363636]]

其实，直接直观看值，我们就可以发现，拟合效果已经很好了，误差都比较小。但为了更清晰的看到结果，我们通过MAPE(%)来看一下拟合效果。

def mape(y,y_fitted):
    e = (sum(abs((y - y_fitted)/y)) / len(y)) * 100
    return e
print(mape(y,y_fitted)[0],"%")

14.69667260842274 %

画图对比一下原始数据图像与拟合效果图

plt.plot(x,y)
plt.plot(x,y,'ro', label = 'Original data')
plt.plot(x,y_fitted,'-b', label ='Fitted')
plt.legend()
plt.show()

可以看到，拟合效果还是比较好的,当然，现在的工具是非常多的。如果不是要了解原理都话，直接调用命令也是非常不错的。下面给出一些拟合的直接调用方法。也根据标准的对照参考一下，刚刚我用理论求解得到的结果正确性。

numpy ployfit包，多项式拟合

import numpy as np
import matplotlib.pyplot as plt 
x = np.array([1,3,4,5,6,7,8,9,10])
y = np.array([10,5,4,2,1,1,2,3,4])
plt.plot(x,y)

params = np.polyfit(x, y, deg = 2)
print(params)

def fun(x):
    y = params[0]*x**2 + params[1]*x + params[2]
    return y
y_fitted = fun(x)

print(mape(y,y_fitted),"%")
plt.plot(x,y,'ro', label = 'Original data')
plt.plot(x,y_fitted,'-b', label ='Fitted data')
plt.legend()
plt.show()

[ 0.26757066 -3.6053094  13.45966387]
14.696672608437163 %

可以看到，和我前面写的结果是一模一样的。但是这里有一点细节要注意，这里的拟合函数形式是
$\phi (x)=a_0{x}^2+a_{1}x+a_2$ 和我前面写的$\phi (x)=a_0+a_{1}x+a_2{x}^2$要注意区别哪个参数是哪项的参数。

Scipy的leastsq和curve_fit拟合

直接掉包的事情，直接看代码,这里参考了知乎大佬的，

Python-最小二乘法曲线拟合 - 知乎 (zhihu.com)

我自己懒没有写，只是稍微改了丢丢哈哈哈哈

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
 
def Fun(p,x):                        
    a1,a2,a3 = p
    return a1*x**2+a2*x+a3
def error (p,x,y):                   
    return Fun(p,x)-y 
def fit():
    x = np.array([1,3,4,5,6,7,8,9,10])
    p_value = [1,1,50] # 原始数据的参数,随便设定
    y = np.array([10,5,4,2,1,1,2,3,4])
    plt.plot(x,y)
    p0 = [0.1,-0.01,100] # 拟合的初始参数设置，随便设定
    para =leastsq(error, p0, args=(x,y)) 
    y_fitted = Fun (para[0],x) 
 
    plt.plot(x,y,'r', label = 'Original data')
    plt.plot(x,y_fitted,'-b', label ='Fitted data')
    plt.legend()
    plt.show()
    return (para[0])
print(fit())

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
 
def Fun(x,a1,a2,a3):                   # 定义拟合函数形式
    return a1*x**2+a2*x+a3
def error (p,x,y): # 拟合残
    return Fun(p,x)-y
def fit():
    x = x = np.array([1,3,4,5,6,7,8,9,10])
    p_value = [1,1,50] # 原始数据的参数,随便设定
    y = np.array([10,5,4,2,1,1,2,3,4])
    p0 = [0.1,-0.01,100] # 拟合的初始参数设置
    para,pcov=curve_fit(Fun,x,y)
    y_fitted = Fun(x,para[0],para[1],para[2]) # 画出拟合后的曲线
 
    plt.figure
    plt.plot(x,y,'r', label = 'Original curve')
    plt.plot(x,y_fitted,'-b', label ='Fitted curve')
    plt.legend()
    plt.show()
    return para
print(fit())

两者实现的结果和之前的都是一样的。如图所示。但是考虑的求解角度稍微不同，后面两者是基于scipy随便设定参数值一步一步通过迭代找最最优参数实现的。其中，

1. curve_fit

是在method =‘lm’（底层应该是lapack) 的情况下，该算法通过Minimumsq使用Levenberg-Marquardt算法。 注意，该算法只能处理不受约束的问题。

2. leastsq

leastsq 是MINPACK的lmdif和lmder算法的包装，cov_x是最小二乘目标函数的Hessian的雅可比近似。 该近似值假设目标函数基于某些观测目标数据（ydata）与参数f（xdata，params）的（非线性）函数之间的差 。