统计预测与决策实验报告

总成绩：

江苏师范大学科文学院

实 验 报 告

课程： 统计预测与决策 班级： 姓名： 学号： 教师：

实验一：多元线性回归模型

实验目的与要求：

熟练掌握建立多元线性回归模型的方法。 实验内容：

问题：国际旅游外汇收入是国民经济发展的重要组成部分，影响一个国家或地区旅游收入的因素包括自然、文化、社会、经济、交通等多方面的因素，本例研究第三产业对旅游外汇收入的影响。《中国统计年鉴》把第三产业划分为12个组成部分，分别为x1农林牧渔服务业，x2地质勘查水利管理业，x3交通运输仓储和邮电通信业，x4批发零售贸易和餐饮业，x5金融保险业，x6房地产业，x7社会服务业，x8卫生体育和社会福利业，x9教育文化艺术和广播，x10科学研究和综合艺术，x11党政机关，x12其他行业。选取1998年我国31个省、市、自治区的数据（见实验一数据.xls）。自变量单位为亿元人民币，以国际旅游外汇收入为因变量y（百万美元）。试建立线性回归模型。

（要求用MATLAB的stepwise函数解决问题。取进0.05,出0.1。）

解：实验步骤：

1.File—Import Data—将data 重命名为A（建立数组A）—Finish 2.File-Save workspace as -shiyan1_1_1.mat ；

3.File—New—M-File—输入代码—Debug—Save As—文件名（Untitled）—保存； 程序代码：

load shiyan1_1_1.mat [n,m]=size(A);

X=A(:,1:m-1);Y=A(:,m);

stepwise(X,Y,[1 2 3 4 5 6 7 8 9 10 11 12],0.05,0.1)

运行结果:

由图可以看出所得的线性回归方程为：

y184.7634.32121x320.202x817.3653x911.6183x1013.0049x11

实验二：时间序列分解法建模

实验目的与要求：

熟练掌握运用时间序列分解法建模。 实验内容：

问题：当将时间序列分解成长期趋势、季节变动、周期变动和不规则变动四个因素后，可以认为时间序列Yt是这四个因素的函数，即：

Ytf(Tt,St,Ct,It)

时间序列分解方法有很多，相对而言，乘法模型（YtTtStCtIt）应用得比较广泛。试就文件（实验二数据.xls）提供的数据，将实际销售额（Y）分解为T、S、C和I的乘积。（只需给出T、S和C即可） 解：实验步骤：

1、统计—时间序列—时间序列图—简单—确定—序列（C3：销售额Y）—数据增量（1）—确定；

2、统计—时间序列—移动平均—变量（C3：销售额Y）—移动平均长度（4）—移动平均居中—时间（指数）—存储（移动平均）—确定；

3、计算—计算器—将结果存储在变量中（C5）—表达式（C3/C5）；

4、数据—拆分列—拆分的数据在（C6）—使用的下标在（C1）—在最后使用的一列之后—确定；

5、数据—转置列—转置以下列（C10-C13）—确定； 6、计算—行统计量—均值—输入变量（C15-C26）—将结果存储在（C28）—确定—得到C28：同季平均；

7、计算—计算器—将结果存储在变量中（C29）—表达式（C28 /SUM(C28) * 4）—得到C29：季节指数（%）；

8、在C2列依次输入1,2…48；

9、统计—时间序列—趋势分析—变量（C3）—模型类型（线性）—时间（标记：t:C2）—确定—存储—拟合值—确定—得到拟合一：Tt；

10、计算—计算器—将结果存储在变量中（C9）—表达式（C5 / C8）—得到C（周期变动） 实验结果与分析： 销售额时序图为：

销售额Y（3） 的时间序列图550050004500销售额Y（3）400035003000250020001510152025指数30354045

由时序图可以看出销售额Y有长期趋势，且周期的长度为4。

（1）季节指数S的计算

季节指数的计算是先用移动平均法剔除长期趋势和周期变动，然后按月(季)平均法求出季节指数

