数学建模社区-数学中国

标题: python时间序列预测(ARIMA模型案例代码) [打印本页]

作者: zhangtt123    时间: 2020-5-23 14:43
标题: python时间序列预测(ARIMA模型案例代码)
1、模型识别0 ]# h/ Q# T0 k
01 主要的模型
% Y2 \& W% v4 mAR(P)模型(Autoregressive Model)" T/ L  ^$ i4 u. b- J; T; h

/ X- `7 X& E6 P& Q    自回归模型描述的是当前值与历史值之间的关系/ m8 u' |& F, N; b6 |
1 ~$ N" P4 i2 h4 R
MA(q)模型(Moving Average Model), _# h. c  v& W
     移动平均模型描述的是自回归部分的误差累计
4 f8 X5 Y3 s( \' ?1 o; J1 x
0 h/ j) I8 P) T8 l7 } ARIMA模型(Autoregressive Integrated Moving Average Model)$ H3 o. a9 k, A5 x) B$ z
+ t( K0 L4 O1 w3 Y; X, k
    所谓ARIMA模型,是指将非平稳时间序列转化为平稳时间序列,然后将因变量仅对它的之后值以及随机误差项的现值和滞后值进行回归所建立的模型
5 h4 w% Y8 ~0 L3 X, G$ {% z( I
" C& T" ^: B6 Q# a+ J! k    Xt=自回归AR+移动平均MA模型
$ ~1 t+ M- R* e0 k1 b4 m 4.png * V# ^0 z0 \/ V
02 截尾和拖尾
& C2 _# t/ S3 K: |' B(1)p阶自回归模型 AR(P) + ]" S) }+ p& h
AR(p)模型的偏自相关函数PACF在p阶之后应为零,称其具有截尾性; : T6 L  B- c3 E6 ]! f6 [+ o, N
AR(p)模型的自相关函数ACF不能在某一步之后为零(截尾),而是按指数衰减(或成正弦波形式),称其具有拖尾性。
' @& h! m+ {  Q( B* j$ K- U# t% H5 J
(2)q阶移动平均模型 MA(q) / A6 Z6 i& b& x2 G
MA(q)模型的自相关函数ACF在q阶之后应为零,称其具有截尾性;
3 i1 P. m) b$ r; A6 TMA(q)模型的偏自相关函数PACF不能在某一步之后为零(截尾),而是按指数衰减(或成正弦波形式),称其具有拖尾性。- k! a* E. ?% N& r, B! u) r

- W% ~, C7 u: r8 _6 L& F# n7 d03 如何判断拖尾和截尾6 L( N. P7 a; b0 c. T6 {: l
, D2 P3 [9 A3 R9 v$ a
(1)如果样本自相关系数(或偏自相关系数)在最初的d阶明显大于2倍标准差范围,而后几乎95%的样本自相关(偏自相关)系数都落在2倍标准差范围以内,而且由非零自相关(偏自相关)系数衰减为小值波动的过程非常突然,这时,通常视为自相关(偏自相关)系数截尾。
$ I) {8 o% T0 i* V2 [' W1 P4 J3 f
(2)如果有超过5%的样本相关系数落在2倍标准差范围以外,或者是由显著非零的相关函数衰减为小值波动的过程比较缓慢或者非常连续,这时,通常视为相关系数不截尾。) m9 ]9 w' T  q3 ^/ }+ L5 L
3 p% k. m# O4 U( i$ {
2、时间序列算法公式5 Y% _! T7 g0 H  g
重要的几种为:AR、MA、ARMA、ARIMA模型,具体公式见下图:
! ?' {& w+ D3 M- n* W0 B7 B/ ~& `  p) l& y

/ M1 X+ C  s( F/ T$ M# n1 ^; y0 W$ z: K5 N2 j3 z
3、详细步骤
6 m' x2 p0 o& n01 平稳性检验(adf检验)1 R3 @0 P; D7 i: d9 \
+ q+ D3 L, n4 o' U# \; Z1 U
#  a  时序图检验
7 W" Z7 [4 g. K1 f* r. x5 `
# O8 K, L0 i' q+ v根据平稳时间序列的均值和方差都为常数的性质,平稳序列的时序图显示该序列值始终在一个常熟附近随机波动,而且波动的范围有界;如果有明显的趋势性或者周期性,那他通常不是平稳序列
/ K. Z# B3 d' d- c9 ]- F
7 q* Z0 C+ ^) G# X: X#  b 自相关图检验。( g" m0 n5 R( G# o8 ?0 m( ?0 O

