多项式函数曲线拟合——最小二乘法

        多项式函数拟合的任务是假设给定数据由M次多项式函数生成，选择最有可能产生这些数据的M次多项式函数，即在M次多项式函数中选择一个对已知数据以及未知数据都有很好预测能力的函数。

        最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。

设训练数据集为： 

                           
其中，是输入x的观测值，是相应的输出y的观测值，i=1,2,⋯,N

设M次多项式为

                      
其中，x是单变量输入，是M+1个参数。

用平方损失作为损失函数（即最小二乘法），系数是为了方便计算，将模型与训练数据代入，有：

                                              
对wj求偏导并令其为0
                     

以上公式最后一步存在一处错误：等式左边x指数为（j+k）

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

                      
所以计算出

                

然后将这些值带入上述线性方程组求解即可。

Python实现：

'''
   多项式：yi = w0 + w1*xi^1 + w2*xi^2 + ... + wn*xi^m
   N为数据点个数，M为阶数
   先用数据点（xa、ya）求出未知参数，然后用参数带入后的公式求解给定值(xxa)
'''
import matplotlib.pyplot as plt
import numpy
import random

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

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

# 1. 生成曲线上的各个点
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]   

# y=[(x^2-1)^3]*sin(2x)，阶数为6？？？

# ax.plot(x,y,color='r',linestyle='-',marker='')
# ,label="(a*a-1)*(a*a-1)*(a*a-1)+0.5"
plt.plot(x,y,c='red')

# 2. 生成的曲线上的各个点偏移一下，并放入到x_data,y_data中去
i=0
x_data=[]
y_data=[]
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
    x_data.append(xx*d)
    y_data.append(yy*d)

ax.plot(x_data,y_data,color='m',linestyle='',marker='.')

# 3. 计算Ax=b中，矩阵A、b
# 存储从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)

# 计算线性方程组系数矩阵：A
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)

matW=numpy.linalg.solve(matA,matB)
# numpy.linalg中的函数solve可以求解形如 Ax = b 的线性方程组，其中 A 为矩阵，b 为一维或二维的数组，x 是未知变量

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

ax.legend()
plt.show()

拟合结果曲线：

      阶数=9

                   

   阶数 = 5 ：

                  

阶数 = 3：