销售额Y（3） 的移动平均图550050004500销售额Y（3）变量实际拟合值移动平均长度4准确度度量平均百分误差 (MAPE)平均绝对误差 (MAD)平均偏差平方和40003500300025002000151015202530指数35404514485309627

（2）长期趋势T的计算

销售额Y（3） 的趋势分析图线性趋势模型Yt = 2736 + 39.0*t550050004500销售额Y（3）变量实际拟合值准确度度量平均百分误差 (MAPE)平均绝对误差 (MAD)平均偏差平方和40003500300025002000151015202530指数35404514478316945 销售额Y具有比较明显的上升趋势，可以用直线趋势拟合，以时间t为自变量，销售额Y为因变量，回归方程:

(3)周期变动的因素C 的计算

将序列TC除以T，即可得到周期变动因素C。

Tt 2775.06 2814.01 2852.96 2891.92 2930.87 2969.83 3008.78 3047.74 3086.69 3125.64 3164.60 C * * 0.97214 0.97534 0.96833 0.96551 0.96412 0.96750 0.98185 1.00128 1.02149

3203.55 3242.51 3281.46 3320.42 3359.37 3398.32 3437.28 3476.23 3515.19 3554.14 3593.10 3632.05 3671.01 3709.96 3748.91 3787.87 3826.82 3865.78 3904.73 3943.69 3982.64 4021.59 4060.55 4099.50 4138.46 4177.41 4216.37 4255.32 4294.28 4333.23 4372.18 4411.14 4450.09 4489.05 4528.00 4566.96 4605.91 1.02957 1.02130 1.00564 0.99886 0.99489 0.99466 0.99757 0.99916 1.00355 1.00634 1.01943 1.03483 1.04095 1.04113 1.03308 1.01909 1.00061 0.98450 0.97199 0.96459 0.97025 0.98118 0.98651 0.99292 1.00363 1.00936 1.01557 1.02474 1.03310 1.03707 1.03000 1.02775 1.03159 1.03071 1.00755 * *

实验三：皮尔曲线模型

实验目的与要求：

熟练掌握皮尔曲线模型。 实验内容：

问题：已知某地区1992~2012年的人口资料（见实验三数据.xls），试用皮尔曲线模型预测该地区2013~2017年的人口总量。

解：实验步骤：

1.将数据录入在c1,c2中

2.统计-时间序列—时间序列图-简单-确定-序列（人口总量yt）-确定 3.计算-计算器-将结果存储在变量中（c4）-表达式（1/yt）-确定

4.统计-时间序列-差分-系列（1./yt）-将差分存储在c5(zt)中-滞后（1）-确定 5. 统计-时间序列-差分-系列（zt）-将差分存储在c6(zt-1)中-滞后（1）-确定 6. 计算-计算器-将结果存储在变量中（c7）-表达式（zt/zt-1）-确定

7. 统计-时间序列-趋势分析-变量（人口总量yt）-模型类型(s曲线)-生成预

测（5）-存储（拟合值）

实验结果：

1. 模型的识别

人口总量yt 的时间序列图400037503500人口总量yt3250300027502500246810指数1214161820 由时序图可以看出人口总量基本符合皮尔曲线模型。

2. 计算yt的倒数及其一阶差分的差比率

1/yt xt/xt-1 1/yt xt/xt-1 1/yt 0.0003808 * 0.0003223 0.87772 0.0002808 0.0003789 * 0.0003153 0.94252 0.0002766 0.0003694 5.07425 0.0003074 1.11978 0.0002730 0.0003570 1.30240 0.0003016 0.74400 0.0002698 0.0003478 0.74123 0.0002950 1.12713 0.0002666 0.0003382 1.04965 0.0002893 0.86848 0.0002637 0.0003297 0.87855 0.0002847 0.80656 0.0002613 xt/xt-0.83211 1.09326 0.86414 0.89392 0.97654 0.91093 0.83671 1 该时间序列1/yt的一阶差分的差比率大致相等。综合散点图和差分分析，最后确定选用皮尔曲线作为预测模型。

