多项式函数拟合sin函数（最小二乘法求解参数及其正则化）

杨利霞

多项式函数拟合sin函数（最小二乘法求解参数及其正则化）
( X" R! \" w* z0 e* I, t0 {  y
. h) U  H# a/ A: _& Y8 P% M0 M2 v  h1．统计学习是关于计算机基于数据构建概率统计模型并运用模型对数据进行分析与预测的一门学科。统计学习包括监督学习、非监督学习、半监督学习和强化学习。
* S# X/ ^3 |1 i8 q2．统计学习方法三要素——模型、策略、算法，对理解统计学习方法起到提纲挈领的作用。
0 S, _) J+ S* G3．本书主要讨论监督学习，监督学习可以概括如下：从给定有限的训练数据出发， 假设数据是独立同分布的，而且假设模型属于某个假设空间，应用某一评价准则，从假设空间中选取一个最优的模型，使它对已给训练数据及未知测试数据在给定评价标准意义下有最准确的预测。
4．统计学习中，进行模型选择或者说提高学习的泛化能力是一个重要问题。如果只考虑减少训练误差，就可能产生过拟合现象。模型选择的方法有正则化与交叉验证。学习方法泛化能力的分析是统计学习理论研究的重要课题。
5．分类问题、标注问题和回归问题都是监督学习的重要问题。本书中介绍的统计学习方法包括感知机、K近邻法、朴素贝叶斯法、决策树、逻辑斯谛回归与最大熵模型、支持向量机、提升方法、EM 算法、隐马尔可夫模型和条件随机场。这些方法是主要的分类、标注以及回归方法。它们又可以归类为生成方法与判别方法。
* r- G$ E0 F* H8 u: [; v" u

; }: y9 f# [6 ]


import numpy as np
- c3 G- n, r& p% g+ j! A/ Fimport matplotlib.pyplot as plt
( g( y0 \5 D7 e$ v5 ~8 ifrom scipy.optimize import leastsq
$ n8 W. h" m5 o9 r+ e
/ Z6 w# Z# ?- A2 P$ G7 M
% j4 k; ]4 f: D9 `6 ]0 @% ^# 我们要拟合的目标函数
def real_func(x):
6 Q* u3 D4 l1 i; e% q. M+ ~7 q6 @  b    return np.sin(2*np.pi*x)

0 A+ x; ]$ x5 w2 F
# 我们自己定义的多项式函数
5 I! f) b& P3 E* edef fit_func(p, x):
    f = np.poly1d(p)  # np.poly1d([2,3,5,7])返回的是函数，2x3 + 3x2 + 5x + 7
- O( W4 w; d+ W# @    ret = f(x)
    return ret
+ E* H1 Y) T& Q8 Z
5 l: N& n4 M6 l7 x2 \0 O# O1 i! T
# 计算残差
def residuals_func(p, x, y):
    ret = fit_func(p, x) - y
3 F  d0 C& M3 b2 U    return ret
, p. l9 P4 ?& U1 O+ {7 u# O) G

def fitting(M=0):
    """
        M    为 多项式的次数
    """
    # 随机初始化多项式参数
    p_init = np.random.rand(M + 1)  # 返回M+1个随机数作为多项式的参数
    # 最小二乘法:具体函数的用法参见我的博客：残差函数，残差函数中参数一，其他的参数
    p_lsq = leastsq(residuals_func, p_init, args=(x, y))
    # 求解出来的是多项式当中的参数，就是最小二乘法中拟合曲线的系数
" U6 s( n2 W9 Z( _+ ^* C    # print('Fitting Parameters:', p_lsq[0])
: e4 Y. E7 [! s( }4 f$ w* z    return p_lsq[0]
; v! `' W3 `) A# H; N3 F" n6 U
: m8 |/ N- c6 g; i
6 F) b! b$ n) w) _, _! }1 k# 书中10个点,对y加上了正态分布的残差
4 A0 j. K7 b0 V# G3 @1 [3 ox = np.linspace(0, 1, 10)
y_old = real_func(x)
y = [np.random.normal(0, 0.1) + yi for yi in y_old]
* o& ~2 }# s5 L7 ^
8 a7 x! m& t' n, Q3 [
6 ?; C% f& x+ Ux_real = np.linspace(0, 1, 1000)
y_real = real_func(x_real)
  N& _0 ?; K9 R* c

; b2 T! A# B0 W) lplt.plot(x_real, y_real, label="real")
plt.plot(x, y, 'bo', label='point')
#  fiitting函数中args=（x， y）是条用的是上面定义的10个点的全局变量x,y
3 Y6 \0 o, S2 r6 Vplt.plot(x_real, fit_func(fitting(9), x_real), label="fitted curve")
% W" L. z1 J) z$ G4 wplt.legend()
plt.show()
$ T# y6 R/ Q2 k5 N) S
- J# B, D3 t# e! f( NM=0
% i; ]; A4 z! \1 W. y' ^
$ @& A+ V$ n8 _9 s* T9 T2 ~
6 ?7 {9 j' f( B: i" lM=1

% m3 y% g" q  J6 B0 eM=3
4 V* v) T$ K# ]2 H( q3 A# [
6 _: A" }* r3 ?: i
3 O/ [3 k1 w8 m& h
2 j$ `0 ^% O8 a# Q5 `8 `1 fM=9
6 }3 C& E0 `& ~8 L1 w1 Q. E9 }
. H  u+ V2 M- }

