她的眼睛同灯火重叠的那一瞬间,就像在夕阳的余晖里飞舞的妖艳而美丽的夜光虫。

# 预测

# 灰色预测

灰色预测的主要特点是模型使用的不是原始数据序列,而是生成的数据序列。(譬如累加生成的数据序列)

优点是不需要很多数据(4 个左右),就能解决历史数据少、序列完整性及可靠性低的问题。能利用微分方程来充分挖掘系统本质,精度高,可以将无规律的原始数据生成得到规律性较强的生成序列。

缺点是只适用于中短期预测,只适合指数增长的预测

# GM(1,1)

:北方某城市连续若干年道路交通平均噪声等级为:

x(0)=(71.1,72.4,72.4,72.1,71.4,72.0,71.6)x^{(0)}=(71.1,72.4,72.4,72.1,71.4,72.0,71.6)

  • 第一步,级比检验

    求级比λ(k)=x(0)(k1)x(0)(k),k=2,...,7\lambda(k)=\frac{x^{(0)}(k-1)}{x^{(0)}(k)},k=2,...,7

    λ=(λ(2),...,λ(7))=(0.982,1,1.0042,1.0098,0.9917,1.0056)\lambda=(\lambda(2),...,\lambda(7))=(0.982,1,1.0042,1.0098,0.9917,1.0056)

    发现所有的级比都落在(e2n+1,e2n+1).n=7(e^{-\frac{2}{n+1}},e^{\frac{2}{n+1}}).n=7 中,故认为可用x(0)x^{(0)} 作令人满意的GM(1,1)GM(1,1) 建模。

    注:若级比有没落在区间内的,可以做数据的平移转换,即x(0)+cx^{(0)}+c,可以收缩级比。

  • 第二步,GM(1,1)GM(1,1) 建模

    对原始数据作一次累加,有

    x(1)=(71.1,143.5,215.9,288,359.4,431.4,503)x^{(1)}=(71.1,143.5,215.9,288,359.4,431.4,503)

    构造数据矩阵BB 和数据向量YY

    B=(12(x(1)(1)+x(1)(2))112(x(1)(2)+x(1)(3))1......12(x(1)(6)+x(1)(7))1)Y=(x(0)(2)x(0)(3)...x(0)(7))B=\begin{pmatrix}-\frac{1}{2}(x^{(1)}(1)+x^{(1)}(2))&1\\ -\frac{1}{2}(x^{(1)}(2)+x^{(1)}(3))&1\\ ...&...\\ -\frac{1}{2}(x^{(1)}(6)+x^{(1)}(7))&1\\ \end{pmatrix}Y=\begin{pmatrix} x^{(0)}(2)\\ x^{(0)}(3)\\ ...\\ x^{(0)}(7) \end{pmatrix}

  • 第三步,计算

u^=(a^b^)=(BTB)1BTY=(0.002372.6573)\hat{u}=\begin{pmatrix}\hat{a}\\ \hat{b}\end{pmatrix}=(B^TB)^{-1}B^TY=\begin{pmatrix}0.0023\\ 72.6573 \end{pmatrix}

  • 第四步,建立模型:

    dx(1)dt+a^x(1)=b^x^(1)(k+1)=(x(0)(1)b^a^)ea^k+b^a^=30929e0.0023k+31000\frac{dx^{(1)}}{dt}+\hat{a}x^{(1)}=\hat{b}\\ \hat{x}^{(1)}(k+1)=(x^{(0)}(1)-\frac{\hat{b}}{\hat{a}})e^{-\hat{a}k}+\frac{\hat{b}}{\hat{a}}=-30929e^{-0.0023k}+31000

  • 第五步,令x^(1)(1)=x^(0)(1)=x(0)(1)=71.1\hat{x}^{(1)}(1)=\hat{x}^{(0)}(1)=x^{(0)}(1)=71.1,然后就可以取k=1,2,3,4,5,6k=1,2,3,4,5,6,用x^(1)(k+1)\hat{x}^{(1)}(k+1) 的表达式可以求出x^(1)\hat{x}^{(1)},再求差分还原x^(0)\hat{x}^{(0)}

  • 第六步,进行模型分析。主要指标:

    • 预测值和原始值的残差
    • 相对误差
    • 级比偏差,计算式为λ=110.5a1+0.5a×λ\lambda'=1-\frac{1-0.5a}{1+0.5a}\times\lambda。针对级比偏差,越小越好,小于 0.2 认为满意,小于 0.1 认为较高完成。其中aa 为发展系数,即a^\hat{a}

# GM(2,1)

对于非单调的、摆动发展序列,或有饱和的 S 形序列,可以考虑使用 GM(2,1)模型。

已知序列

x(0)=(41,49,61,78,96,104)x^{(0)}=(41,49,61,78,96,104)

  • 第一步,计算x(0)x^{(0)} 的 1-AGO 序列(一次累加序列)和 1-IAGO 序列(一次累减序列)为:

x(1)=(41,90,151,229,325,429)α(1)x(0)=(8,12,17,18,8)x^{(1)}=(41,90,151,229,325,429)\\ \alpha^{(1)}x^{(0)}=(8,12,17,18,8)

​ 计算x(1)x^{(1)} 的均值序列为(即相邻两个数求均值)

