最小二乘拟合-python scipy实现

 

题目论述

 

用下面的模型产生50个数，构成一个序列：

$$y = \frac{{\sin (\pi x)}}{{\pi x}} + 0.1x + 0.05r{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} x \in \left[ { - 3,3} \right]$$1-1

其中r为[0,1]分布的均匀随机数。

使用下列基函数对该模型进行近似：$${\bf{\varphi }}\left( x \right) = \left[ {1,\sin \frac{x}{2},\cos \frac{x}{2},\sin \frac{{2x}}{2},\cos \frac{{2x}}{2}, \cdots ,\sin \frac{{15x}}{2},\cos \frac{{15x}}{2}} \right]$$

1-2

 

原理分析

 

这是一个函数拟合问题，已知离散数据点，要求使用给定的基函数进行拟合，最后得到基函数的系数向量。基函数有多种形式，比如多项式、指数、对数、三角函数等。常见的计算工具中提供了多项式拟合的函数，如matlab中的polyfit(x,y,n), python的numpy.polyfit(x,y,n); 而对于非线性拟合，matlab中可以使用fittype()和fit(), python scipy库中有curve_fit(). 但针对本题，基函数有31项， 实测curve_fit()无法收敛。事实上，各种基函数的拟合问题都可以归类为同一种问题，使用同样的算法求解基函数系数。

评价拟合结果的好坏，应当是拟合所得函数对原始数据的逼近程度，并将以此为约束，寻求最适合的基函数系数。定量描述误差的参数有：

偏差绝对值之和：\(\sum\limits_{i = 1}^n {\left| {{f_{\bf{\theta }}}\left( {{x_i}} \right) - {y_i}} \right|} \)

最大的偏差绝对值：\(\max \left( {\left| {{f_{\bf{\theta }}}\left( {{x_i}} \right) - {y_i}} \right|} \right)\)

偏差的平方和：\(\sum\limits_{i = 1}^n {{{\left( {{f_{\bf{\theta }}}\left( {{x_i}} \right) - {y_i}} \right)}^2}} \)

其中，\({y_i}\)为第i个原始数据，\({f_{\bf{\theta }}}\left( {{x_i}} \right)\)是在基函数系数\({\bf{\theta }}\)下拟合得到的方程，\({x_i} = {y_i}\)，n是原始数据的总长度。

其中，第三种约束条件最为常用，称为最小二乘法。在该约束条件下的基函数系数向量\({\bf{\theta }}\)的求解推导如下：

平方误差可以写为：\({J_{LS}}\left( {\bf{\theta }} \right) = \frac{1}{2}\sum\limits_{i = 1}^n {\left( {{f_{\bf{\theta }}}({x_i}) - {y_i}} \right)} \) 2-1

用矩阵可以表示为：

\({J_{LS}}\left( {\bf{\theta }} \right) = \frac{1}{2}{\left| {{\bf{D\theta }} - {\bf{y}}} \right|^2}\)   2-2

其中，\({\bf{y}} = {[{y_1},{y_2}, \cdots ,{y_n}]^T}\)，为原始数据的向量，\({\bf{\theta }} = {[{\theta _1},{\theta _2}, \cdots ,{\theta _b}]^T}\)，为基函数的系数向量，

，是由基函数\({\bf{\varphi }}(x)\)中, x取n个不同的值\({x_i} = {y_i}\)，所构成的n*b维矩阵，称之为设计矩阵。

    对式2-2求\({\bf{\theta }}\)的偏导，可得：$${\nabla _{\bf{\theta }}}{J_{LS}} = \left( {\frac{{\partial {J_{LS}}}}{{\partial {\theta _1}}},\frac{{\partial {J_{LS}}}}{{\partial {\theta _2}}}, \cdots ,\frac{{\partial {J_{LS}}}}{{\partial {\theta _b}}}} \right) = {{\bf{D}}^T}{\bf{D\theta }} - {{\bf{D}}^T}{\bf{y}}$$

2-3

当\({\nabla _{\bf{\theta }}}{J_{LS}}\)=0时，认为此时的系数\({\bf{\theta }}\)可以使估计偏差最小，所以\({{\bf{D}}^T}{\bf{D\theta }} = {{\bf{D}}^T}{\bf{y}}\)。当\({\bf{D}}\)可逆时，

\({\bf{\theta }} = {{\bf{D}}^T}{\bf{y}}\)    2-4

这就是求基函数系数的方法。一旦求出\({\bf{\theta }}\)后,\({f_{\bf{\theta }}}(x)\) 就可以表示为：

$${f_{\bf{\theta }}}(x) = {\bf{\varphi }}(x)$$  2-5

在本题中，基函数为三角函数，原始数据长度n=50，基函数项数b=31因而设计矩阵为50×31维，即：

2-6

 

程序框图

 

根据上节所述，要计算基函数的系数，首先要构造2-6式所述的设计矩阵D，然后计算2-4表达式，就可以得到基函数的系数向量。程序框图表示如下：

 

Fig. 31程序框图

 

编程实现

 

 

运行环境

 

Python2.7.12, 32Bits

Numpy(1.11.2+mkl)

Matplotlib(1.5.3)

 

变量定义与含义说明

n=50#原始数据的个数

b=31#基向量个数

designMatrix=np.ones((n,b)) #设计矩阵，50×31，初始化为1

xIndex=np.linspace(-3,3,n) #产生了-3～3的50个数作为x轴

yIndex=np.linspace(-3,3,500)#产生了-3～3的500个数作为y轴，用于拟合后的函数绘图

 

yRaw：存放50个原始数据的y值

