第1章 统计学习方法概论（第2版）

使用最小二乘法拟和曲线

高斯于1823年在误差 e 1 , … , e n e_1,…,e_n e1​,…,en​独立同分布的假定下,证明了最小二乘方法的一个最优性质: 在所有无偏的线性估计类中,最小二乘方法是其中方差最小的！
对于数据 ( x i , y i ) ( i = 1 , 2 , 3... , m ) (x_i, y_i) (i=1, 2, 3...,m) (xi​,yi​)(i=1,2,3...,m)

拟合出函数 h ( x ) h(x) h(x)

有误差，即残差： r i = h ( x i ) − y i r_i=h(x_i)-y_i ri​=h(xi​)−yi​

此时 L 2 L2 L2范数(残差平方和)最小时， h ( x ) h(x) h(x) 和 y y y 相似度最高，更拟合

一般的 H ( x ) H(x) H(x)为 n n n次的多项式， H ( x ) = w 0 + w 1 x + w 2 x 2 + . . . w n x n H(x)=w_0+w_1x+w_2x^2+...w_nx^n H(x)=w0​+w1​x+w2​x2+...wn​xn

w ( w 0 , w 1 , w 2 , . . . , w n ) w(w_0,w_1,w_2,...,w_n) w(w0​,w1​,w2​,...,wn​)为参数

最小二乘法就是要找到一组 w ( w 0 , w 1 , w 2 , . . . , w n ) w(w_0,w_1,w_2,...,w_n) w(w0​,w1​,w2​,...,wn​) ，使得 ∑ i = 1 n ( h ( x i ) − y i ) 2 \sum_{i=1}^n(h(x_i)-y_i)^2 ∑i=1n​(h(xi​)−yi​)2 (残差平方和) 最小即，求 m i n ∑ i = 1 n ( h ( x i ) − y i ) 2 min\sum_{i=1}^n(h(x_i)-y_i)^2 min∑i=1n​(h(xi​)−yi​)2

例题1.1 P21
使用0~9次多项式函数对数据进行拟合，我们用目标函数 y = s i n 2 π x y=sin2{\pi}x y=sin2πx, 加上一个正态分布的噪音干扰，用多项式去拟合

一、引入所需要的相关库

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

二、创建所需要的函数

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

三、方法代码主体

# 在0~1之间均匀生成十个点
x = np.linspace(0, 1, 10)# linspace(start,end,num)
# print(x)
# 在0~1之间均匀生成一千个点
x_points = np.linspace(0, 1, 1000)
#print(x_points)
# 加上正态分布噪音的目标函数的值
y_ = real_func(x)
# print(y_)
# 先使用for循环读出y_每一个元素，赋值给y1,之后让y1与均值为0，方差为1的随机数进行相加，最终得到数组y
y = [np.random.normal(0, 0.1) + y1 for y1 in y_]


def fitting(M=0):
    """
    M    为 多项式的次数
    """
    # 随机初始化多项式参数
    # no.random.rand()通过本函数可以返回一个或一组服从“0~1”均匀分布的随机样本值。随机样本取值范围是[0,1)，不包括1。
    p_init = np.random.rand(M + 1)
    print('返回M+1个0-1之间的随机数：',p_init)
    # 最小二乘法
    p_lsq = leastsq(residuals_func, p_init, args=(x, y))#最小二乘函数包含三个参数：
    #（1）func：误差函数 （2）x0：表示函数的参数 （3）args（）表示数据点
    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时

# M=0 调用fitting函数
p_lsq_0 = fitting(M=0)

实验结果：

M=1时

# M=1 调用fitting函数
p_lsq_1 = fitting(M=1)

实验结果：

M=3时

# M=3 调用fitting函数
p_lsq_3 = fitting(M=3)

实验结果：

M=9时

# M=9 调用fitting函数
p_lsq_9 = fitting(M=9)

实验结果：

当M=9时，多项式曲线通过了每个数据点，但是造成了过拟合

五、使用正则化

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

Q ( x ) = ∑ i = 1 n ( h ( x i ) − y i ) 2 + λ ∣ ∣ w ∣ ∣ 2 Q(x)=\sum_{i=1}^n(h(x_i)-y_i)^2+\lambda||w||^2 Q(x)=∑i=1n​(h(xi​)−yi​)2+λ∣∣w∣∣2。

回归问题中，损失函数是平方损失，正则化可以是参数向量的L2范数,也可以是L1范数。
定义残差函数正则化方法

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()

实验结果

六、课后习题答案（P33）

习题1.1
  说明伯努利模型的极大似然估计以及贝叶斯估计中的统计学习方法三要素。伯努利模型是定义在取值为0与1的随机变量上的概率分布。假设观测到伯努利模型 n n n次独立的数据生成结果，其中 k k k次的结果为1，这时可以用极大似然估计或贝叶斯估计来估计结果为1的概率。