z(1)=(65.5,120.5,190,277,377)z^{(1)}=(65.5,120.5,190,277,377)

  • 第二步,构造矩阵

    B=(x(0)(2)z(1)(2)1x(0)(3)z(1)(3)1.........x(0)(6)z(1)(6)1)Y=(α(1)x(0)(2)α(1)x(0)(3)...α(1)x(0)(6))B=\begin{pmatrix} -x^{(0)}(2)&-z^{(1)}(2)&1\\ -x^{(0)}(3)&-z^{(1)}(3)&1\\ ...&...&...\\ -x^{(0)}(6)&-z^{(1)}(6)&1\\ \end{pmatrix} Y=\begin{pmatrix} \alpha^{(1)}x^{(0)}(2)\\ \alpha^{(1)}x^{(0)}(3)\\ ...\\ \alpha^{(1)}x^{(0)}(6)\\ \end{pmatrix}

  • 第三步,求解

    u^=(a^1a^2b^)=(BTB)1BTY\hat{u}=\begin{pmatrix}\hat{a}_1\\ \hat{a}_2\\ \hat{b} \end{pmatrix}=(B^TB)^{-1}B^TY

  • 第四步,建立 GM(2,1)白化模型

    d2x(1)dt2+a^1dx(1)dt+a^2x(1)=b^\frac{d^2x^{(1)}}{dt^2}+\hat{a}_1\frac{dx^{(1)}}{dt}+\hat{a}_2x^{(1)}=\hat{b}

    利用边界条件x(1)(1),x(1)(6)x^{(1)}(1),x^{(1)}(6) 解得x^(1)(t)\hat{x}^{(1)}(t)

  • 第五步,利用x^(1)(t)\hat{x}^{(1)}(t) 还原x^(0)\hat{x}^{(0)},并分析残差和相对误差。

# DGM(2,1)

B=(x(0)(2)1x(0)(3)1......x(0)(6)1)Y=(α(1)x(0)(2)α(1)x(0)(3)...α(1)x(0)(6))u^=(a^b^)=(BTB)1BTYB=\begin{pmatrix} -x^{(0)}(2)&1\\ -x^{(0)}(3)&1\\ ...&...\\ -x^{(0)}(6)&1\\ \end{pmatrix} Y=\begin{pmatrix} \alpha^{(1)}x^{(0)}(2)\\ \alpha^{(1)}x^{(0)}(3)\\ ...\\ \alpha^{(1)}x^{(0)}(6)\\ \end{pmatrix}\\ \hat{u}=\begin{pmatrix}\hat{a}\\ \hat{b} \end{pmatrix}=(B^TB)^{-1}B^TY

白化方程为

d2x(1)dt2+a^dx(1)dt=b^\frac{d^2x^{(1)}}{dt^2}+\hat{a}\frac{dx^{(1)}}{dt}=\hat{b}

边界条件是x(1)(1)x^{(1)}(1)dx(1)dt=x(1)\frac{dx^{(1)}}{dt}=x^{(1)} 求解方程。

# 灰色 Verhulst 模型

B=(z(1)(2)(z(1)(2))2z(1)(3)(z(1)(3))2......z(1)(3)(z(1)(3))2)Y=(x(0)(2)x(0)(3)...x(0)(6))u^=(a^b^)=(BTB)1BTYB=\begin{pmatrix} -z^{(1)}(2)&(z^{(1)}(2))^2\\ -z^{(1)}(3)&(z^{(1)}(3))^2\\ ...&...\\ -z^{(1)}(3)&(z^{(1)}(3))^2\\ \end{pmatrix} Y=\begin{pmatrix} x^{(0)}(2)\\ x^{(0)}(3)\\ ...\\ x^{(0)}(6)\\ \end{pmatrix}\\ \hat{u}=\begin{pmatrix}\hat{a}\\ \hat{b} \end{pmatrix}=(B^TB)^{-1}B^TY

白化方程为

dx(1)dt+a^x(1)=b^(x(1))2\frac{dx^{(1)}}{dt}+\hat{a}x^{(1)}=\hat{b}(x^{(1)})^2

边界条件x(1)(1)x^{(1)}(1)

# 微分方程预测

# 回归分析预测

# 马尔可夫预测

# 时间序列预测

将预测对象按照时间顺序排列起来,构成一个所谓的时间序列,根据时间序列过去的变化规律,推断今后变化的可能性及变化趋势、变化规律,这就是时间序列预测法。

时间序列预测本质也是一种回归模型,它一方面考虑到事物发展的持续性,运用过去的时间序列事件推测事物的发展趋势;另一方面又充分考虑到偶然因素产生的随机性。

优点是简单易行,便于掌握,能够充分运用原时间序列的数据,计算速度快,对模型参数有动态确定能力。

缺点是不能反映事物的内在联系,不能分析两个因素的相关关系,只适用于短期预测。

# 确定性时间序列分析法

一个时间序列往往是以下几类变化形式的组合:

  • 长期趋势变动TtT_t
  • 季节变动StS_t
  • 循环变动CtC_t
  • 不规则变动RtR_t,也分为突然变动和随机变动

常见的确定性时间序列模型有以下几种:

  • 加法模型

    yt=Tt+St+Ct+Rty_t=T_t+S_t+C_t+R_t

  • 乘法模型

    yt=TtStCtRty_t=T_t\cdot S_t\cdot C_t\cdot R_t

  • 混合模型

    yt=TtSt+Rtyt=St+TtCtRty_t=T_t\cdot S_t+R_t\\ y_t=S_t+T_t\cdot C_t\cdot R_t

