最小二乘法多项式曲线拟合及其python实现¶

import numpy as np
import math
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
%matplotlib inline
path = '/usr/share/fonts/truetype/SimHei.ttf'
fontprop = fm.FontProperties(fname=path, size=13)

SAMPLE_NUM=100#要生成的sample个数
M=9 #多项式阶数

def Fun(N, deg):
  ''' 产生原曲线 sin函数
  '''
  #产生带有高斯噪声的信号
  mid, sigma = 0, 0.3 # 设置均值和方差
  noise = np.random.normal(mid, sigma, N).reshape(N,1) #生成SAMPLE_NUM个数据
  #产生SAMPLE_NUM个序号(范围是2pi)
  x = np.arange(0, N).reshape(N,1)/(N-1)*(2*math.pi)
  #generate y and y_noise, and both y's and y_noise's shape is (SAMPLE_NUM*1)
  y=np.sin(x)
  y_noise=np.sin(x)+noise
  return (x,y,y_noise)

x, y, y_noise = Fun(SAMPLE_NUM, M)
#绿色曲线显示x - y，散点显示x - y_noise
f1, = plt.plot(x,y,'g',lw=4.0) # , label = r'$y=\sin(x)$'
f2, = plt.plot(x,y_noise,'bo') # ,label = r'$y=\sin(x)$ with noise'

#generate Matrix X which has M order
X=x
for i in range(2,M+1):
         X = np.column_stack((X, pow(x,i)))

#add 1 on the first column of X, now X's shape is (SAMPLE_NUM*(M+1))
X = np.insert(X,0,[1],1)
#print(X)

#calculate W, W's shape is ((M+1)*1)#
#W=np.linalg.inv((X.T.dot(X))).dot(X.T).dot(y_noise)#have no regularization
W=np.linalg.inv((X.T.dot(X))+np.exp(-8) * np.eye(M+1)).dot(X.T).dot(y_noise)#introduce regularization
y_estimate=X.dot(W)

#红色曲线显示x - y_estimate
f3, = plt.plot(x,y_estimate,'r',lw=4.0)
plt.suptitle(r'最小二乘法多项式拟合原函数 $y=\sin(x)$', fontsize=10,fontproperties=fontprop)
plt.grid(color='gray')
plt.grid(linewidth='1')
plt.grid(linestyle='--')
plt.legend([f1,f2,f3],[r'原函数 $y=\sin(x)$',r'$y=\sin(x)$ 叠加噪音',r'多项式拟合'],fontsize=10,prop=fontprop)
plt.show()

import numpy as np
import math
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
from scipy.optimize import leastsq
%matplotlib inline
path = '/usr/share/fonts/truetype/SimHei.ttf'
fontprop = fm.FontProperties(fname=path, size=13)

SAMPLE_NUM=100#要生成的sample个数
M=9 #多项式阶数

def Fun(N, deg):
  ''' 产生原曲线 sin函数
  '''
  #产生带有高斯噪声的信号
  mid, sigma = 0, 0.3 # 设置均值和方差
  noise = np.random.normal(mid, sigma, N).reshape(N,1) #生成SAMPLE_NUM个数据
  #产生SAMPLE_NUM个序号(范围是2pi)
  x = np.arange(0, N).reshape(N,1)/(N-1)*(2*math.pi)
  #generate y and y_noise, and both y's and y_noise's shape is (SAMPLE_NUM*1)
  y=np.sin(x)
  y_noise=np.sin(x)+noise
  return (x,y,y_noise)

def Error(p, x, y):
  '''定义误差函数
  例如输入：[1,2,3]
  返回：x^2 + 2x + 3
  多项式的次数是从0开始记的，要注意这个地方
  '''
  fun = np.poly1d(p)    # poly1d（）函数可以按照输入的列表p返回一个多项式函数
  return (y - fun(x)).flatten()   # 返回真实值 与我们拟合的曲线上对应的值的差

def fit(deg, x, y):
    '''拟合函数  
    返回：注意返回值中的第一行才是拟合曲线的参数列表 para = fit(deg, x, y)[0]
    '''
    pars = np.random.rand(deg+1)  # 生成p+1个随机数的列表，这样poly1d函数返回的多项式次数就是p
    r = leastsq(Error, pars, args=(x, y))   # 三个参数：误差函数、函数参数列表、数据点
    return r

x, y, y_noise = Fun(SAMPLE_NUM, M)
#绿色曲线显示x - y，散点显示x - y_noise
f1, = plt.plot(x,y,'g',lw=4.0) # , label = r'$y=\sin(x)$'
f2, = plt.plot(x,y_noise,'bo') # ,label = r'$y=\sin(x)$ with noise'

W=fit(M, x, y_noise)[0]
y_estimate = np.poly1d(W)(x)
#红色曲线显示x - y_estimate
f3, = plt.plot(x,y_estimate,'r',lw=4.0)

# 使用 numpy.polynomial.polynomial.polyfit
from numpy.polynomial import polynomial as P
c, stats = P.polyfit(x.flatten(),y_noise,M,full=True)
y_np_estimate = np.poly1d(c.flatten()[::-1])(x)
#黄色曲线显示x - y_np_estimate, 注意，为能观测到曲线， y_np_estimate上移了0.1
f4, = plt.plot(x,y_np_estimate+0.1,'y',lw=4.0)

plt.suptitle(r'最小二乘法多项式拟合原函数 $y=\sin(x)$', fontsize=10,fontproperties=fontprop)
plt.grid(color='gray')
plt.grid(linewidth='1')
plt.grid(linestyle='--')
plt.legend([f1,f2,f3,f4],[r'原函数 $y=\sin(x)$',r'$y=\sin(x)$ 叠加噪音','Scipy leastsq拟合','np polyfit拟合'],fontsize=10,prop=fontprop)
plt.show()