解答：
伯努利模型的极大似然估计以及贝叶斯估计中的统计学习方法三要素如下：

1. 极大似然估计
   模型： F = { f ∣ f p ( x ) = p x ( 1 − p ) ( 1 − x ) } \mathcal{F}=\{f|f_p(x)=p^x(1-p)^{(1-x)}\} F={f∣fp​(x)=px(1−p)(1−x)}
   策略： 最大化似然函数
   算法： arg ⁡ min ⁡ p L ( p ) = arg ⁡ min ⁡ p ( n k ) p k ( 1 − p ) ( n − k ) \displaystyle \mathop{\arg\min}_{p} L(p)= \mathop{\arg\min}_{p} \binom{n}{k}p^k(1-p)^{(n-k)} argminp​L(p)=argminp​(kn​)pk(1−p)(n−k)
2. 贝叶斯估计
   模型： F = { f ∣ f p ( x ) = p x ( 1 − p ) ( 1 − x ) } \mathcal{F}=\{f|f_p(x)=p^x(1-p)^{(1-x)}\} F={f∣fp​(x)=px(1−p)(1−x)}
   策略： 求参数期望
   算法：
   E π [ p ∣ y 1 , ⋯   , y n ] = ∫ 0 1 p π ( p ∣ y 1 , ⋯   , y n ) d p = ∫ 0 1 p f D ( y 1 , ⋯   , y n ∣ p ) π ( p ) ∫ Ω f D ( y 1 , ⋯   , y n ∣ p ) π ( p ) d p d p = ∫ 0 1 p k + 1 ( 1 − p ) ( n − k ) ∫ 0 1 p k ( 1 − p ) ( n − k ) d p d p \begin{aligned} E_\pi\big[p \big| y_1,\cdots,y_n\big] & = {\int_0^1}p\pi (p|y_1,\cdots,y_n) dp \ & = {\int_0^1} p\frac{f_D(y_1,\cdots,y_n|p)\pi(p)}{\int_{\Omega}f_D(y_1,\cdots,y_n|p)\pi(p)dp}dp \ & = {\int_0^1}\frac{p^{k+1}(1-p)^{(n-k)}}{\int_0^1 p^k(1-p)^{(n-k)}dp}dp \end{aligned} Eπ​[p∣∣​y1​,⋯,yn​]​=∫01​pπ(p∣y1​,⋯,yn​)dp=∫01​p∫Ω​fD​(y1​,⋯,yn​∣p)π(p)dpfD​(y1​,⋯,yn​∣p)π(p)​dp=∫01​∫01​pk(1−p)(n−k)dppk+1(1−p)(n−k)​dp​

伯努利模型的极大似然估计：
定义 P ( Y = 1 ) P(Y=1) P(Y=1)概率为 p p p，可得似然函数为： L ( p ) = f D ( y 1 , y 2 , ⋯   , y n ∣ θ ) = ( n k ) p k ( 1 − p ) ( n − k ) L(p)=f_D(y_1,y_2,\cdots,y_n|\theta)=\binom{n}{k}p^k(1-p)^{(n-k)} L(p)=fD​(y1​,y2​,⋯,yn​∣θ)=(kn​)pk(1−p)(n−k)方程两边同时对 p p p求导，则： 0 = ( n k ) [ k p k − 1 ( 1 − p ) ( n − k ) − ( n − k ) p k ( 1 − p ) ( n − k − 1 ) ] = ( n k ) [ p ( k − 1 ) ( 1 − p ) ( n − k − 1 ) ( m − k p ) ] \begin{aligned} 0 & = \binom{n}{k}[kp^{k-1}(1-p)^{(n-k)}-(n-k)p^k(1-p)^{(n-k-1)}]\ & = \binom{n}{k}[p^{(k-1)}(1-p)^{(n-k-1)}(m-kp)] \end{aligned} 0​=(kn​)[kpk−1(1−p)(n−k)−(n−k)pk(1−p)(n−k−1)]=(kn​)[p(k−1)(1−p)(n−k−1)(m−kp)]​可解出 p p p的值为 p = 0 , p = 1 , p = k / n p=0,p=1,p=k/n p=0,p=1,p=k/n，显然 P ( Y = 1 ) = p = k n \displaystyle P(Y=1)=p=\frac{k}{n} P(Y=1)=p=nk​

