请选择 进入手机版 | 继续访问电脑版

QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 914|回复: 0

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

[复制链接]
字体大小: 正常 放大
杨利霞        

5250

主题

81

听众

16万

积分

  • TA的每日心情
    开心
    2021-8-11 17:59
  • 签到天数: 17 天

    [LV.4]偶尔看看III

    网络挑战赛参赛者

    网络挑战赛参赛者

    自我介绍
    本人女,毕业于内蒙古科技大学,担任文职专业,毕业专业英语。

    群组2018美赛大象算法课程

    群组2018美赛护航培训课程

    群组2019年 数学中国站长建

    群组2019年数据分析师课程

    群组2018年大象老师国赛优

    发表于 2020-4-25 16:12 |显示全部楼层
    |招呼Ta 关注Ta
    多项式函数拟合sin函数(最小二乘法求解参数及其正则化)- q, a  [& g3 t* B  c

    ; O; e0 q/ ]) ]5 d* E* H1.统计学习是关于计算机基于数据构建概率统计模型并运用模型对数据进行分析与预测的一门学科。统计学习包括监督学习、非监督学习、半监督学习和强化学习。
    , X7 ?9 q8 e' O; C! g, W2.统计学习方法三要素——模型、策略、算法,对理解统计学习方法起到提纲挈领的作用。6 q! P9 |5 L6 d- S
    3.本书主要讨论监督学习,监督学习可以概括如下:从给定有限的训练数据出发, 假设数据是独立同分布的,而且假设模型属于某个假设空间,应用某一评价准则,从假设空间中选取一个最优的模型,使它对已给训练数据及未知测试数据在给定评价标准意义下有最准确的预测。
    * p' a% i, A1 V7 f4.统计学习中,进行模型选择或者说提高学习的泛化能力是一个重要问题。如果只考虑减少训练误差,就可能产生过拟合现象。模型选择的方法有正则化与交叉验证。学习方法泛化能力的分析是统计学习理论研究的重要课题。8 @3 K8 f. c/ c( E2 a
    5.分类问题、标注问题和回归问题都是监督学习的重要问题。本书中介绍的统计学习方法包括感知机、K近邻法、朴素贝叶斯法、决策树、逻辑斯谛回归与最大熵模型、支持向量机、提升方法、EM 算法、隐马尔可夫模型和条件随机场。这些方法是主要的分类、标注以及回归方法。它们又可以归类为生成方法与判别方法。. E3 e4 V9 M' W' F

    * ^) n2 z: X; T, l: q. _1 l  l& F4 c! z% q" {6 \+ [
    1.png

    8 @2 |$ v+ D: t) X4 q, B
    * @( o- Q7 K& {# P9 \ 2.png
    # C; Q. e; I* {" z; ^
    import numpy as np
    & ], j# d# o' f4 wimport matplotlib.pyplot as plt
    ; c3 q3 n0 [0 A2 @5 zfrom scipy.optimize import leastsq. g4 Z# E# X- K  _2 n6 N
    8 T! L! d# G# ?& v8 X4 X
    % a( e8 G9 t3 s& K4 M, V
    # 我们要拟合的目标函数3 U& b; o* D- K8 c" Y* ?* b' Z$ G) q6 I
    def real_func(x):
      x1 U/ u, h. D) `1 j0 ?    return np.sin(2*np.pi*x)
    / w- L+ l8 U- q8 ~6 Z
    . X( e1 N3 H! y6 T  @, j! ?! Q3 c: j% ]
    # 我们自己定义的多项式函数
    ' s  U. U6 G' p4 D- Kdef fit_func(p, x):3 h! `+ B, A; q- j: J. `) Q
        f = np.poly1d(p)  # np.poly1d([2,3,5,7])返回的是函数,2x3 + 3x2 + 5x + 7# `' t6 c1 }5 v& b* R  N/ d
        ret = f(x)* a" U: X8 v: P( w
        return ret
    . O5 E4 [. S/ t5 u$ r7 v% A
    3 V3 |+ ]8 ?/ L7 J8 I$ E
    5 [& ?9 t6 n+ H6 M! t& t6 R% l" f# 计算残差7 J. m9 j  D( b/ r/ d
    def residuals_func(p, x, y):! Z6 U% J5 c. p* P! C* S0 E
        ret = fit_func(p, x) - y
    1 O3 [5 a- u( G7 |0 l+ h( A    return ret- u- b6 Q! S  N, o& U+ f

    5 H  d3 g, U- R% \& I9 n# u/ V$ W+ i! T# i0 ]
    def fitting(M=0):# b$ V2 i4 a' T$ y. b# k' M
        """
    " I+ E8 ^1 b1 `) `8 s        M    为 多项式的次数4 b) D- g( B: b$ I4 I( y( u
        """
    ' @2 e" M/ E$ D' z2 z    # 随机初始化多项式参数
    ; n4 q" W$ g, Z" Z+ c    p_init = np.random.rand(M + 1)  # 返回M+1个随机数作为多项式的参数
    , N9 _: `# W4 A7 \+ i    # 最小二乘法:具体函数的用法参见我的博客:残差函数,残差函数中参数一,其他的参数
    8 [* C  r6 W- Z! \# o  I* ]    p_lsq = leastsq(residuals_func, p_init, args=(x, y))' R" L, S$ b* ~
        # 求解出来的是多项式当中的参数,就是最小二乘法中拟合曲线的系数
    + E. a1 K4 f1 U& ?3 L0 ~. `    # print('Fitting Parameters:', p_lsq[0])
    8 e* i8 E6 Q* Y$ L    return p_lsq[0]
    0 Z) [3 [! @1 D! g( w& i& F& q+ R9 e( @7 A" V! i7 r9 X
    6 q! P1 E8 |6 N' r
    # 书中10个点,对y加上了正态分布的残差
    * ~; _% c6 [) R" @, N/ Jx = np.linspace(0, 1, 10)
    0 u. k1 S( V$ o& g1 m) @y_old = real_func(x)! |6 z' j: R0 m5 v! V4 F- F; w
    y = [np.random.normal(0, 0.1) + yi for yi in y_old]
    " Y5 e( h9 ~* Q3 L. B& M# g
    ) D2 F8 ?* ~/ Z1 C4 ]  c) e
    & h5 W: g" x5 Z# gx_real = np.linspace(0, 1, 1000)
    2 `' D" _: M4 L9 r; a. q) w6 Oy_real = real_func(x_real)
    / }  i( Y0 H0 T0 J/ P0 r/ ^6 `4 z! }! @8 N; ^& A9 f. M" F1 ?
    % ]% a+ i; C) G
    plt.plot(x_real, y_real, label="real")
    : n1 u/ |$ ?0 W) F( J% Q( rplt.plot(x, y, 'bo', label='point') 5 ^& r$ u! j1 {4 N+ E
    #  fiitting函数中args=(x, y)是条用的是上面定义的10个点的全局变量x,y
    1 O& n$ G8 j8 v9 C6 i; j+ Bplt.plot(x_real, fit_func(fitting(9), x_real), label="fitted curve")( I+ S" _- Z& T# G0 z7 Z- j
    plt.legend()
    - p3 `2 a, Y2 p! O, F; [plt.show()
    8 m' n, }- n  a! R8 X7 a# h1 G( Q8 m2 Z2 g, q$ v- f
    M=0  I6 Q1 |7 D4 d4 O! _
    - f3 t/ v/ b" I3 I' X" a
    3.png
    % d  G3 E! N7 h4 ?9 v6 F6 Q4 B/ i! q
    M=16 r5 J5 T' ]4 N9 d. a
    4.png

    * i& B. ?& s9 U+ e4 A' g* IM=38 B$ X! H8 M* X& e( g6 ^( p; ~6 U

    ( S& w+ i5 M4 ^ 5.png

    : U# `. J- [: }
    7 ?1 q1 G$ q$ f; }  Z  KM=99 B0 y5 K5 ~% C% q4 K
    6.png
    " A0 v* U" D2 [9 T# y5 C, s
    7.png

    + \/ n) H9 [" j7 X) q
    : D4 H% j$ A: D6 s% nW是参数,就是最小二乘法求得到的系数/ S9 l# x0 Y' j% }& W7 T; @  W& e
    lambda是regularization,是自定义的系数。/ a( d; l: O/ \+ n/ q
    import numpy as np/ I( V% ^6 k' Z
    import matplotlib.pyplot as plt
    * F0 p7 w' @* |: c2 ?7 afrom scipy.optimize import leastsq7 o# ]. T6 I( c0 h4 m

    $ I$ ~$ \' z3 Z8 Z2 \1 N
    2 r# w/ m2 G/ f" _1 c6 m; u* |# 我们要拟合的目标函数
    4 F9 }9 k/ l+ I$ ^def real_func(x):) [4 W/ T4 c  B% N7 M
        return np.sin(2*np.pi*x)4 m& l/ I  H6 b' H9 D
    + i" u+ t/ m# L$ R

    4 W, \4 n5 z( I$ P; X! W# }# 我们自己定义的多项式函数
    + X2 l# k5 H6 s$ A* w. Kdef fit_func(p, x):0 M! U1 L8 k* a( f
        f = np.poly1d(p)  # np.poly1d([2,3,5,7])返回的是函数,2x3 + 3x2 + 5x + 7
    0 _: n  W+ M4 M0 _- J- Z    ret = f(x)
    $ c  @" m! z( l1 i    return ret
    + i8 K2 z. H9 E* B. B" y6 d9 \# t
    9 y. i8 M" Z, F$ P# L, e
    : h/ R  Q3 U: B. v7 N, O' A/ C' T1 V# 计算残差. z+ i( j- d9 T5 r) [
    def residuals_func(p, x, y):. Y+ b2 N5 O) q9 I- Z8 {6 ~
        ret = fit_func(p, x) - y  T0 U* B% y# E# [" p" I
        return ret3 @, }, d2 {+ G1 o& j
    7 F1 f( S0 M  ]& g, Y) U
    & l# [7 s: X9 ?4 a! M0 R% Q
    # 返回残差和正则项  S9 ^( U$ Q( [/ p: l7 [. @
    def residuals_func_regularization(p, x, y):
    ' n& d8 n" X1 D5 C- e: f3 ~    ret = fit_func(p, x) - y
    ) M* }3 `: }" l1 W1 ?6 e/ f    ret = np.append(ret,8 m0 H) k8 L9 f! q
                        np.sqrt(0.5 * regularization * np.square(p)))  # L2范数作为正则化项* m4 B/ Z4 Q7 c0 ?7 B/ `1 N, }
        return ret2 S4 q6 y) c& m/ J! u# O0 p
    : [/ G0 }: Z' V: a& `. h
    0 b1 A0 u7 P0 D0 ~9 P
    def fitting(M=0):
    & y: w& l* k1 t/ V! U    """
    ! v' k5 l2 |4 e  v' v        M    为 多项式的次数6 F9 q& {; R( i
        """% {: E$ c: n3 `5 Q9 v1 i# O
        # 随机初始化多项式参数
    ( q$ z! n/ X4 F9 T4 B# _    p_init = np.random.rand(M + 1)  # 返回M+1个随机数作为多项式的参数% }, `$ k- t: l; h( l
        # 最小二乘法:具体函数的用法参见我的博客:残差函数,残差函数中参数一,其他的参数0 d0 j. T0 [8 o6 y4 l; G
        p_lsq = leastsq(residuals_func, p_init, args=(x, y))
    / T* B* [7 y4 m# j* N) Z7 J    # 求解出来的是多项式当中的参数,就是最小二乘法中拟合曲线的系数8 g# \+ z0 ]# w) ^. J( m
        # print('Fitting Parameters:', p_lsq[0]): S" @/ K5 p: w: H$ }* B+ w7 e
        return p_lsq[0]: [& e  Z3 \' }6 m( g1 p

    ; j8 W) }; l, D8 C8 Z8 H7 |
    ! S& ~8 x5 @, D. S8 k% S/ M# 书中10个点,对y加上了正态分布的残差
    2 D1 _1 A2 x; w. Lx = np.linspace(0, 1, 10)& N# X7 d4 }+ m! L3 U
    y_old = real_func(x)
    4 Y; C# b8 |& G$ Ty = [np.random.normal(0, 0.1) + yi for yi in y_old]- X6 u. C* C# S2 K1 @

    2 _9 d4 p9 a( L7 g1 c1 \
    # h* n2 x' t6 E& `x_real = np.linspace(0, 1, 1000)
    ; N0 j$ k. F2 uy_real = real_func(x_real)
    6 m/ ^5 I  g0 t( P5 h+ _$ f( n1 c+ [+ m* H" p7 W! O/ @
    , x; N2 k4 B3 X
    # # 画出10个散点,sin图像,和拟合的曲线
    , W0 b7 U- n; {  \* }' n0 K# plt.plot(x_real, y_real, label="real"): i4 P, c( b  q# M
    # plt.plot(x, y, 'bo', label='point')6 T, `% y! |6 c' I
    # plt.plot(x_real, fit_func(fitting(9), x_real), label="fitted curve")8 o; i5 ?/ [' n4 ~
    # plt.legend()
    ' Q/ B$ @, Q0 A# plt.show()
    * V0 I4 }$ D! R9 U) X% y( L+ w" y  ]6 T9 l% e% S

    # `. @: x' _) r! U* D2 g) H# 画出添加正则项的曲线1 n& [' D5 r5 Y3 N; F
    regularization = 0.0001
    6 n* b( {& {6 y; z7 H' N) A5 C, zp_init = np.random.rand(9 + 1)  y" P6 N( C  Z# e4 s8 A
    p_lsq_regularization = leastsq(/ I* i4 {7 _& ~( k/ ^
        residuals_func_regularization, p_init, args=(x, y))8 a0 p( i2 J9 L+ @2 I0 g
    ; R: [. }# @9 p- F8 ^4 f2 A

    5 v3 @  V2 B& y: H# 画出原sin图像,不加正则项的图像,加上正则项的图像,10个点的散点图
    9 ]) m$ v) V, [/ `. {! r# 不加正则项和加上正则项都是9次方,10个系数$ O  b2 L6 I- l" j& [2 k
    plt.plot(x_real, real_func(x_real), label='real')
    % f/ e# \. y5 }) Gplt.plot(x_real, fit_func(fitting(9), x_real), label='fitted curve')/ l0 }# B" W. r& K
    plt.plot(0 i; U5 ?" S5 o( G* U
        x_real,
      r1 m( b3 N; o+ f, T    fit_func(p_lsq_regularization[0], x_real),/ G: u+ N4 n# `; ]) Z$ t. V: k( r
        label='regularization')5 h' q- e! e( \' }) B
    plt.plot(x, y, 'bo', label='noise')# Q2 g7 @  A! J# A) }
    plt.legend()
    0 _9 I  }- P  [* ?. D# `* }plt.show()# m: d* ^0 w3 }2 L3 L2 p% a
    % q0 F- @/ q; i. I
    8.png

    1 m) g' a$ A- J3 S- m" y9 k6 ]1 U3 Z+ y

    & P7 Z  ~  [! e7 n( {/ ?; Y5 m$ h7 i4 [9 [: ]# C; @% f
    zan
    您需要登录后才可以回帖 登录 | 注册地址

    qq
    收缩
    • 电话咨询

    • 04714969085
    fastpost

    关于我们| 联系我们| 诚征英才| 对外合作| 产品服务| QQ

    手机版|Archiver| |繁體中文 手机客户端  

    蒙公网安备 15010502000194号

    Powered by Discuz! X2.5   © 2001-2013 数学建模网-数学中国 ( 蒙ICP备14002410号-3 蒙BBS备-0002号 )     论坛法律顾问:王兆丰

    GMT+8, 2024-3-28 21:22 , Processed in 0.353717 second(s), 54 queries .

    回顶部