最近在看李航的统计学习方法P11时发现一个多项式函数拟合问题觉得公式的推导有问题，于是看了一些资料发现这里的推倒是有错误的，用python编程验证后发现按书上的求导结果拟合后的函数图像完全不对，下面给出正确的推导结果和对应的python实现与拟合效果。

（第一次写了一版这篇博客，公式自己敲的，但是从别人那里贴的代码，后面感觉那个代码思路没问题，但就是进行了很多次重复计算，数据量大的时候效率低，而且很多地方对python的语言特性利用不足，写的有点类似C语言风格，比较晦涩。于是乎，把代码重写了一次，准备把文章撤下来重新贴自己的代码，这时候发现没有保存，前面用Tex敲的一堆公式就没了，心中顿时万马奔腾。。。。好在今天考试考完，强行克服拖延症把这篇文章重写一遍并附上重写的代码。）

下面开始正文

问题描述

假定给定一个训练数据集：

$$T=\{(x_1,y_1),(x_2,y_2),\cdots,(x_N,y_N)\}$$

其中，$x_i\in R$是输入$x$的观测值，$y_i\in R$是相应的输出y的观测值，$i=1,2,\cdots,N$，多项式函数拟合的任务是假设给定数据由$M$次多项式函数生成，选择最有可能产生这些数据的$M$次多项式函数，即在$M$次多项式函数中选择一个对已知数据以及未知数据都有很好预测能力的函数。

设$M$次多项式为

$$f_M(x,w)=w_0+w_1x+w_2x^2+\cdots+w_Mx^M=\sum\limits_{j=0}^Mw_jx^j$$

式中$x$式单变量输入，$w_0,w_1,\cdots,w_m$是$M+1$个参数。

用平方损失作为损失函数，系数$\frac{1}{2}$是为了方便计算，将模型与训练数据代入，有

$$L(w)=\frac{1}{2}\sum\limits_{i=1}^{N} \big(\sum\limits_{j=0}^Mw_jx_i^j-y_i\big)^2$$

对$w_j$求偏导并令其为$0$

书里这里的求导是错误的，就不写错误的推导了，下面给出正确的推导

$$set \frac{\partial\ L(w)}{\partial {w_k}}=0\Rightarrow\\frac{1}{2}\sum \limits_{i=1}^{N}2\big(\sum\limits_{j=0}^Mw_jx_i^j-y_i\big)\times x_i^k=0\Rightarrow\sum \limits_{i=1}^{N}\sum\limits_{j=0}^Mw_jx_i^j＝\sum \limits_{i=1}^{N}x_i^ky_i(k=0,1,2,\cdots,M)$$

所以要求拟合多项式系数$w_0^*,w_1^*,\cdots,w_M^*$需要解下面这个线性方程组，下面的求和符号上下限都是$i=1$到$N$,为了方便略去不写。

$$\begin{bmatrix}N&\sum x_i &\sum x_i^2&\cdots& \sum x_i^M\\\sum x_i&\sum x_i^2&\sum x_i^3&\cdots&\sum x_i^{M+1}\\\sum x_i^2&\sum x_i^3&\sum x_i^4&\cdots&\sum x_i^{M+2}\\\vdots&\vdots&\vdots&\ddots&\cdots\\\sum x_i^M&\sum x_i^{M+1}&\sum x_i^{M+2}&\cdots&\sum x_i^{2M}\end{bmatrix}\begin{pmatrix}w_0\\w_1\\w_2\\\vdots\\w_m\end{pmatrix}=\begin{bmatrix}\sum y_i\\\sum x_iy_i\\\sum x_i^2y_i\\\vdots\\\sum x_i^My_i\end{bmatrix}$$

所以计算出$\sum\limits_{i=1}^{N} x_i^j(j=0,1,2,\cdots,2M)$和$\sum\limits_{i=1}^{N}x_i^jy_i(j=0,1,2,\cdots,M)$然后将这些值带入上述线性方程组求解即可。

下面给出python实现

# coding=utf-8

‘‘‘
作者:Xiaole Wen
程序:多项式曲线拟合算法
‘‘‘
import matplotlib.pyplot as plt
import math
import numpy
import random

fig = plt.figure()
ax = fig.add_subplot(111)

#在这里给出拟合多项式的阶数
order=9

#生成曲线上的各个点
x = numpy.arange(-1,1,0.02)
y = [((a*a-1)*(a*a-1)*(a*a-1)+0.5)*numpy.sin(a*2) for a in x]
#ax.plot(x,y,color=‘r‘,linestyle=‘-‘,marker=‘‘)
#,label="(a*a-1)*(a*a-1)*(a*a-1)+0.5"
plt.plot(x,y)
#生成的曲线上的各个点偏移一下，并放入到xa,ya中去
i=0
xa=[]
ya=[]
for xx in x:
    yy=y[i]
    d=float(random.randint(60,140))/100
    #ax.plot([xx*d],[yy*d],color=‘m‘,linestyle=‘‘,marker=‘.‘)
    i+=1
    xa.append(xx*d)
    ya.append(yy*d)

ax.plot(xa,ya,color=‘m‘,linestyle=‘‘,marker=‘.‘)
#存储从0次到m次的所有冥方和
bigMat=[]
for j in range(0,2*order+1):
    sum=0
    for i in range(0,len(xa)):
        sum+=(xa[i]**j)
    bigMat.append(sum)

#计算线性方程组系数矩阵
matA=[]
for rowNum in range(0,order+1):
    row=bigMat[rowNum:rowNum+order+1]
    matA.append(row)

matA=numpy.array(matA)

matB=[]
for i in range(0,order+1):
    ty=0.0
    for k in range(0,len(xa)):
        ty+=ya[k]*(xa[k]**i)
    matB.append(ty)

matB=numpy.array(matB)

matAA=numpy.linalg.solve(matA,matB)

#画出拟合后的曲线
#print(matAA)
xxa= numpy.arange(-1,1.06,0.01)
yya=[]
for i in range(0,len(xxa)):
    yy=0.0
    for j in range(0,order+1):
        dy=(xxa[i]**j)
        dy*=matAA[j]
        yy+=dy
    yya.append(yy)
ax.plot(xxa,yya,color=‘g‘,linestyle=‘-‘,marker=‘‘)

ax.legend()
plt.show()

下面给出阶叔分别取3和取9的时候的拟合结果

图中蓝色的线代表原始数据生成函数，绿色代表拟合函数