3. 求模型的参数

人口总量yt 的趋势分析图S 曲线趋势模型Yt = (10**5) / (22.6405 + 17.6926*(0.924807**t))40003750人口总量yt变量实际拟合值预测曲线参数截距2479.35渐近线4416.86渐近率0.92准确度度量平均百分误差 (MAPE)平均绝对误差 (MAD)平均偏差平方和350032503000275025001994199720002003200620092012年份0.36210.747275.293 105 由图得皮尔曲线方程为：yt

22.640517.6926*0.924807t 4. 模型的预测

用所求曲线方程对22—26期人口总量作预测：

年份 2013 2014 2015 2016 2017

人口总量 3874.54 3910.64 3944.64 3976.61 4006.64

实验四：修正指数曲线模型

实验目的与要求：

熟练掌握修正指数曲线模型。 实验内容：

问题：根据实验四数据.xls中的数据资料，用修正指数曲线模型预测2013年取暖器的销售量，并说明其最高限度。

解：修正指数曲线预测模型为: ytabct (011a 将修正指数曲线化为皮尔曲线模型：tytabc1bcta1实验步骤：1.将数据录入在c1,c2中

2.统计-时间序列—时间序列图-简单-确定-序列（销售量yt）-确定 3.计算-计算器-将结果存储在变量中（c3）-表达式（1/yt）-确定

4.统计-时间序列-差分-系列（1./yt）-将差分存储在c4(zt)中-滞后（1）-确定 5. 统计-时间序列-差分-系列（zt）-将差分存储在c5(zt-1)中-滞后（1）-确定 6. 计算-计算器-将结果存储在变量中（c6）-表达式（zt/zt-1）-确定

7. 统计-时间序列-趋势分析-变量（人口总量1/yt）-模型类型(s曲线)-生成预测（1）-存储（拟合值）

8. 计算——计算器——将结果存储在变量中：c8，表达式：1/0.0000169——确定 实验结果： 1. 时序图： 销售量yt 的时间序列图600005750055000销售量yt5250050000475004500012345指数6789 由时序图可以看出取暖器销售量基本符合修正指数曲线模型。 2. 计算一阶差的一阶比率

yt xt xt/xt-1 46000 49000 51400 53320 54856 * * 3000 * 2400 0.8 1920 0.8 1536 0.8 56085 1229 57088 1003 57900 812 58563 663 0.800130 0.816110.809571 0.816501 2 由表可以看出该时间序列yt的一阶差分的差比率大致相等。

3. 计算yt的倒数及一阶差分的差比率。 1/yt 0.0000217 0.0000204 * 0.0000195 0.715953 0.0000188 0.735184 0.0000182 0.749599 0.0000178 0.760684 0.0000175 0.784203 0.0000173 0.784194 0.0000171 0.795938 zt/zt-1 * 由表可得1/yt的一阶差分的差比率大致相等，选用皮尔曲线作为预测模型。

4. 求模型的参数

将修正指数曲线化为皮尔曲线模型的结果为：

1/yt 的趋势分析图S 曲线趋势模型Yt = (10**-4) / (6.12068 - 1.88831*(0.804211**t))0.000022变量实际拟合值预测曲线参数截距0.000012渐近线0.000016渐近率0.804211准确度度量平均百分误差 (MAPE)0.0200733平均绝对误差 (MAD)0.0000000平均偏差平方和0.00000000.0000210.0000201/yt0.0000190.0000180.000017200420052006200720082009201020112012年份

拟合趋势方程

1/yt = (10**-4) / (6.12068 - 1.88831*(0.804211**t))

1104/6.12则皮尔曲线模型为： yt11.89*0.804t/6.125．模型的预测

周期 预测 10 0.0000169

即2013年取暖器的销售量的预测值为：1/0.0000169=59171（台） 其最高限度L=1/( 10/6.12)=61200 (台)

4

实验五：时间序列平滑预测法建模

实验目的与要求：

熟练掌握运用二次曲线指数平滑法建模。 实验内容：

问题：某地区统计了从1989年到2012年每年的消费品销售额（见实验五数据.xls），请用二次曲线指数平滑法预测2013年的消费品销售额。（取平滑常数0.5）

解：

实验步骤1)将数据录入工作表中，c1为年度，c2为销售额

