题目论述
用下面的模型产生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