. X5 z9 @$ B4 d% D平稳序列具有短期相关性民政性质对平稳序列而言通常只有近期的序列值对现时值的影响比较明显,间隔越远的过去只对现时值得影响越小。随着延迟期数K的增加,平稳序列的自相关系数Pk(延迟K期)会比较快的衰减趋向于零,并在零附近随机波动,而非平稳序列的自相关系数衰减的速度比较慢,这就是利用自相关图进行平稳性检验的标准$ ^8 N- M% {! E  e0 d

; L8 l. u1 D$ A
: P+ X* t6 v6 c" [, B" _#c 单位根检验  
3 l) w$ x! h' t& e" a# K6 X5 H; e8 ?* h/ V! M0 r
单位根检验是指检验序列中是否存在单位根,如果存在单位根就是非平稳时间序列了7 @6 E! J+ ^" ]. |  Q
4 B7 ~8 F: K/ ~$ X( O" _1 [  S

) E$ |& W& X4 ]6 ?# ADF单位根检验
" _5 D; J" Q- x$ u+ R6 ]( `& O5 ?
/ ]7 m" }) x) P& c4 J其中第二中的ADF检验中,如果p值显著大于0.05,统计量化大于三个或者两个水平值,差距越大,越不平稳。
' \1 o; W+ S( P9 e4 u0 `5 U' {8 K! u$ O. [1 X* Q& E+ Y- P
若统计量显著小于三个置信度且p值接近0,为平稳序列' X! f2 U3 |/ O

) W& `4 F( g: O. v6 v其他情况,可能是数据量不够的原因没有展现趋势' k, H* z0 l  m$ b9 J8 \0 i
' n9 x3 M5 |' x+ ?
02 对数据进行差分构造平稳序列. D: v! C% ]! V
差分运算:
( Y) {% v2 W+ k6 C7 o  A$ A
6 s: P' \1 b7 _p阶差分, u5 s  ~' F# y2 {4 g# d
# I4 W6 b" I( h/ V& m; I7 G
相聚一期的两个序列值之间的减法运算称为I阶差分运算1 v  E2 X2 A' L* O8 Y1 P( x; q

( T/ a& P+ h$ N$ O+ i! Mk步差分0 w9 A, Q! O1 h

  g& J# `8 p% r相聚k期的两个序列值之间的减法运算称为k步chafenyunsuan. g9 Q& r1 Q3 }; k8 N- I2 Q* T
