最小二乘法拟合正弦函数

$y=sin2\pi x$

$y+ϵ$, $ϵ$~N(0,1)
#给y加上高斯分布的噪点

$(X,Y)$

attempt to fitting the function : $y=a_0+a_{1}x+a_2{x}^2+..+a_n{x}^n$

#导入包
import numpy as np
import scipy as sp
from scipy.optimize import leastsq
import matplotlib.pyplot as plt

#定义正弦函数
def sine_func(x):
    return np.sin(2*np.pi*x)
#定义多项式函数
def fit_func(p,x):    #p为数组
    f=np.poly1d(p)    #poly1d函数能根据输入的数组p生成一个多项式
    return f(x)

np.poly1d([1,2,3])

$1\ast x^2+2\ast x^1+3\ast x^0$

policy: ${\min }_a(y-(ax+b{x}^0){)}^2$

#评价差 error
def error_func(p,x,y):
    error=fit_func(p,x)-y   # error=将x值带入多项式后减去y值
    return error            #返回真实值与我们拟合的曲线上对应的值的差
x=np.linspace(0,1,10)  #linspace函数产生0-1之间的步长为0.1的数组
y_=sine_func(x)        
y=[np.random.normal(0,0.1)+y1 for y1 in y_] #得到加上噪点后的y值
print(y)

[0.018863818084578765, 0.643352435292383, 0.9995270889071429, 0.8859076268379011, 0.3123492536633993, -0.2445361290110839, -0.857851035405584, -0.7992334756421094, -0.6318390256462963, -0.15169614809124077]

x_points=np.linspace(0,1,1000)
#leastsqura..
def fitting(M=0):                 #得到拟合M次的曲线
    p_init=np.random.rand(M+1)    #生成M+1个服从[0,1)均匀分布的随机样本值
    p_lsq=leastsq(error_func,p_init,args=(x, y))  # 最小二乘法函数 三个参数：误差函数、函数参数列表、数据点
    print('parameters:', p_lsq[0])              #p_lsp函数返回值的第一个为拟合曲线的参数
    plt.plot(x_points,sine_func(x_points),label='real')  # real曲线
    plt.plot(x_points,fit_func(p_lsq[0],x_points),label='fitting') #拟合曲线  
    plt.plot(x,y,'bo', label = 'noise')    #训练值样本点   'b0'画出蓝色的圆圈
    plt.legend()
    return p_lsq
fitting(M=0)

parameters: [0.01748444]
(array([0.01748444]), 3)

fitting(M=1)

parameters: [-1.37690968  0.70593928]
(array([-1.37690968,  0.70593928]), 1)

fitting(M=5)

parameters: [-5.67180721e+01  1.37260155e+02 -9.56060001e+01  9.15613970e+00
  5.73126111e+00  1.54421206e-02]
(array([-5.67180721e+01,  1.37260155e+02, -9.56060001e+01,  9.15613970e+00,
  5.73126111e+00,  1.54421206e-02]),1)

fitting(M=9)

parameters: [ 2.44048348e+04 -1.07709860e+05  1.99353327e+05 -2.00963213e+05
  1.19785217e+05 -4.27681566e+04  8.80453207e+03 -9.53138175e+02
  4.62865992e+01  1.88638181e-02]
(array([ 2.44048348e+04, -1.07709860e+05,  1.99353327e+05, -2.00963213e+05,
  1.19785217e+05, -4.27681566e+04,  8.80453207e+03, -9.53138175e+02,
  4.62865992e+01,  1.88638181e-02]),2)

policy: ${\min }_a(y-(ax+b{x}^0){)}^2+\lambda \Vert a{\Vert }_2$

penalty parameter: $\lambda$

regularization: $\Vert a{\Vert }_2$

lam=0.001
def error_plus_regu(p,x,y):
    error=y-fit_func(p,x)
    error=np.append(error,np.sqrt(0.5*lam*np.square(p)))   #sqrt函数 求算数平方根  square函数 平方
    return error 

M=3
p_init=np.random.rand(M+1) 
p_lsq = leastsq(error_func, p_init, args=(x,y))
p_lsq_regu=leastsq(error_plus_regu,p_init,args=(x,y))  #添加正则项后的参数学习结果：p_lsq_regu
plt.plot(x_points,sine_func(x_points),label='real')   #follow sin函数的plot
plt.plot(x_points,fit_func(p_lsq_regu[0],x_points),label='regulation') #follow lasso的plot
plt.plot(x_points,fit_func(p_lsq[0],x_points),label='fitting')  # 原始fitting plot
plt.plot(x,y,'bo', label = 'noise')#训练值样本点
plt.legend()