（2）统计-时间序列-单指数平滑-变量（c2）-使用（0.5）-存储（修匀数据）-选项（K=1）-将C3命名为St(1)

(3) 统计-时间序列-单指数平滑-变量（c3）-使用（0.5）-存储（修匀数据）-选项（K=1）-将C4命名为St(2)

(4) 统计-时间序列-单指数平滑-变量（c4）-使用（0.5）-存储（修匀数据）-选项（K=1）-将C5命名为St(3)

(5) 计算-计算器-将结果存储在变量中（c6）-表达式（3*c3-3*c4+c5）-确定-将C6命名为At

（6）计算-计算器-将结果存储在变量中（c7）-表达式（3.5*c3-6*c4+2.5*c5）-确定-将C7命名为Bt

(7) 计算-计算器-将结果存储在变量中（c8）-表达式（c3-2*c4+ c5）- 确定-将

C8命名为Ct

(8) 计算-计算器-将结果存储在变量中（c9）-表达式（c6+c7+ 1/2*c8）- 确定-将C9命名为Ft+m

(1)实验结果：（1）计算t时期的单指数平滑值St

St(1)xt(1)St1

销售总额(亿元) 的单指数平滑

数据 销售总额(亿元) 长度 24 平滑常量

Alpha 0.5

准确度度量

平均百分误差 (MAPE) 10.0449 平均绝对误差 (MAD) 2.9602 平均偏差平方和 19.2161

平滑图为：