伯努利模型的贝叶斯估计：
定义 P ( Y = 1 ) P(Y=1) P(Y=1)概率为 p p p， p p p在 [ 0 , 1 ] [0,1] [0,1]之间的取值是等概率的，因此先验概率密度函数 π ( p ) = 1 \pi(p) = 1 π(p)=1，可得似然函数为： L ( p ) = f D ( y 1 , y 2 , ⋯   , y n ∣ θ ) = ( n k ) p k ( 1 − p ) ( n − k ) L(p)=f_D(y_1,y_2,\cdots,y_n|\theta)=\binom{n}{k}p^k(1-p)^{(n-k)} L(p)=fD​(y1​,y2​,⋯,yn​∣θ)=(kn​)pk(1−p)(n−k)
根据似然函数和先验概率密度函数，可以求解 p p p的条件概率密度函数： π ( p ∣ y 1 , ⋯   , y n ) = f D ( y 1 , ⋯   , y n ∣ p ) π ( p ) ∫ Ω f D ( y 1 , ⋯   , y n ∣ p ) π ( p ) d p = p k ( 1 − p ) ( n − k ) ∫ 0 1 p k ( 1 − p ) ( n − k ) d p = p k ( 1 − p ) ( n − k ) B ( k + 1 , n − k + 1 ) \begin{aligned}\pi(p|y_1,\cdots,y_n)&=\frac{f_D(y_1,\cdots,y_n|p)\pi(p)}{\int_{\Omega}f_D(y_1,\cdots,y_n|p)\pi(p)dp}\ &=\frac{p^k(1-p)^{(n-k)}}{\int_0^1p^k(1-p)^{(n-k)}dp}\ &=\frac{p^k(1-p)^{(n-k)}}{B(k+1,n-k+1)} \end{aligned} π(p∣y1​,⋯,yn​)​=∫Ω​fD​(y1​,⋯,yn​∣p)π(p)dpfD​(y1​,⋯,yn​∣p)π(p)​=∫01​pk(1−p)(n−k)dppk(1−p)(n−k)​=B(k+1,n−k+1)pk(1−p)(n−k)​​所以 p p p的期望为： E π [ p ∣ y 1 , ⋯   , y n ] = ∫ p π ( p ∣ y 1 , ⋯   , y n ) d p = ∫ 0 1 p ( k + 1 ) ( 1 − p ) ( n − k ) B ( k + 1 , n − k + 1 ) d p = B ( k + 2 , n − k + 1 ) B ( k + 1 , n − k + 1 ) = k + 1 n + 2 \begin{aligned} E_\pi[p|y_1,\cdots,y_n]&={\int}p\pi(p|y_1,\cdots,y_n)dp \ & = {\int_0^1}\frac{p^{(k+1)}(1-p)^{(n-k)}}{B(k+1,n-k+1)}dp \ & = \frac{B(k+2,n-k+1)}{B(k+1,n-k+1)}\ & = \frac{k+1}{n+2} \end{aligned} Eπ​[p∣y1​,⋯,yn​]​=∫pπ(p∣y1​,⋯,yn​)dp=∫01​B(k+1,n−k+1)p(k+1)(1−p)(n−k)​dp=B(k+1,n−k+1)B(k+2,n−k+1)​=n+2k+1​​
∴ P ( Y = 1 ) = k + 1 n + 2 \therefore \displaystyle P(Y=1)=\frac{k+1}{n+2} ∴P(Y=1)=n+2k+1​

习题1.2
  通过经验风险最小化推导极大似然估计。证明模型是条件概率分布，当损失函数是对数损失函数时，经验风险最小化等价于极大似然估计。

解答：

假设模型的条件概率分布是 P θ ( Y ∣ X ) P_{\theta}(Y|X) Pθ​(Y∣X)，现推导当损失函数是对数损失函数时，极大似然估计等价于经验风险最小化。
极大似然估计的似然函数为： L ( θ ) = ∏ D P θ ( Y ∣ X ) L(\theta)=\prod_D P_{\theta}(Y|X) L(θ)=D∏​Pθ​(Y∣X)两边取对数： ln ⁡ L ( θ ) = ∑ D ln ⁡ P θ ( Y ∣ X ) arg ⁡ max ⁡ θ ∑ D ln ⁡ P θ ( Y ∣ X ) = arg ⁡ min ⁡ θ ∑ D ( − ln ⁡ P θ ( Y ∣ X ) ) \ln L(\theta) = \sum_D \ln P_{\theta}(Y|X) \ \mathop{\arg \max}_{\theta} \sum_D \ln P_{\theta}(Y|X) = \mathop{\arg \min}_{\theta} \sum_D (- \ln P_{\theta}(Y|X)) lnL(θ)=D∑​lnPθ​(Y∣X)argmaxθ​D∑​lnPθ​(Y∣X)=argminθ​D∑​(−lnPθ​(Y∣X))
反之，经验风险最小化等价于极大似然估计，亦可通过经验风险最小化推导极大似然估计。

原文地址链接（注：本文章是用于个人学习）