其中,yty_t 是观测目标的观测记录,期望E(Rt)=0E(R_t)=0 方差Var(Rt)=σ2Var(R_t)=\sigma^2。可以使用一些经验方法来预测

# 移动平均法

移动平均法只适用于短期预测,而且是预测目标的发展趋势变化不大的情况。如果目标的发展趋势存在其他的变化,采用简单移动平均法就会产生较大的预测偏差和滞后。

取移动平均项数NN,有

yt^=1N(yt1+...+ytN)\hat{y_t}=\frac{1}{N}(y_{t-1}+...+y_{t-N})

这个方法非常简单,譬如数据:

y=(1,2,3,4,5)y=(1,2,3,4,5)

那么取N=3N=3,就可以预测:

y^=(1,2,3,1+2+33=2,2+3+23=73,3+2+7/33=229,...)\hat{y}=(1,2,3,\frac{1+2+3}{3}=2,\frac{2+3+2}{3}=\frac{7}{3},\frac{3+2+7/3}{3}=\frac{22}{9},...)

显然可以算一个预测标准误差:

S=t=N+1T(yt^yt)2TN=12[(42)2+(573)2]S=\sqrt{\frac{\sum_{t=N+1}^T(\hat{y_t}-y_t)^2}{T-N}}=\sqrt{\frac{1}{2}[(4-2)^2+(5-\frac{7}{3})^2]}

# 一次指数平滑法

设时间序列为y1,y2,...,yt,...y_1,y_2,...,y_t,...α\alpha 为加权系数,那么一次指数平滑法就是要:

yt+1^=αj=0(1α)jytj=αyt+(1α)yt^\hat{y_{t+1}}=\alpha\sum_{j=0}^\infty(1-\alpha)^jy_{t-j}\\ =\alpha y_{t}+(1-\alpha)\hat{y_t}

由于加权系数复合指数要求,而且具有平滑数据的功能,故称作指数平滑法。

关于加权系数的选择,可以发现,α\alpha 越大,新数据所占的比重就越大。反之,α\alpha 越小,历史数据所占比重越大。具体如何选择可采用以下一般规律:

  • 如果时间序列波动不大,则α\alpha 小一点,约在0.10.50.1-0.5 之间,以减少修正幅度,包含更多时间的信息
  • 如果时间序列具有迅速且明显的变化,则α\alpha 大一点,约在0.60.80.6-0.8 之间,使预测模型灵敏度高一点,迅速跟上数据变化。

# 二次指数平滑法

一次指数平滑法虽然克服了移动平均法的缺点,但当时间序列的变动出现直线趋势时,仍存在明显的滞后偏差。

因此必须加以修正,再作二次指数平滑,利用滞后偏差规律建立直线趋势模型,这就是二次指数平滑法。

计算公式为:

yt+1^(1)=αyt+(1α)yt^(1)yt+1^(2)=αyt^(1)+(1α)yt^(2)\hat{y_{t+1}}^{(1)}=\alpha y_t+(1-\alpha)\hat{y_{t}}^{(1)}\\ \hat{y_{t+1}}^{(2)}=\alpha \hat{y_{t}}^{(1)}+(1-\alpha)\hat{y_{t}}^{(2)}

其中,yt^(1)\hat{y_t}^{(1)} 是一次指数平滑值,和之前的计算公式相同。yt^(2)\hat{y_t}^{(2)} 是二次指数平滑值。

当时间序列开始呈现直线趋势时,可以用直线趋势模型进行预测:

yt+m^=at+btm{at=2yt^(1)yt^(2)bt=α1α(yt^(1)yt^(2))\hat{y_{t+m}}=a_t+b_tm\\ \begin{cases} a_t=2\hat{y_t}^{(1)}-\hat{y_t}^{(2)}\\ b_t=\frac{\alpha}{1-\alpha}(\hat{y_t}^{(1)}-\hat{y_t}^{(2)}) \end{cases}

然后就可以得到预测yt+1^=at+bt×1=at+bt\hat{y_{t+1}}=a_t+b_t\times 1=a_t+b_t

# 三次指数平滑法

当数据显示出二次曲线趋势时,而不是直线趋势,则可以考虑三次指数平滑法。

yt+1^(1)=αyt+(1α)yt^(1)yt+1^(2)=αyt^(1)+(1α)yt^(2)yt+1^(3)=αyt^(2)+(1α)yt^(3)\hat{y_{t+1}}^{(1)}=\alpha y_t+(1-\alpha)\hat{y_{t}}^{(1)}\\ \hat{y_{t+1}}^{(2)}=\alpha \hat{y_{t}}^{(1)}+(1-\alpha)\hat{y_{t}}^{(2)}\\ \hat{y_{t+1}}^{(3)}=\alpha \hat{y_{t}}^{(2)}+(1-\alpha)\hat{y_{t}}^{(3)}

三次指数平滑法预测模型为

yt+m^=at+btm+ctm2{at=3yt^(1)3yt^(2)+yt^(3)bt=α2(1α)2[(65α)yt^(1)2(54α)yt^(2)+(43α)yt^(3)]ct=α22(1α)2[yt^(1)2yt^(2)+yt^(3)]\hat{y_{t+m}}=a_t+b_tm+c_tm^2\\ \begin{cases} a_t=3\hat{y_t}^{(1)}-3\hat{y_t}^{(2)}+\hat{y_t}^{(3)}\\ b_t=\frac{\alpha}{2(1-\alpha)^2}[(6-5\alpha)\hat{y_t}^{(1)}-2(5-4\alpha)\hat{y_t}^{(2)}+(4-3\alpha)\hat{y_t}^{(3)}]\\ c_t=\frac{\alpha^2}{2(1-\alpha)^2}[\hat{y_t}^{(1)}-2\hat{y_t}^{(2)}+\hat{y_t}^{(3)}] \end{cases}

然后就可以预测yt+1^=at+bt+ct\hat{y_{t+1}}=a_t+b_t+c_t 了。

# 一阶差分指数平滑法

yt=ytyt1yt+1^=αyt+(1α)yt^yt+1^=yt+1^+yt\nabla y_t = y_t - y_{t - 1}\\ \nabla \hat{y_{t + 1}}=\alpha\nabla y_t + (1 - \alpha)\nabla\hat{y_t}\\ \hat{y_{t+1}}=\nabla\hat{y_{t+1}}+y_t

前面已经分析过,指数平滑值实际上是一种加权平均数,因此把序列中逐期增量的加权平均数(指数平滑值)加上当前值的实际数进行预测,比一次指数平滑法只用变量以往的取值的加权平均数作为下一期预测更合理。从而使预测值始终围绕实际值上下波动,从根本上解决了在有直线增长趋势的情况下,用一次指数平滑法所得出的结果始终落后于实际值的问题。

# 二阶差分指数平滑法

yt=ytyt12yt=ytyt12yt+1^=α2yt+(1α)2yt^yt+1^=2yt+1^+yt+yt\nabla y_t = y_t - y_{t - 1}\\ \nabla^2y_t=\nabla y_t-\nabla y_{t-1}\\ \nabla^2 \hat{y_{t + 1}}=\alpha\nabla^2 y_t + (1 - \alpha)\nabla^2\hat{y_t}\\ \hat{y_{t+1}}=\nabla^2\hat{y_{t+1}}+\nabla y_t+y_t

# 具有季节性特点的时间序列预测

这里的季节性,不仅是自然季节,也可以是某种产品的销售季节等。季节系数法的具体计算步骤为:

  • 收集mm 年的每年各年度或各月份(每年nn 个季度)的时间序列样本aija_{ij}ii 表示年份,jj 表示月份。

  • 计算每年所有季度的算术平均值:

    aˉ=1nmi=1mj=1naij\bar{a}=\frac{1}{nm}\sum_{i=1}^m\sum_{j=1}^na_{ij}

  • 计算同季度或同月份数据的算术平均值

    a.jˉ=1mi=1maij\bar{a_{.j}}=\frac{1}{m}\sum_{i=1}^ma_{ij}

  • 计算季度系数或月份系数

    bj=a.jˉ/aˉb_j=\bar{a_{.j}}/\bar{a}

  • 预测计算,当时间序列是按季度列出时,先求出预测年份(下一年)的年加权平均

    ym+1=i=1mwiyii=1mwiy_{m+1}=\frac{\sum_{i=1}^mw_iy_i}{\sum_{i=1}^mw_i}

    其中,yi=j=1naijy_i=\sum_{j=1}^n a_{ij} 为第ii 年的年合计数,wiw_i 为第ii 年的权数,按自然数列取值,即wi=iw_i=i

    再计算预测年份的季度平均值yˉm+1=ym+1/n\bar{y}_{m+1}=y_{m+1}/n。最后预测年份第jj 季度的预测值为

    ym+1,j=bjyˉm+1y_{m+1,j}=b_j\bar{y}_{m+1}

# 平稳时间序列模型

# AR§ 序列

{Xt,t=0,±1,±2,...}\{X_t,t=0,\pm 1,\pm 2,...\} 是零均值平稳序列,满足以下模型

Xt=φ1Xt1+φ2Xt2+...+φpXtp+εtX_t=\varphi_1 X_{t-1}+\varphi_2X_{t-2}+...+\varphi_pX_{t-p}+\varepsilon_t

其中,εt\varepsilon_t 是零均值、方差是σε2\sigma_{\varepsilon}^2 的平稳白噪声。称XtX_t 是阶数为pp自回归序列,简记为AR(p)AR(p) 序列。而

φ=[φ1,φ2,...,φp]T\varphi=[\varphi_1,\varphi_2,...,\varphi_p]^T

称为自回归参数向量,其分量φj\varphi_j 称为自回归系数。

引入后移算子B:BkXtXtkB:B^kX_t\equiv X_{t-k},记算子多项式

φ(B)=1φ1Bφ2B2...φpBp\varphi(B)=1-\varphi_1B-\varphi_2B^2-...-\varphi_pB^p

则原式也可改写为

φ(B)Xt=εt\varphi(B)X_t=\varepsilon_t

# ARMA (p,q) 序列

{Xt,t=0,±1,±2,...}\{X_t,t=0,\pm 1,\pm 2,...\} 是零均值平稳序列,满足以下模型

Xtφ1Xt1φ2Xt2...φpXtp=εtθ1εt1...θqεtqX_t-\varphi_1X_{t-1}-\varphi_2X_{t-2}-...-\varphi_pX_{t-p}=\varepsilon_t-\theta_1\varepsilon_{t-1}-...-\theta_q\varepsilon_{t-q}

其中,εt\varepsilon_t 是零均值、方差是σε2\sigma_\varepsilon^2 的平稳白噪声。则称XtX_t 是阶数为p,qp,q自回归移动平均序列,简记为ARMA(p,q)ARMA(p,q).

考虑后移算子BB 定义为Bkεt=εtkB^k\varepsilon_t=\varepsilon_{t-k}BkXt=XtkB^kX_t=X_{t-k},和多项式

φ(B)=1φ1Bφ2B2...φpBpθ(B)=1θ1Bθ2B2...θpBp\varphi(B)=1-\varphi_1B-\varphi_2B^2-...-\varphi_pB^p\\ \theta(B)=1-\theta_1B-\theta_2B^2-...-\theta_pB^p

于是原式改为φ(B)Xt=θ(B)εt\varphi(B)X_t=\theta(B)\varepsilon_t

对于一般的平稳序列{Xt,t=0,±1,±2,...}\{X_t,t=0,\pm 1,\pm 2,...\} 设其均值为E(Xt)=μE(X_t)=\mu 满足下列模型:

(Xtμ)φ1(Xt1μ)φ2(Xt2μ)...φpXtp=εtθ1εt1...θqεtq(X_t-\mu)-\varphi_1(X_{t-1}-\mu)-\varphi_2(X_{t-2}-\mu)-...-\varphi_pX_{t-p}=\varepsilon_t-\theta_1\varepsilon_{t-1}-...-\theta_q\varepsilon_{t-q}

于是该式可以重写为

φ(B)(Xtμ)=θ(B)εt\varphi(B)(X_t-\mu)=\theta(B)\varepsilon_t

关于算子多项式φ(B)\varphi(B)θ(B)\theta(B) 往往还要做以下假定

  • φ(B)\varphi(B)θ(B)\theta(B) 无公共因子,又φp0,θq0\varphi_p\neq 0,\theta_q\neq 0
  • φ(B)=0\varphi(B)=0 的根全在单位圆外,这一条件称为模型的平稳性条件。
  • θ(B)=0\theta(B)=0 的根全在单位圆外,这一条件称为模型的可逆性条件。

在实际问题的建模中,首先要估计阶数p,qp,q,然后对参数θ,φ\theta,\varphi 进行估计。

定阶和估计结束后,还要对模型进行检验,即检验εt\varepsilon_t 是否为平稳白噪声。若检验通过,则认为建模完成。

  1. ARMA 模型定阶的 AIC 准则:选p,qp,q 使得

    min  AIC=nlnσ^ε2+2(p+q+1)min\;AIC=nln\hat{\sigma}_\varepsilon^2+2(p+q+1)

    其中,nn 为样本容量,σ^ε2\hat{\sigma}_\varepsilon^2σϵ2\sigma_\epsilon^2 的一个估计,与ppqq 有关。若当p=p^,q=q^p=\hat{p},q=\hat{q} 时,AICAIC 达到最小值,那么就认为阶为p^,q^\hat{p},\hat{q}

    当序列还含有未知参数μ\mu 时,需最小化

    min  AIC=nlnσ^ε2+2(p+q+2)min\;AIC=nln\hat{\sigma}_\varepsilon^2+2(p+q+2)

    事实上和前式有相同的最小值点p^,q^\hat{p},\hat{q}

  2. ARMA 模型的参数估计:矩估计,逆函数估计,最小二乘估计,条件最小二乘估计,最大似然估计等方法。

  3. ARMA 模型的χ2\chi^2 检验

  4. ARMA 模型的预报

MATLAB 来全做了(

# ARIMA 序列与季节性序列

考虑一个序列的例子:

Xt1.5Xt1+0.5Xt2=εtX_t-1.5X_{t-1}+0.5X_{t-2}=\varepsilon_t

其中,φ(B)=11.5B+0.5B2\varphi(B)=1-1.5B+0.5B^2,它的两个根B1=1,B2=2B_1=1,B_2=2,并非都在单位圆外,故不满足平稳性条件。

若改写式子为

(XtXt1)0.5(Xt1Xt2)=εt(X_t-X_{t-1})-0.5(X_{t-1}-X_{t-2})=\varepsilon_t

Xt=XtXt1=Wt\nabla X_t = X_t - X_{t-1}=W_t,则有

Wt0.5Wt1=εtW_t-0.5W_{t-1}=\varepsilon_t

此时它是平稳的AR(1)AR(1) 序列,即把不稳定的序列XtX_t 转化成为了稳定的序列WtW_t

一般地,可以考虑qq 阶差分:

Xt=(1B)XtqXt=(1B)qXt\nabla X_t = (1-B)X_t\\ \nabla^qX_t = (1-B)^q X_t

若令WtqXtW_t\equiv\nabla^qX_t,然后WtW_t 是一个ARMA(p,q)ARMA(p,q) 序列的话,那么XtX_t 就被称作ARIMA(p,q,d)ARIMA(p,q,d) 序列。

# 小波分析预测

# 混沌序列预测

# 评价与决策

# 模糊综合评判

# 主成分分析法

# 层次分析法

# 数据包络分析法

# 秩和比综合评价法

# 投影寻踪综合评价法

# 方差、协方差分析

# 分类与判别

# 距离(系统)聚类

# 关联性聚类

# 层次聚类

# 密度聚类

# 贝叶斯分类器

以一个例子说朴素贝叶斯分类器,数据集:

纹理 (x1x_1) 色泽 (x2x_2) 敲声 (x3x_3) 类别 (yy)
1 清晰 清绿 清脆 好瓜
2 模糊 乌黑 浊响 坏瓜
3 模糊 清绿 浊响 坏瓜
4 清晰 乌黑 沉闷 好瓜
5 清晰 清绿 浊响 好瓜
6 模糊 乌黑 沉闷 坏瓜
7 清晰 乌黑 清脆 好瓜
8 模糊 清绿 沉闷 好瓜
9 清晰 乌黑 浊响 坏瓜
10 模糊 清绿 清脆 好瓜

就是有三个特征:纹理、色泽、敲声。然后需要判断是不是好瓜,

  • 其中,纹理有清晰、模糊两种水平
  • 其中,色泽有清绿、乌黑两种水平
  • 其中,敲声有浊响、沉闷、清脆三种水平

现在,我买了各瓜是 “清晰 + 清绿 + 沉闷” 的,那它是好瓜的概率是多少?


做法:

由贝叶斯公式:

P(\mbox{好瓜}|\mbox{纹理清晰、色泽清绿、敲声沉闷})=\frac{P(\mbox{纹理清晰、色泽清绿、敲声沉闷}|\mbox{好瓜})P(好瓜)}{P(\mbox{纹理清晰、色泽清绿、敲声沉闷})}

其中,根据全概率公式:

P(\mbox{纹理清晰、色泽清绿、敲声沉闷})=P(\mbox{纹理清晰、色泽清绿、敲声沉闷}|\mbox{好瓜})P(\mbox{好瓜})\\ +P(\mbox{纹理清晰、色泽清绿、敲声沉闷}|\mbox{坏瓜})P(\mbox{坏瓜})

然后贝叶斯分类有个非常重要的假设:不同特征间独立,即:

P(\mbox{纹理清晰、色泽清绿、敲声沉闷}|\mbox{好瓜})=P(\mbox{纹理清晰}|\mbox{好瓜})P(\mbox{色泽清绿}|\mbox{好瓜})P(\mbox{敲声沉闷}|\mbox{好瓜})

记 “纹理清晰” 为aa,“色泽清绿” 为bb,敲声沉闷为cc,好瓜为xx,坏瓜为yy,则有:

P(xabc)=P(abcx)P(x)P(abcx)P(x)+P(abcy)P(y)=P(ax)P(bx)P(cx)P(x)P(ax)P(bx)P(cx)P(x)+P(ay)P(by)P(cy)P(y)P(x|abc)=\frac{P(abc|x)P(x)}{P(abc|x)P(x)+P(abc|y)P(y)}\\ =\frac{P(a|x)P(b|x)P(c|x)P(x)}{P(a|x)P(b|x)P(c|x)P(x)+P(a|y)P(b|y)P(c|y)P(y)}

然后我们开始计算其中需要的概率:

  1. 先验概率P(x)P(x),只看数据集最后一列,有

    P(x)=6/10=3/5P(x)=6/10=3/5

  2. 先验概率P(y)P(y)

    P(y)=4/10=2/5P(y)=4/10=2/5

  3. 条件概率P(ax)P(a|x):看数据集中,最后一列是好瓜时,纹理清晰的概率:

    P(ax)=4/6=2/3P(a|x)=4/6=2/3

  4. 条件概率P(bx)P(b|x)

    P(bx)=4/6=2/3P(b|x)=4/6=2/3

  5. 条件概率P(cx)P(c|x)

    P(cx)=2/6=1/3P(c|x)=2/6=1/3

  6. 条件概率P(ay),P(by),P(cy)P(a|y),P(b|y),P(c|y)

    P(ay)=1/4P(by)=1/4P(cy)=1/4P(a|y)=1/4\\ P(b|y)=1/4\\ P(c|y)=1/4

所以最终要求的值为:

P(xabc)=P(ax)P(bx)P(cx)P(x)P(ax)P(bx)P(cx)P(x)+P(ay)P(by)P(cy)P(y)=2/32/31/33/52/32/31/33/5+1/41/41/42/5=12/13512/135+2/320=128137P(x|abc)=\frac{P(a|x)P(b|x)P(c|x)P(x)}{P(a|x)P(b|x)P(c|x)P(x)+P(a|y)P(b|y)P(c|y)P(y)}\\ =\frac{2/3\cdot 2/3\cdot 1/3\cdot 3/5}{2/3\cdot 2/3\cdot 1/3\cdot 3/5+1/4\cdot1/4\cdot1/4\cdot 2/5}\\ =\frac{12/135}{12/135+2/320}=\frac{128}{137}

所以说,纹理清晰、色泽清绿、敲声沉闷的瓜为好瓜的概率为128137\frac{128}{137},认为极大概率是好瓜。


说明:

  1. 当特征不是离散的,是连续的怎么办?譬如计算P(P(根蒂长度| 好瓜)=?)=?

​ 这时,我们可以假设该特征是正态分布的。然后计算好瓜情况下,数据集中根蒂长度的均值标准差,记为μ,σ\mu,\sigma。然后把数据集样本的均值和标准差,作为正态分布的参数估计。故有:

P(ax)=12πσ2exp((aμ)22σ2)P(a|x)=\frac{1}{\sqrt{2\pi\sigma^2}}exp(\frac{-(a-\mu)^2}{2\sigma^2})

​ 其中,aa 表示要预测的瓜的根蒂长度,μ,σ\mu,\sigma 是根据数据集算出来的。

  1. P(ax)=0P(a|x)=0 时怎么办?即数据集中没出现 “纹理清晰 + 好瓜” 的情况,这时会对分类器造成很大影响。

    我们引入 Laplace 校准,思想非常简单就是对每个类别(好瓜、坏瓜)的计数都加一。譬如P(ax)=(4+1)/(6+2)=5/8P'(a|x)=(4+1)/(6+2)=5/8。分子加一就是因为要平滑地加一个数字,分母加二是因为纹理有两个选择:清晰、模糊。

    Laplace 平滑公式:

    P(ax)=I(ax)+λI(x)+SλP(a|x)=\frac{I(a|x)+\lambda}{I(x)+S\cdot\lambda}

    其中,I(ax)I(a|x) 表示类别为xx 情况下,特征的选择为aa 的个数;I(x)I(x) 表示类别为xx 的个数;SS 为特征有几个选择(譬如纹理特征有清晰、模糊两个选择);λ\lambda 为平滑系数,选 1 就好。

    这样做后,当数据集较大时,其实对分类器没有影响。

# 模糊识别

# 关联与因果

# 灰色关联分析方法

# Spearman 或 Kendall 等级相关分析

# Person 相关

# Copula 相关

# 典型相关性分析

# 标准化回归分析

# 生存分析

# 格兰杰因果检验

# 优化

# 规划

# 排队论

排队论也称随机服务系统理论,它研究的内容有三部分:

  • 形态问题:各种排队系统的概率规律性,如队长分布、等待时间分布、忙期分布等;
  • 最优化问题:分静态最优 (最优设计) 和动态最优 (最优运营);
  • 排队系统的统计推断:判断一个给定的排队系统符合于哪种模型。

一般的排队过程由三部分组成:

  • 输入过程:顾客到来时间的规律性
    • 顾客组成可能有限无限
    • 顾客到达可能逐个成批
    • 顾客到达可以是相互独立相关
    • 输入过程可以平稳 (与时间无关)非平稳
  • 排队规则:顾客按怎样的规则排队等待
    • 损失制 (消失制):顾客到达时所有服务台均被占用,顾客随机离去;
    • 等待制:顾客到达时所有服务台均被占用,顾客排队等待,直到接受完服务才离去;
    • 混合制:既有等待又有损失,有队列长度有限排队等待时间有限两种情况。
  • 服务过程
    • 先到先服务;
    • 后到先服务;
    • 随机服务;
    • 优先服务。

排队模型用六个符号表示:X/Y/Z/A/B/CX/Y/Z/A/B/C

  • XX 表示顾客到达流或顾客到达间隔时间分布;
  • YY 表示服务时间分布;
  • ZZ 表示服务台数目;
  • AA 表示系统容量限制;
  • BB 表示顾客源数目;
  • CC 表示服务规则。

如略去后三项,即指X/Y/Z///FCFSX/Y/Z/\infty/\infty/FCFS

当顾客到达时间符合如下条件时,可以看作泊松分布:

  • 将时间段无限分割成若干小的时间段,在每个接近于零的小时间段里,事件发生一次的概率与时间段的长度接近成正比;
  • 在这个极小时间段内,事件发生二次及以上的概率恒等于 0;
  • 事件在不同的时间段内发生与否相互独立。

# M/M/1 模型

顾客到达时间间隔服从参数为 λ 的指数分布,服务时间服从参数为 μ 的指数分布,服务台数量为 1。

主要指标:

  • 服务强度ρ=λμ\rho=\frac{\lambda}{\mu}
  • 平稳时,系统中顾客人数为 n 的概率pn=(1ρ)ρnp_n=(1-\rho)\rho^n
  • 平均顾客数Ls=n=1npn=λμλL_s=\sum_{n=1}^\infty np_n=\frac{\lambda}{\mu-\lambda}
  • 平均排队长度Lq=n=1(n1)pn=λ2μ(μλ)L_q=\sum_{n=1}^\infty(n-1)p_n=\frac{\lambda^2}{\mu(\mu-\lambda)}
  • 平均逗留时间 TP(T>t)=e(μλ)tP(T>t)=e^{-(\mu-\lambda)t},TT 的期望为Ws=1μλW_s=\frac{1}{\mu-\lambda}
  • 平均等待时间Wq=λμ(μλ)W_q=\frac{\lambda}{\mu(\mu-\lambda)}

某医院急诊室同时只能接待一个病人,诊治时间服从指数分布,每个病人平均需要 15 分钟。病人按照泊松分布到达,平均每小时到达 3 人。试对此排队系统进行分析。

确定参数:λ=3\lambda =3 人 / 小时,μ=6015=4\mu=\frac{60}{15}=4 人 / 小时

故服务强度为ρ=3/4=0.75\rho=3/4=0.75,稳态概率为P0=(1ρ)ρ0=0.25P_0=(1-\rho)\rho^0=0.25,这就是急诊室空闲的概率,也是病人不必等待立刻能就诊的概率。

故病人需要等待的概率为0.250.25,也就是急诊室繁忙的概率。

急诊室内外病人总数的期望为Ls=343=3L_s=\frac{3}{4-3}=3

急诊室外等待病人数的期望为Lq=324(43)=2.25L_q=\frac{3^2}{4(4-3)}=2.25

病人在急诊室内外逗留时间期望为Ws=143=1W_s=\frac{1}{4-3}=1 小时

病人平均等待时间为Wq=34(43)=0.75W_q=\frac{3}{4(4-3)}=0.75 小时,也等于总逗留时间减去接诊时间 15 分钟。

# M/M/S 模型

顾客到达时间间隔服从参数为 λ 的指数分布,服务时间服从参数为 μ 的指数分布(服务率都相同),服务台数量为 S。

主要指标

  • 总服务强度ρ=λSμ\rho^*=\frac{\lambda}{S\mu}

  • 平稳时,系统中顾客人数为 n 的概率

P0=[k=0S11k!(λμ)k+1S!11ρ(λμS)]1Pn={1n!(λμ)nP00<nS1S!SnS(λμ)nP0SnP_0=[\sum_{k=0}^{S-1}\frac{1}{k!}(\frac{\lambda}{\mu})^k+\frac{1}{S!}\frac{1}{1-\rho^*}(\frac{\lambda}{\mu}^S)]^{-1}\\ P_n=\begin{cases} \frac{1}{n!}(\frac{\lambda}{\mu})^nP_0&0<n\leq S\\ \frac{1}{S!S^{n-S}}(\frac{\lambda}{\mu})^nP_0&S\leq n \end{cases}

  • 平均排队长度Lq=(Sρ)SρS!(1ρ)2P0L_q=\frac{(S\rho^*)^S\rho^*}{S!(1-\rho^*)^2}P_0
  • 平均顾客数Ls=Lq+SρL_s=L_q+S\rho^*
  • 平均等待时间:W_q=\frac{L_q}
  • 平均逗留时间 T:W_s=W_q+\frac{1}
  • 排队人数过多的概率P(Nk)=(Sρ)kk!(1ρ)P0P(N\geq k)=\frac{(S\rho^*)^k}{k!(1-\rho^*)}P_0

:某医院急诊室能同时接待两个个病人,诊治时间服从指数分布,每个病人平均需要 15 分钟。病人按照泊松分布到达,平均每小时到达 3 人。试对此排队系统进行分析。

类似前例,有λ=3\lambda=3 人 / 小时,μ=4\mu=4 人 / 小时,S=2S=2ρ=3/8\rho^*=3/8

急诊室空闲的概率为P0=[1+0.75+0.7522!(10.375)]1=511P_0=[1+0.75+\frac{0.75^2}{2!(1-0.375)}]^{-1}=\frac{5}{11}

病人需要等待,即P(N2)P(N\geq 2) 的概率为P(N2)=(0.75)22!(10.375)5/110.2P(N\geq 2)=\frac{(0.75)^2}{2!(1-0.375)}*5/11\approx 0.2

急诊室外等待病人数的期望为Lq=0.275/110.12L_q=0.27*5/11\approx 0.12

急诊室内外总人数为Ls=Lq+Sρ=0.87L_s=L_q+S\rho^*=0.87

病人在急诊室内外逗留时间期望为Ws=Lsλ=0.29W_s=\frac{L_s}{\lambda}=0.29 小时

病人平均等待时间为Wq=Lqλ=0.04W_q=\frac{L_q}{\lambda}=0.04 小时=0.29=0.29 小时 - 15 分钟(接诊时间)


下面可以证明一下,如果一个事件单位时间内发生次数服从泊松分布,那么这个它发生间隔的时间服从指数分布。

建模为,事件AA 在任意长度时间tt 内发生的次数N(t)N(t) 服从参数为λt\lambda t 的泊松分布,其中λ\lambda 为单位时间内发生次数:

P{N(t)=k}=(λt)kk!eλtP\{N(t)=k\}=\frac{(\lambda t)^k}{k!}e^{-\lambda t}\\

然后关键来了,事件{T>t}\{T>t\} 和事件N(t)=0N(t)=0 等价,其中TT 为事件发生的时间间隔。为什么呢?就是两次事件间隔时间大于tt,就等价于,时间长度为tt 中没发生事件。所以

F(T)=1P{T>t}=1eλtf(T)=F(T)=λeλtF(T)=1-P\{T>t\}=1-e^{-\lambda t}\\ f(T)=F'(T)=\lambda e^{-\lambda t}

TT 是服从指数分布的。


补充说明:泊松事件流的要求

  • 增量独立性,对于tn>tn1>...>t1t_n>t_{n-1}>...>t_1N(tn)N(tn1),...,N(t2)N(t1)N(t_n)-N(t_{n-1}),...,N(t_2)-N(t_1) 之间互相独立。

  • 增量平稳性,即一段时间tt 内发生的次数只与tt 有关,而与tt 的起始时刻和终止时刻无关。即[0,t)[0,t)[a,a+t)[a,a+t) 内发生次数为kk 的概率相同。

  • 普遍性,即当tt 充分小时,P(N(t)=1)=λt+o(t),P(N(t)>1)=o(t)P(N(t)=1)=\lambda t+o(t),P(N(t)>1)=o(t),即在充分小的时间段内事件发生一次的概率相对于时间为常数:

    limt0P(N(t)=1)t=λlimt0P(N(t)>1)t=0limt0P(N(t)=0)t=1λ\lim_{t\rightarrow 0}\frac{P(N(t)=1)}{t}=\lambda\\ \lim_{t\rightarrow 0}\frac{P(N(t)>1)}{t}=0\\ \lim_{t\rightarrow 0}\frac{P(N(t)=0)}{t}=1-\lambda\\

# 模糊约束

# 灰色规划

# 网格优化

# 柔性约束