销售总额(亿元） 的平滑图单指数法变量实际拟合值平滑常量Alpha0.5准确度度量平均百分误差 (MAPE)平均绝对误差 (MAD)平均偏差平方和50销售总额(亿元）403010.04492.960219.2161201024681012141618202224指数

（2）计算t时期的双指数平滑值St(2)

St(2)St(1)(1)St1(2)

修匀数据1 的单指数平滑

数据 修匀数据1 长度 24 平滑常量

Alpha 0.5

准确度度量

平均百分误差 (MAPE) 8.9837 平均绝对误差 (MAD) 2.4715 平均偏差平方和 12.9108

平滑图为：

修匀数据1 的平滑图单指数法50变量实际拟合值平滑常量Alpha0.5准确度度量平均百分误差 (MAPE)平均绝对误差 (MAD)平均偏差平方和40修匀数据1308.98372.471512.9108201024681012141618202224指数 （3）计算t时期的三重指数平滑值St

(3)St(3)St(3)(1)St1(3)

修匀数据2 的单指数平滑

数据 修匀数据2 长度 24 平滑常量

Alpha 0.5

准确度度量

平均百分误差 (MAPE) 8.26700 平均绝对误差 (MAD) 2.11037 平均偏差平方和 8.85952

平滑图为：

修匀数据2 的平滑图单指数法454035修匀数据2变量实际拟合值平滑常量Alpha0.5准确度度量平均百分误差 (MAPE)平均绝对误差 (MAD)平均偏差平方和302520151024681012141618202224指数8.267002.110378.85952 （4）计算t时期的水平值At

At3St(1)3St(2)St(3)

（5）计算t时期的线性增量Bt

2(1)(2)(3)Bt(65)S(108)S(43)S ttt2(1)（6） 计算t时期的抛物线增量Ct

2Ct(St(1)2St(2)St(3)) 2(1)

(7)预测m时期以后，即(t+m)时期的数值Ftm

1FtmAtBtmCtm2

2

当0.5，m=1时，F2012161.7675，即2013年消费品的销售额的预测值为61.77亿元

实验六：自适应过滤法建模

实验目的与要求：

熟练掌握运用自适应过滤法建模。 实验内容：

问题：假设有以下10个样本的时间序列数据，见下表1。

表1 10个样本的时间序列数据

期数 数据xi 1 1.97 2 1.89 3 2.4 4 3.2 5 2.78 6 3.45 7 1.9 8 2.34 9 2.02 10 2.5 试用自适应过滤法预测第11期的数值。

解：实验步骤：

一．确定加权平均的权数个数

一般来说，如果时间序列原始数据具有明显的季节特征，则可确定权数个数p为季节周期长度。当数据以1年为周期进行季节变动时，若数据以月份统计，则取P=12，若数据以季度统计，则取p=4。一般地，p取在两者之间即可。 二．确定初始权数

一般情况下，初始权数取为12...p1，即以简单的算术平均数作为p初始的加权平均数。 三．计算预测值

利用xt11xt2xt1...pxtp1，选取数据样本中前p个数据x1,x2,...xp，计算第p+1期的预测值

xp11xp2xp1...px1

四．计算预测误差

第p+1期的预测误差的计算公式为：

ˆp1 ep1xp1x五．权数调整

利用ii2ket1xti1对权数进行调整：

ii2kep1xpi1

11k的取值范围为：kminp,

px2ii1max六．进行迭代调整

利用第五步得到的新一组权数i(i=1,2…p),返回第三步，进行第p+2期的预测值的计算，并产生预测误差ep2，再按ii2ket1xti1进行再一次的权数调整。 这样反复调整下去，当预测误差降为0时，系数的调整过程即结束。然而，大多

数情况下，由于序列不是随机的，最终的预测误差无法降到0，此时使用的衡量标准为均方误差

ˆt1)2/(np) MSE(xt1xtpn

当继续迭代而MSE没有进一步的改善时，即认为MSE达到最小，系数的调整过程结束。这时的权数就是最佳权数，可以用它们来计算第n+1期的预测值。

程序代码(见example6_2_1d1)：

%数据经过标准化处理的自适应过滤法的应用－－用于迭代调整的函数 function [b,SE0,SE]=example6_2_1d1(p,n,y,Y,b0,k,SE0) %b0和y必须为列向量。 yhat=zeros(n,1); for t=p+1:n

yhat(t)=Y(t,2:p+1)*b0; e=Y(t,1)-yhat(t);

b0=b0+(2*k*e).*(Y(t,2:p+1))'; end b=b0; %下面计算SE for t=p:n-1

yhat(t+1)=b'*y(t:-1:t-p+1); end

SE=(mean((y(p+1:n)-yhat(p+1:n)).^2)).^0.5;

程序代码(见example6_2_1d2)：

%自适应过滤法的应用 clear,clc

y=[1.97 1.89 2.4 3.2 2.78 3.45 1.9 2.34 2.02 2.5]';%写入数据向量y n=numel(y); p=2;

M=20;xdSE0=0.001;

b0=ones(p,1)./p;%给出初始权数 %下面求标准化常数,计算数据的标准化值。 Y=zeros(n,p+1); w=zeros(n,1); for i=p+1:n

w(i)=sum(y(i-p:i-1).^2).^0.5; Y(i,:)=y(i:-1:i-p)'/w(i); end

%计算学习常数 zz=zeros(n-2*p+1,1); for j=1:n-2*p+1

zz(j)=sum(Y(j+p:j+2*p-1,1).^2); end

k=min(1/max(zz),1/p); %迭代 SE0=inf; for i=1:M

[b,SE0,SE]=example6_2_1d1(p,n,y,Y,b0,k,SE0); xdSE=abs(SE0-SE)/SE; if xdSE<=xdSE0

break else

SE0=SE;b0=b; end end

%计算并输出预测值 yhat_1=b'*y(n:-1:n-p+1) %输出结果 i,xdSE,SE,b

运行结果： yhat_1 =

2.3907 i =

4

xdSE =

1.2582e-004 SE =

0.7059 b =

0.3075

0.8029

从运行结果可以看出，经过4次迭代，使得均方误差MSE达到最小，为 1.2582e-004 最后得到的最佳权数为：

10.3075，20.8029 第11期的预测值为：x111x102x90.3075*2.50.8092*2.022.39

实验七：ARMA模型的建模

实验目的与要求：

熟练掌握ARMA模型的建模方法。 实验内容：

问题：某市统计了2003~2012年各月的工业生产总值（见实验七数据.xls），请预测2013年各月的工业生产总值。

解：实验步骤：

2. 将数据录入工作表中

3. 数据-堆叠-列-堆叠以下列（c1-c10）-将堆叠的数据存储在当前工作表列（c11-xt） 4. 统计-时间序列-时间序列图-简单-确定-序列(c11)-确定

5. 统计-时间序列-差分-系列（xt）-将差分存储在（c12-yt）-滞后（12）-确定 6. 统计-时间序列-时间序列图-简单-确定-序列(c12-yt)-确定 7. 统计-时间序列-自相关-序列（yt）-默认滞后数-确定 7 统计-时间序列-偏自相关-序列（yt）-默认滞后数-确定

8 统计-时间序列-综合自回归移动平均-序列(yt)-自回归 (1)- 移动平均（0）-差分（0）-存储：残差，拟合值

9. 计算-概率分布-卡方-累积概率-自由度（6）-输入常量（2.16）-确定 10. 计算-概率分布-卡方-累积概率-自由度（12）-输入常量（8.48）-确定 11. 计算-概率分布-卡方-累积概率-自由度（18）-输入常量（13.66）-确定

实验结果： 时序图： 某市2003-2012年各月的工业生产总值3025工业总产值20151011224364860指数728496108120 从图中可以看出数据具有明显的周期性，作一次季节差分，ytxtxt12，差分后的结果为：

季节差分后的工业生产总值432C1210-1-211224364860指数728496108120 数据趋于平稳，平稳化后得到的序列记为{ yt}， {yt}的自相关图为：

yt 的自相关函数（包含自相关的 5% 显著限）1.00.80.60.4自相关0.20.0-0.2-0.4-0.6-0.8-1.02468101214滞后161820222426 从自相关图可以看出，自相关系数呈拖尾现象，初步判定此时间序列{yt}适合AR模型

{yt}的自相关图为：

yt 的偏自相关函数（包含偏自相关的 5% 显著限）1.00.80.60.4偏自相关0.20.0-0.2-0.4-0.6-0.8-1.02468101214滞后161820222426 从偏自相关图可以看出，偏自相关系数具有1阶截尾，考虑用AR（1）模型拟合

AR（1）模型为：

yt01yt1t

得yt0.7960.4846yt1t，

对

yt的残差进行白噪声检验

残差1的自相关图为;

残差1 的自相关函数（包含自相关的 5% 显著限）1.00.80.60.4自相关0.20.0-0.2-0.4-0.6-0.8-1.01102030405060滞后708090100 从残差1的自相关图可以看出， 自相关系数均落在2倍标准差之内，且在0附近小值波动，是白噪声序列

延迟阶数 6 12 18 LB统计量的值 2.16 8.48 13.66 P值 0.9044 0.7466 0.751 结论 拟合模型 显著有效 对{yt}再拟合AR（p）,发现也可以考虑AR（2），对AR（2）进行建模

AR（2）模型为：

yt 01yt12yt2t

由上得20.0755，对应的p值=0.441>0.05，未通过检验， 对AR（2）模型不包括常数项进行拟合：

各个参数对应的p值均小于0.05，通过检验，得AR（2）模型为：

yt0.63yt10.2456yt2t

残差的自相关图：

残差4 的自相关函数（包含自相关的 5% 显著限）1.00.80.60.4自相关0.20.0-0.2-0.4-0.6-0.8-1.02468101214滞后161820222426 从残差4的自相关图可以看出， 自相关系数均落在2倍标准差之内，且在0附近小值波动，是白噪声序列 ，优化模型

模型 AR(1) AR(2) AIC 502.057 517.830 SBC 510.103 525.877

最小信息量检验显示无论是使用AIC准则还是使用SBC准则，AR（1）模型都优于AR（2）模型，所以AR（1）是相对最优模型。 用AR（1）模型进行预测：

从周期 120 后开始的预测

95% 限制 周期 预测 下限 上限 121 31.1337 26.7807 35.4866 122 31.3160 25.3811 37.2509 123 31.4845 24.2952 38.6738 124 31.6540 23.3996 39.9085 125 31.8235 22.6264 41.0206 126 31.9930 21.9413 42.0446 127 32.1624 21.3233 43.0016 128 32.3319 20.7588 43.9050 129 32.5013 20.2381 44.7646 130 32.6708 19.7542 45.5874 131 32.8403 19.3019 46.3787 132 33.0097 18.8768 47.1426

实验八：灰色预测法建模

实验目的与要求：

熟练掌握GM(1,1)模型的建模方法。 实验内容：

问题：有如下表所示的时间序列，试建立GM(1,1)预测模型，并预测第8期的预测值。 序号 1 26.7 2 31.5 3 32.8 4 34.1 5 35.8 6 37.5 X(0)(i) 解：实验步骤：

一．GM（1,1）模型的建立

设时间序列X(0)有n个观察值，X(0)X(0)(1),X(0)(2),....X(0)(n)，通过累加生成新序列X(1)X(1)(1),X(1)(2),....X(1)(n)，则GM（1,1）模型相应的微分方程为

dX(1)aX(1) dtˆ为待估参数向量，ˆ设其中：

a，利用最小二乘法求解，可得; ˆ(BTB)1BTYn 1(1)(1)[X(1)X(2)]121[X(1)(2)X(1)(3)]1 B2``````1(1)(1)[X(n1)X(n)]12X(0)(2)(0)X(3)Yn.

.X(0)(n)求解微分方程，即可得预测模型：

X(1)(k1)X(0)(1)eak （k=0.1.2..n）

aa程序代码如下：

x0=[26.7 31.5 32.8 34.1 35.8 37.5]'; n=numel(x0);

x1=zeros(n,1); x1(1)=x0(1); for j=2:n

x1(j)=x0(j)+x1(j-1); end

B=[(-0.5)*(x1(1:n-1)+x1(2:n)) ones(n-1,1)]; Yn=x0(2:n); A1phahat=B\\Yn

x0(1)-A1phahat(2)/A1phahat(1) A1phahat(2)/A1phahat(1) 运行结果： A1phahat =

-0.0438 29.5412

ans =

701.0883

ans =

-674.3883

ˆ可得：0.0438a = 29.5412 预测模型为：

dX(1)0.0438X(1)29.5412 dt X(1)(k1)701.0883e0.0438k674.3883

二．模型检验 （一） 残差检验

ˆ按预测模型计算X(1)ˆ(0)(i)，然后计算原始序列X(0)(i)与ˆ(1)(i)累减生成X(i)，并将Xˆ(0)(i)的绝对误差及相对误差序列 Xˆ(0)(i)，i=1,2…..n (0)(i)X(0)(i)X(0)(i)(i)(0) i=1,2…..n

X(i)程序代码如下：

x0=[26.7 31.5 32.8 34.1 35.8 37.5]';

n=numel(x0); x1=zeros(n,1); x1(1)=x0(1); for j=2:n

x1(j)=x0(j)+x1(j-1); end

B=[(-0.5)*(x1(1:n-1)+x1(2:n)) ones(n-1,1)]; Yn=x0(2:n); A1phahat=B\\Yn;

x1hat=((x0(1)-A1phahat(2)./A1phahat(1)).*exp(-A1phahat(1).*(0:n-1))+A1phahat

(2)/A1phahat(1));

x0hat=zeros(n,1); x0hat(1)=x1hat(1); for k=2:n

x0hat(k)=x1hat(k)-x1hat(k-1); end

phi=abs(x0-x0hat)./x0 max(phi) 运行结果： phi =

0.0000 0.0034 0.0000 0.0049 0.0001 0.0025

ans =

0.0049

可得，相对误差小于0.5%，模型的精确度较高。

（二） 关联度检验

ˆ(0)(i)与原始序列X(0)(i)的关联系数，然后计算出关联度，根据经验，当计算出X0.5时，关联度大于0.6便满意了。

程序代码如下：

x0=[26.7 31.5 32.8 34.1 35.8 37.5]'; n=numel(x0); x1=zeros(n,1); x1(1)=x0(1); for j=2:n

x1(j)=x0(j)+x1(j-1); end

B=[(-0.5)*(x1(1:n-1)+x1(2:n)) ones(n-1,1)]; Yn=x0(2:n);

A1phahat=B\\Yn;

x1hat=((x0(1)-A1phahat(2)./A1phahat(1)).*exp(-A1phahat(1).*(0:n-1))+A1phahat(2)/A1phahat(1)); x0hat=zeros(n,1); x0hat(1)=x1hat(1); for k=2:n

x0hat(k)=x1hat(k)-x1hat(k-1); end

rho=0.5;

delta=abs(x0hat-x0);

eta=(min(delta)+rho*max(delta))./(delta+rho*max(delta)); r=mean(eta) 运行结果： r =

0.7012

1nr(i)=0.7012 ，满足0.5时检验准则r>0.60。

ni1(三)后验差检验

（1）计算原始序列标准差：

S1（2）计算绝对误差序列标准差：

[X(0)(i)Xn1(0)2] S2（3）计算方差比：

[(0)(i)]2n1CS1 S2(0)(0) （4）计算S00.6745S12.5462,ei(0)(i) 则Pp{eiS0}=1

程序代码如下：

x0=[26.7 31.5 32.8 34.1 35.8 37.5]'; n=numel(x0); x1=zeros(n,1); x1(1)=x0(1); for j=2:n

x1(j)=x0(j)+x1(j-1); end

B=[(-0.5)*(x1(1:n-1)+x1(2:n)) ones(n-1,1)]; Yn=x0(2:n); A1phahat=B\\Yn;

x1hat=((x0(1)-A1phahat(2)./A1phahat(1)).*exp(-A1phahat(1).*(0:n-1))+A1phahat(2)/A1phahat(1));

x0hat=zeros(n,1); x0hat(1)=x1hat(1); for k=2:n

x0hat(k)=x1hat(k)-x1hat(k-1); end

s1=std(x0)

s2=std(abs(x0hat-x0)) c=s2/s1

s0=0.6745*s1

ehat=abs(abs(x0hat-x0)-mean(abs(x0hat-x0))); ehat=sort(ehat); for i=1:n

if ehat(i)break end end

运行结果： s1 =

3.7750 s2 =

0.0713 c =

0.0189 s0 =

2.5462 p = 1 四、预测

模型经检验后可用于预测，预测公式为：

X(1)(k1)701.0883e0.0438k674.3883

程序代码如下：

x0=[26.7 31.5 32.8 34.1 35.8 37.5]'; n=numel(x0); x1=zeros(n,1); x1(1)=x0(1); for j=2:n

x1(j)=x0(j)+x1(j-1); end

B=[(-0.5)*(x1(1:n-1)+x1(2:n)) ones(n-1,1)]; Yn=x0(2:n); A1phahat=B\\Yn;

x1hat=((x0(1)-A1phahat(2)./A1phahat(1)).*exp(-A1phahat(1).*(n:n+1))+A1phahat(2)/A1phahat(1));

x8hat=x1hat(2)-x1hat(1)

运行结果： x8hat =

40.8303

即第8期的预测值为40.8303