W是参数，就是最小二乘法求得到的系数
lambda是regularization，是自定义的系数。
import numpy as np
# p& k2 ]$ Y7 f' f8 ?- ~3 [import matplotlib.pyplot as plt
# ~$ J% M3 k/ S3 \5 m- Y+ kfrom scipy.optimize import leastsq


# 我们要拟合的目标函数
def real_func(x):
/ e  w. U# P4 W( h. s- l    return np.sin(2*np.pi*x)
  }) Z) \; w5 J! r

# 我们自己定义的多项式函数
def fit_func(p, x):
    f = np.poly1d(p)  # np.poly1d([2,3,5,7])返回的是函数，2x3 + 3x2 + 5x + 7
( }) P- s# h; _3 ]    ret = f(x)
0 N% B0 }* E9 q2 I4 \! s    return ret

) {. r% X8 G2 M' ]& f$ r+ D
# 计算残差
9 `; H/ J4 _5 ]0 pdef residuals_func(p, x, y):
    ret = fit_func(p, x) - y
    return ret
" ]8 f& w" Q6 Y
' l! {  q5 u  m0 d1 O
# B7 d# E7 J8 }2 S# C, ^. k/ g# 返回残差和正则项
0 O! F% \. L, q" a1 n2 z! W! ddef residuals_func_regularization(p, x, y):
    ret = fit_func(p, x) - y
3 ]& n, E6 I. g$ G4 B9 {) i. U6 Q    ret = np.append(ret,
                    np.sqrt(0.5 * regularization * np.square(p)))  # L2范数作为正则化项
    return ret
, Y) \1 R* `6 j( v0 C* m

def fitting(M=0):
    """
        M    为 多项式的次数
    """
/ K- ?& x8 B! M- D4 V9 j    # 随机初始化多项式参数
7 z! ]* T( W/ X6 e4 {    p_init = np.random.rand(M + 1)  # 返回M+1个随机数作为多项式的参数
    # 最小二乘法:具体函数的用法参见我的博客：残差函数，残差函数中参数一，其他的参数
    p_lsq = leastsq(residuals_func, p_init, args=(x, y))
    # 求解出来的是多项式当中的参数，就是最小二乘法中拟合曲线的系数
    # print('Fitting Parameters:', p_lsq[0])
$ i5 }2 h" f7 w  v* q1 h5 U    return p_lsq[0]

" N) y! s$ m- v  Q
- {* C! p+ E/ F: z! }6 u: e+ X6 _3 s# 书中10个点,对y加上了正态分布的残差
: a8 u0 U5 i8 f- {4 {; Qx = np.linspace(0, 1, 10)
y_old = real_func(x)
y = [np.random.normal(0, 0.1) + yi for yi in y_old]


0 K& f+ w/ O8 B/ d' ix_real = np.linspace(0, 1, 1000)
  L1 J: L; G7 O6 Qy_real = real_func(x_real)

$ v) a" b& \) V& c, B
# # 画出10个散点，sin图像，和拟合的曲线
# plt.plot(x_real, y_real, label="real")
* k6 M. S( T/ v# plt.plot(x, y, 'bo', label='point')
# plt.plot(x_real, fit_func(fitting(9), x_real), label="fitted curve")
# plt.legend()
+ B( c! W- W6 [+ n# plt.show()

+ y7 Q" V* o2 S8 _4 F/ k
# 画出添加正则项的曲线
$ i: f0 D1 J3 B9 j+ oregularization = 0.0001
) ^6 I" b, Y) T  q" g1 Gp_init = np.random.rand(9 + 1)
p_lsq_regularization = leastsq(
    residuals_func_regularization, p_init, args=(x, y))
4 O  p# A, `3 E( v  _9 C+ b; X1 h
5 w5 {, ~5 G: }3 R6 O
# 画出原sin图像，不加正则项的图像，加上正则项的图像，10个点的散点图
# 不加正则项和加上正则项都是9次方，10个系数
plt.plot(x_real, real_func(x_real), label='real')
plt.plot(x_real, fit_func(fitting(9), x_real), label='fitted curve')
plt.plot(
* W8 k) f* q2 S) {$ B    x_real,
- M  x4 M* {1 X7 Z    fit_func(p_lsq_regularization[0], x_real),
2 \+ D! e* _. ~2 l    label='regularization')
plt.plot(x, y, 'bo', label='noise')
plt.legend()
" E: y) G3 P+ ~2 ?8 |- `) |, z% m2 Hplt.show()
1 R# ]7 F% m" y+ j, J( E
, F0 E2 |, R0 m) `( Y


' m' @+ W2 w% V( i0 `2 c! O3 ~