1 F' a( P1 M. H& q
#差分后的结果# n/ K+ I# e& W/ q/ A
D_data = data.diff(k).dropna()
" M  n( m# O" v* p) r+ g* v4 e: n& }3 [+ o
k 相距k时间
& {" ?! U" @: D2 i) I# U- b* t( c/ _# r& R' J
一般在一阶差分后就是比较平稳的序列
( E2 v% E: Z) k1 z' f6 L+ Y( W# e! j# j1 `: l6 r
03 平稳性监测! t. {  E% |. T
04 白噪声检验
; O) S+ x" @1 U" Z5 u& D! B3 X#白噪声检验
1 {8 r/ p2 q! J) Y* K: P7 xfrom statsmodels.stats.diagnostic import acorr_ljungbox
8 \3 p$ H! X+ j! E
$ S( e$ O5 m4 j% Z7 ]#返回统计量和p值
- \" r. y( j8 q4 M) ~# ~6 A. n0 W' t, h2 u  _" Z
print(u'差分序列的白噪声检验结果为:', acorr_ljungbox(D_data, lags=1))  # 分别为stat值(统计量)和P值
( b" O0 J5 A, g" v! M! S2 _: R+ b* ]6 g5 s7 B) S( i
# P值小于0.05,所以一阶差分后的序列为平稳非白噪声序列。
. `8 P, Y7 J9 Q! X" _. O! ]$ y) T
3 M. @' m4 Z; B05 定阶) U% q, p$ P& i3 E* v5 y( j
第一种方法人为识别5 |/ L3 y9 w' A8 y) x; h5 A
3 t3 R. b, u' D( B
一阶差分后自相关图显示1阶截尾,
- [) v5 T; @: I/ H! g
; R1 S3 Q) I- U3 b偏自相关显示托尾性,所以建立ARIMA(0,1,1)
  \& [+ D# W0 P" Z% }* L% Z7 B" m9 v! D) A% H! L4 a# R
第二种:相对最优模型识别5 u0 q% I8 |% C( ?

$ U3 \# W- R3 [计算ARMA(p,q)。当p,q均小于所组合BIC信息量,取其中BIC信息量达到最小的模型阶数' z5 v# q: n! E; @5 f

6 E$ k4 v3 d% i确定pq值# ^& C; M$ w4 _. T

2 l6 `1 y# t) W( bfrom statsmodels.tsa.arima_model import ARIMA8 t* x2 w$ |, l' j3 T. O  @
#定阶
% L7 @9 n9 F) x( M+ V6 d
' I4 F: g; z( I. Y2 M#一般阶数不超过length/10
1 O" p2 _8 w# Z" j# G1 [! R2 [" r1 {: I* D. X' d
pmax = int(len(D_data)/10)
7 \$ Y% _, U4 Q( {) B
7 i4 |7 y5 G3 {' ]1 u#一般阶数不超过length/109 G' ~. f3 p! k- O' C0 T! @

# b* y7 |% L; q5 ^$ tqmax = int(len(D_data)/10) " W: f: \; q0 O, ^; ]
/ A1 V9 u# l8 b& I8 X, d

4 r5 c+ F% V8 B9 v6 {#bic矩阵
# P2 H% ^; s+ A% R  T7 G9 j" n( d. s! \; c6 n  S% G" w
bic_matrix = [] $ n8 A9 _' T6 o
for p in range(pmax+1):
% o" Y  _/ [5 s3 Z9 q. M: H/ h  tmp = []# x9 a" V& B( e" v
  for q in range(qmax+1):
7 m6 Y' `2 `' c4 a7 v" t! f4 q#存在部分报错,所以用try来跳过报错。' A! S" o, s2 U0 j0 Z/ _5 N
    try:
" \2 a. S$ P2 P$ _2 b; G      tmp.append(ARIMA(data, (p,1,q)).fit().bic)
/ r' A1 Z4 }3 C6 P    except:4 Z  g5 z- L, i4 y; C+ F
      tmp.append(None)  d/ D/ J2 [  Q  w6 p
  bic_matrix.append(tmp)* e) E' H5 z7 O+ I. N' T4 ]; @

; F& D* E! }+ P6 o* _5 {# `0 C#从中可以找出最小值
' s& [& M, y4 n( ~; r% d. l( d- l3 G- h# L8 z. p
bic_matrix = pandas.DataFrame(bic_matrix) & x; a% K0 w3 x) J% w6 j9 @7 O

0 Q2 M- b) t; E4 o- n) H' J#先用stack展平,然后用idxmin找出最小值位置。2 p' R8 f- m7 k

: {' ^9 \* I' G: K, x/ M1 X1 Wp,q = bic_matrix.stack().idxmin()
4 \+ ?, d6 Z& o" L
& D6 U( @2 O- w0 H$ Y( r+ a$ U& pprint(u'BIC最小的p值和q值为:%s、%s' %(p,q)); s' s2 U- }5 F; R0 o6 N; y0 ]$ a3 V

2 y9 \: u0 Y; |) b. {# y# 取BIC信息量达到最小的模型阶数,结果p为0,q为1,定阶完成。6 _' P! Z- j; u

* `6 F1 c" r, o8 x" a06 模型预测2 u. y1 c* A& R) i. ?8 T1 E
7 w6 r1 G: ^) A& W
#建立ARIMA(0, 1, 1)模型# g  k' b* I7 O
4 X' d, C  B0 N3 z! o9 P
model = ARIMA(data, (p,1,q)).fit() 3 e6 ^) @# j3 l6 b

1 @0 z, J% M* U; q' m( Z# Q# q#给出一份模型报告
1 R6 f2 u& A3 S1 H9 x4 f+ e/ A; C0 x, b6 j. p! k, p1 a
model.summary2() 9 R8 f+ H0 a2 T. O+ L: P
5 ^+ F( b: [! E4 ?$ P( @9 }
#作为期5天的预测,返回预测结果、标准误差、置信区间。( M) z1 s* B( K' X. ]/ {( {

! `8 I7 i  Z7 v: t; \model.forecast(5)" D& t  Z4 ?: R' F' \; Z% B: A( N
9 u$ r7 [+ A  x
' o, {8 |) H( i8 c9 E( o+ T
4、案例代码
6 a, G" e$ o. Y& f, v* A1 f: X9 I' W# e. n( m6 ?" r( J/ [) @" k
import pandas
; V  R2 i$ e/ O
2 U( m+ N9 i- H6 k# 读取数据,指定日期为索引列1 K7 A0 Y" b. N+ x. e

! a9 @! M. ^" h' V+ pdata = pandas.read_csv(5 T. h3 Y7 y0 y1 T
    'D:\\DATA\\pycase\\number2\\9.3\\Data.csv' ,. e$ `0 N: c+ U! @' L1 d- a( T
    index_col='日期'
6 r4 a7 x" i, y: B2 U3 x)/ y& B! q( s' z8 Y

2 [2 A8 R" E$ l6 X, [7 ]( O- m$ z# 绘图过程中' r2 B" s8 t# p% y# r

" i/ ^/ I( X5 r( j0 Cimport  matplotlib.pyplot as plt
0 E- o2 g( s% ]) H0 j  w. `3 G8 x( G3 X  O- w$ v2 h- r
# 用来正常显示中文标签
6 [9 p+ I. }4 |9 i# H" s
# z% k- C/ @% ~2 \; M' aplt.rcParams['font.sans-serif']=['SimHei']" `7 H  {. T- z% D, M
' e. p* t! M% k+ g( Z4 d; B/ g
# 用来正常显示负号; T' Y& J/ g: @8 l
3 c0 g2 }1 S/ N& C- ]1 i+ @
plt.rcParams['axes.unicode_minus'] = False ! [6 }0 T% s" [5 O- B4 Y
9 E: @& s" Q" C' x
# 查看趋势图: K  i) L$ D) K* m1 A3 r
data.plot() #有增长趋势,不平稳/ d& q+ M7 W  q0 T

0 W8 I- q6 ]4 C& x( h& f( d7 U# U1 D* y+ m
# 附加:查看自相关系数合片自相关系数(查分之后),可以用于平稳性的检测,也可用于定阶系数预估( D$ r% G6 Z( b
1 r  G$ r7 w" t  y
#自相关图()2 h+ X  b4 o, i4 O) X' {4 F

, L+ ?) {5 n+ z2 j$ Lfrom statsmodels.graphics.tsaplots import plot_acf
4 u& M# |  N+ E- P# j
2 Q* O' z; g1 I9 s4 q/ d. lplot_acf(data).show() #自相关图既不是拖尾也不是截尾。以上的图的自相关是一个三角对称的形式,这种趋势是单调趋势的典型图形,说明这个序列不是平稳序列4 P# H& e' }& D; @3 i+ s6 S

9 t) D' S5 Z: U0 m; A# M0 x# e: S! s

' c5 y  V# y% v" r8 S1 m& I* W- G: v
: }9 l# \  C4 G( x3 A8 G
# 1 平稳性检测( p" R* k$ Z  Z

4 D- z3 f* g0 K' O* B& Y+ Mfrom statsmodels.tsa.stattools import adfuller as ADF
, j0 B1 g  g' k4 M# }& g6 m4 L- ?" E( w# @9 r5 W: l: l
5 s: L' o2 x) q: \" d& W
def tagADF(t):: k+ t) ]1 S3 U" P; q
    result = pandas.DataFrame(index=[
! D; ^& [; B0 y* w9 x# U1 s            "Test Statistic Value", "p-value", "Lags Used", - e$ c. |! `# g: @# v) `
            "Number of Observations Used", * k$ p, j' z# \( U$ }* A( }
            "Critical Value(1%)", "Critical Value(5%)", "Critical Value(10%)"/ l7 ?7 V' `( N% }' N4 {
        ], columns=['销量']) ?- A( j* m: a; K
    );
% @4 ], B! h& r' B    result['销量']['Test Statistic Value'] = t[0]( q/ Y+ J+ \( O+ y
    result['销量']['p-value'] = t[1]
, ?" h' S' ^9 v) F3 |( R  V    result['销量']['Lags Used'] = t[2]5 t/ N( W( i' r4 [" ?
    result['销量']['Number of Observations Used'] = t[3]
8 |3 w0 w$ _8 |% z: s    result['销量']['Critical Value(1%)'] = t[4]['1%']
" I! x* x  p$ k+ ~* U    result['销量']['Critical Value(5%)'] = t[4]['5%']
; I) h7 E' \  E/ [$ F3 ?    result['销量']['Critical Value(10%)'] = t[4]['10%']& z6 @# x; I; s6 e
    return result;
0 ^" D# r, D; J% y9 e: M! c! [$ w$ c7 o1 ~: D

1 G" L; w: K1 d; m' P* W4 Iprint('原始序列的ADF检验结果为:',tagADF(ADF(data[u'销量'])))  # 添加标签后展现: n" U% n' u# X6 A# ^+ T

  o" Q1 P8 `% a' A/ f# 平稳判断:得到统计量大于三个置信度(1%,5%,10%)临界统计值,p值显著大于0.05,该序列为非平稳序列。/ T: S0 D! D; d/ ]
# 备注:得到的统计量显著小于3个置信度(1%,5%,10%)的临界统计值时,为平稳 此时p值接近于0 此处不为0,尝试增加数据量,原数据太少( y5 {# I6 k, c) d
, [, |; t! o: r# Q" i- w, l) F
# 2 进行数据差分,一般一阶差分就可以
# [# n& p. K  [/ ?- e  d* k  @1 R9 m6 z7 W
D_data = data.diff(1).dropna()" l" D5 [  \$ f$ E) i: c9 \
D_data.columns = [u'销量差分']1 A- H% _6 ~) E. }
$ x' d% y; ^' B
#差分图趋势查看" F3 H: y9 T+ X0 b2 b0 C' v) n
( y5 ]$ n) N0 \7 X6 `* l0 R! ~
D_data.plot()
5 s+ N  B+ v- e5 p2 m, iplt.show()0 V, a2 e/ c+ g7 {/ T2 x# }

" m6 {5 p1 D3 m  D6 s, o$ z# 附加:查看自相关系数合片自相关系数(查分之后),可以用于平稳性的检测,也可用于定阶系数预估
4 H* j. i! z1 F7 m8 n
% v& D1 K+ f4 L# r& z7 w#自相关图! V5 M& N% Z- f* N0 c
* @" f$ |9 ?  I9 M
plot_acf(D_data).show()/ j- i/ d. D. V% `7 o, v2 W

