【统计学习代码】第一章

1.导读

在这个文章中主要是统计学习方法第二版第一章的代码，主要有使用最小二乘法曲线和正则化

2.最小二乘法

目标函数是sin(2*pi*x) 加上一个正态分布的噪音干扰，用多项式去拟合，具体代码如下：

import numpy as np
import scipy as sp
from scipy.optimize import  leastsq
import matplotlib.pyplot as plt

导入包 第一行 是导入numpy;NumPy（Numerical Python）是Python的一种开源的数值计算扩展 第二行 导入scipy包 ; scipy包含致力于科学计算中常见问题的各个工具箱 第三行 导入最小二乘， 第四行 是画图的库matplotlib ps: numpy.poly1d([1,2,3]) 生成1x^0+2x^1+3x^0

 

 

 

# 目标函数
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

real_func 是目标函数 fit_func 是多项式 resiudals_func 是误差，或者是残差

 

 

# 十个点
x = np.linspace(0, 1, 10)
x_points = np.linspace(0, 1, 1000)
# 加上正态分布噪音的目标函数的值
y_ = real_func(x)
y = [np.random.normal(0, 0.1) + y1 for y1 in y_]

 

def fitting(M=0):
    """
    M    为 多项式的次数
    """
    # 随机初始化多项式参数
    p_init = np.random.rand(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='fitted curve')
    plt.plot(x, y, 'bo', label='noise')
    plt.legend()
    return p_lsq

# M=0
p_lsq_0 = fitting(M=0)

  

# M=1
p_lsq_1 = fitting(M=1)

 

# # M=3
p_lsq_3 = fitting(M=3)

 

# M=9
p_lsq_9 = fitting(M=9)

可以看出M=9时候造成了过拟合 

 

 3.正则化

结果显示过拟合，引入正则化项，降低过拟合

 

• L1: regularization*abs(p)
• L2: 0.5 * regularization * np.square(p)

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)))  # L2范数作为正则化项
    return ret

# 最小二乘法,加正则化项
p_init = np.random.rand(9 + 1)
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[0], 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()