李航《统计学方法chapter one 实践》

import  numpy as np
import  scipy as sp
from scipy.optimize import  leastsq
import matplotlib.pyplot as plt
#目标函数
def real_func(x):
    return np.sin(2*np.pi*x)
#多项式
def fit_func(p,x):
    f=np.poly1d(p)
    return f(x)

def residuals_func(p,x,y):
    ret=fit_func(p,x)-y
    return ret

# 十个点
x=np.linspace(0,1,10)
x_points = np.linspace(0,1,1000)
y_=real_func(x)
#np.random.normal(0,0.1)(噪音符合正太分布)
y =[np.random.normal(0,0.1)+y1 for y1 in y_]

def fitting (M=0):
    # 随机初始化多项式参数
    #返回M+1个在(0,1)之间的值
    p_init= np.random.rand(M+1)
    #最小二乘法 least_sq返回值第一个是(M+1)个系数的结果
    p_lsq=leastsq(residuals_func,p_init,args=(x,y))
    print("fitting parameters:",p_lsq[0])

    #可视化
    plt.plot(x_points,real_func(x_points),label="real")
    plt.plot(x_points,fit_func(p_lsq[0],x_points),label="fit")
    plt.plot(x,y,'bo',label="noise")
    #legend()用于显示 label
    plt.show()
    plt.legend()
    return p_lsq[0]

p_lsq_0 = fitting(M=0)
p_lsq_1 = fitting(M=1)
p_lsq_3 = fitting(M=3)
p_lsq_9 = fitting(M=9)

regularization  = 0.0001
def residuals_func_regularization(p,x,y):
    ret=fit_func(p,x)-y
    #数组拼接    注意这里正则项要开根
    ret=np.append(ret,np.sqrt(0.5*regularization*np.square(p)))
    return ret
p_init = np.random.rand(9)
p_lsq_regularization = leastsq(residuals_func_regularization,p_init,args=(x,y))
plt.plot(x_points,real_func(x_points),label="real")
plt.plot(x_points,fit_func(p_lsq_9,x_points),label="fitted curve")
plt.plot(x_points,fit_func(p_lsq_regularization[0],x_points),label="regularization")
plt.plot(x,y,"bo",label="noise")
plt.legend()
plt.show()