$ g: ]$ c5 V7 C+ @plt.show()* ~  K2 Q8 @0 d$ ?

) }3 I) d& F2 |% V6 \#偏自相关图
/ r4 d. N/ E) m: \! ?! n' _. c! i
from statsmodels.graphics.tsaplots import plot_pacf  N: y  n; r, [

" B: h1 A3 R4 Hplot_pacf(D_data).show()
, d2 P( G# r: _: p/ e
! Z- q( d* y% V: c/ u9 r5 F# 3 平稳性检测
/ W8 G& B& e% W6 w' ~4 e. ~7 |6 V- u
; L. w+ _& {: P  i9 [0 w5 L- aprint(u'差分序列的ADF检验结果为:', tagADF(ADF(D_data[u'销量差分'])))
. `3 M* ~  F7 C1 R4 o/ p3 y+ `& D  y( d
# 解释:Test Statistic Value值小于两个水平值,p值显著小于0.05,一阶差分后序列为平稳序列。  ?0 P6 n, h) n7 x

! E: F1 e; Q. K# 4 白噪声检验; }* @+ }. Y( u3 w; A0 H
from statsmodels.stats.diagnostic import acorr_ljungbox
6 V5 ]" n: [: W# }( m9 f
( [; F6 B$ n- s7 m* h#返回统计量和p值/ O" L, ]( n# m+ m7 @; L
4 O8 M* ~, Z, W$ x( R
print(u'差分序列的白噪声检验结果为:', acorr_ljungbox(D_data, lags=1))  # 分别为stat值(统计量)和P值& J# ~) B, k8 o/ Q

5 H) P9 z4 u- w- J; e) i# P值小于0.05,所以一阶差分后的序列为平稳非白噪声序列。
9 ~# J4 _5 |& s5 |  B& W6 i5 h$ O% n: M* K) l# x5 m

. _* ^, U: c7 P  H9 E0 |# 5 p,q定阶/ w" k  [. A# M) K8 L
% k- k7 f( s+ ?# l" f; R# J
from statsmodels.tsa.arima_model import ARIMA" @; {, F4 ?4 h6 B# @6 R" x& w
! I/ Z+ K5 Y! |& R
#一般阶数不超过length/10
8 e1 i+ _* C5 E' A/ Y7 Z4 f/ h" x# S' }! x  ]
pmax = int(len(D_data)/10)
' A( Q5 i* P6 U0 B! Y( M) U
1 v& E8 ^" A9 M* a. e4 O: x( Z: Q8 ?
#一般阶数不超过length/10
, R* r) {) K0 o  a8 X, }' n% Z
- Z, l6 a1 Q% rqmax = int(len(D_data)/10)
* @. k( i0 ^  A& w* o2 w4 k
3 ~6 m/ k* {% G" c# @- H#bic矩阵
5 D1 X/ d/ K3 g/ b3 X( T% v; \1 F, y& h9 b% I0 f/ i
bic_matrix = []
* f5 g, b9 P" k- g8 b1 Cfor p in range(pmax+1):3 e  ~' n# C. w
  tmp = []
6 ]% k! G" F$ K) S! S, V  ~  for q in range(qmax+1):
  R: E2 d  N! E; P. V#存在部分报错,所以用try来跳过报错。& Q- x! w( f1 o7 R8 o) P- `8 X" q
    try: 5 Q! a# D& B9 o7 [
      tmp.append(ARIMA(data, (p,1,q)).fit().bic)6 J" d& k' G6 g2 S1 i% I8 R
    except:
" B8 ]; d( |$ H8 h. p5 |# j( D1 O. }      tmp.append(None): u5 t1 X2 H2 Q
  bic_matrix.append(tmp); [/ I( r7 C8 l: h* i4 O! H/ ^* m2 k

/ e5 J" x! {4 H0 v4 F2 m#从中可以找出最小值# O4 B- z0 T9 I3 ~8 T7 W

6 D  l. i' F6 T! ?! k5 `+ P; bbic_matrix = pandas.DataFrame(bic_matrix)
  a8 n8 E5 F, P0 I" [: m, H; ~/ a$ E( Y! i) ^1 w( Q
#先用stack展平,然后用idxmin找出最小值位置。
, |9 V* \6 F/ ^8 p& d
& e; U6 K' c. W. K' Op,q = bic_matrix.stack().idxmin()
9 E* b  b7 [! h4 [% P  y+ U3 X/ o) N1 @1 C6 M( z( }. r1 w9 R

4 _( x3 e3 M4 k8 X4 M; z* e6 h3 N# W8 Y2 A3 t, K4 }+ o. M
print(u'BIC最小的p值和q值为:%s、%s' %(p,q))
# b8 T0 F5 R3 j9 k$ {# 取BIC信息量达到最小的模型阶数,结果p为0,q为1,定阶完成。
0 @1 _0 w' v# P: a7 |7 y8 R: f- U# c! O3 A% S
# 6 建立模型和预测
0 V5 ]/ P! k1 Y0 D7 p: C) a# d/ W
8 O- E5 [& M$ i3 I0 k- g! x  x7 ymodel = ARIMA(data, (p,1,q)).fit()
6 p! U: h8 ~+ D$ i2 Z6 I" r* U
, V+ x$ u6 P, f  |. Z#给出一份模型报告% \9 h" H8 y, R, i4 U2 C+ U. M+ a
! `5 O* P9 U4 q; K8 z5 p* X- e
model.summary2()
( k) B4 \% x: g: W! ~, r7 n8 w& q% w
# q/ V$ z2 F1 a0 V/ u#作为期5天的预测,返回预测结果、标准误差、置信区间。
, S  z" _& B; Z: b
$ y* x9 v1 B% [# Vmodel.forecast(5)
) f) H3 B  N/ E8 W( a' @# D& x2 `8 f1 |7 P
( `) Q4 D) x4 a
————————————————/ r7 B! i% z' l5 J, B0 ~
版权声明:本文为CSDN博主「UP Lee」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。- Y. ~% r& A+ [2 K- M
原文链接:https://blog.csdn.net/qq_36327687/article/details/85696152
8 |: Q& W& C( I' C
  q2 ?. E( \6 W  J% o* P  B0 l1 k0 S3 u. f( }

作者: 柠檬草lll    时间: 2020-5-25 19:06
发表回复谢谢分享) ~" R$ P; ?' X; s  G. q; x; w





欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5