coef：计算所得的基函数系数

yPoly：存放拟合函数fx的y值

 

框图1，生成原始数据

 

根据1-1式，首先定义了一个函数 dataGener(x)，然后使用numpy的vectorize()函数对其矢量化，命名为Func_RawArray(). dataGener(x)原本是python中的普通函数，输入参数是标量。矢量化后得到的Func_RawArray()函数，输入变量是array对象，可以实现向量输入，向量返回。调用Func_RawArray(xIndex)，就可以得到50个原始数据，存入yRaw。这里[0,1]分布的随机数使用numpy.random.random得到。

这部分的全部代码如下：

def dataGener (x): #定义原始数据产生的函数
    return np.sin(np.pi*x)/(np.pi*x)+0.1*x+0.05*np.random.random(size=1)
  
#使用vectorize矢量化函数，进而计算产生原始数据的向量,存入yRaw
Func_RawArray=np.vectorize(dataGener)
yRaw=Func_RawArray(xIndex)

 

框图2，构建设计矩阵

 

在一开始时，已经初始化了n*b的矩阵designMatrix，用于保存设计矩阵。设计矩阵由2-6时表示，使用2个循环分别计算sin项和cos项即可。

构建设计矩阵的代码为：

#构建设计矩阵
for i in range (1,b,2):#sin 项
    designMatrix[:,i]=np.sin((i+1.0)/2* xIndex/2 )

for i in range (2,b,2): #cos 项
    designMatrix[:,i]=np.cos(i/2.0*xIndex/2 )

 

框图3，计算基函数系数

 

基函数的系数coef计算使用式2-4. 由于numpy的array对象没有方便的求逆方法，所以使用numpy.matrix()函数将array转换为matrix。matrix.I可以方便得求逆矩阵。为了矩阵运算的维度对应，yRaw用reshape转换维度后也转换到了matrix. 计算得到系数后顺便打印到屏幕。这部分的代码如下：

#计算基函数系数
coef=np.matrix(designMatrix).I *np.matrix(yRaw.reshape(50,1))
print '基函数的系数为'
print coef #打印基函数的系数

 

框图4，根据所得基函数系数，得到拟合函数fx

 

基函数为1-2式，即：

其系数已经计算得到，存储在coef中。所以对于任意的x，都应当可以计算出拟合结果.

    为了验证拟合效果，x的范围也取-3～3，但点数足够多，取500点。一开始已经初始化了这500个样点的横坐标，存在yIndex中。

为计算，定义了函数fxGen()，输入形参为横坐标的向量x，基函数矩阵coef和基函数项数b。使用该函数就可以计算出对于的y值向量。最后 y值的size也是500点，存入到yPoly中。代码如下：

def fxGen(x,coef=coef,b=b):
    #x为array类型的f(x)横坐标向量，使用linspace产生, coef为计算所得
    #的基函数系数，b为基函数项数
    fx=np.zeros(x.shape[0]) #初始化和x同样长度的fx
    fx=coef[0,0]*1   #此处是1!!!,找了七八个小时的bug :(
    for i in range(1,b,2):
        fx+=coef[i,0]*np.sin((i+1.0)/2* x/2 )
        
    for i in range (2,b,2): #cos 项
        fx+=coef[i,0]*np.cos(i/2.0*x /2 )
    
    # for i in range(1,(b+1)/2): #i从1～15
        # fx+=coef[2*i-1,0]*np.sin(i* x /2)
        # fx+=coef[2*i,0]*np.cos(i* x /2)
    return fx
    
#yPoly=designMatrix*coef
yPoly=fxGen (yIndex,coef,b)

 

框图5，绘图

 

绘图使用matplotlib.pyplot即可。

 

plot(xIndex,yRaw,'r.')
plot(yIndex,yPoly)
plt.xlabel('x')
plt.ylabel('y')
show()

 

运行结果

 

Fig. 41 运行结果

 

问题与总结

之前使用过matlab的polyfit作线性拟合，但不知原理。这次本来打算再用现成的函数，得知scipy中curve_fit()可以作非线性拟合，示例中展示的指数拟合，对数拟合效果也确实不错。但这个题目的基函数有31项，使用curve_fit()就不能收敛了，只能计算出前几项的系数，后面的都是0，拟合结果如图5-1，并不理想。后来在网上看到了最小二乘约束条件下的拟合方法，才明白在同样约束条件下，基函数的形式其实不重要，因为算法都是相同的。计算方法就是基于平方误差最小这个条件，对平方误差求系数向量的偏微分，偏微分等于0的解就是所求的系数。至此又学到了新的知识。

    经过尝试，我还认为，基函数也应当根据原始数据适当选取。在本题中，基函数使用，在该基函数下，计算所得的系数很大，在10^8数量级(见fig4-1)，而且受到随机误差的影响很大，两次计算，只是因为0.05的随机误差r，同一项的系数甚至会改变正负号，因而计算结果鲁棒性欠缺，会出现如图5-1中所表现出来的震荡现象，即在数据起始和末尾阶段，误差会突然变大。

而如果使用基函数，其系数在10^-1左右，很接近原始数据的数量级，稳定性就比较好。

    

 

Fig. 51 使用curve_fit()求解的结果并不好

Reference

--------------------------------------------------------------

[1] JairusChan, 最小二乘法多项式曲线拟合原理与实现, 2012-04-27

http://blog.csdn.net/jairuschan/article/details/7517773

[1] lfdanding, 最小二乘法曲线拟合原理与实现, 2016-02-25

http://blog.csdn.net/lfdanding/article/details/50739987