基础的统计学预测方法分析,内容参考国防工业出版社-司守奎,孙玺菁主编-《数学建模算法与应用(第三版)》。本文结合实际应用对文章内容进行了提取,结合matlab算法进行程序编写。

本文所涉及的所有代码内容可通过百度网盘进行下载:下载地址

统计学预测方法Statistical data prediction methods

Contents

1.典型微分方程、指数平滑法和灰色预测

Typical differential equation, exponential smoothing method and grey prediction

1.1.典型微分方程

微分方程是数学建模的重要方法,许多实际问题的数学描述将导致微分方程的定解问 题。列方程常见的类型有3种:(1)按规律直接列方程;(2)微元分析与区域积分;(3) 模拟近似法。 在各类学科都有很多经过实践检验的规律和定律可以列写微分方程,如牛顿运动定律、 基尔霍夫电流电压定律、物质放射性规律、曲线切线性质等。

例(物体冷却过程):牛顿冷却定律指出,当物体表面与周围存在温度差时,单位时间 内从单位面积散失的热量与温度差成正比,比例系数为热传递系数,记为k,试建立 相关传热模型。我们使用2018年全国大学生数学建模国赛A题第一问进行一些分析:

(2018数学建模国赛A.1)高温作业服可避免人们在高温环境下作业受到灼伤。防护服分 为4层,其中Ⅰ、Ⅱ、Ⅲ层为织物制造,第Ⅳ层为织物与皮肤之间的空气间隙;第Ⅰ 层直接与外部环境接触,假人体内恒温为 37℃,外界温度为75℃。试建立非稳态传热 模型,反应不同时间节点的传热情况。 数据见文件:data1。

利用文件中附件1信息可以建立逻辑严密的热传导模型,读者不妨思考,可不可以只 用附件1有关厚度信息,忽略边界效应直接利用牛顿冷却建立微分方程,利用优化算法 去逼近附件2中数据呢?答案是可行的。 设最外层防护服与空气间热传递系数为k0,Ⅰ、Ⅱ、Ⅲ、Ⅳ层内部热传递系数分别为 k1,k2,k3,k4,Ⅳ与皮肤间热传递系数为k5,问题变为建立模型,求Ⅳ靠近皮肤的表面 温度,使其值逼近附件2。 可编写如下程序建立简单的传热模型:

clear,clc
k0=1;k1=0.99;k2=0.98;k3=0.9;k4=0.99;k5=0.01;% 传热系数
d=[0.6,0.6,3.6,0.6];% 厚度设置
d1=cumsum(d);
dx=[1e-2,1e-2,1e-1,1e-2];% 厚度微元设置
tw=75;% 外部温度75
tn=37;% 假人温度37
u1=tn*ones(1,length(0:dx(1):d1(1)));
u2=tn*ones(1,length(d1(1):dx(2):d1(2)));
u3=tn*ones(1,length(d1(2):dx(3):d1(3)));
u4=tn*ones(1,length(d1(3):dx(4):d1(4)));% 防护服温度初始化
tm=5400;% 总时间
u=zeros(1,tm);
for t=1:tm
    tic
    for i=1:length(0:dx(1):d1(1)) % 第一层防护服温度迭代
        if(i==1)
            deltau=-k0*(u1(i)-tw)-k1*(u1(i)-u1(i+1));
            u1(i)=u1(i)+deltau;
        elseif(i==length(0:dx(1):d1(1)))
            deltau=-k1*(u1(i)-u1(i-1))-k2*(u1(i)-u2(1));
            u1(i)=u1(i)+deltau;
        else
            deltau=-k1*(u1(i)-u1(i-1))-k1*(u1(i)-u1(i+1));
            u1(i)=u1(i)+deltau;
        end
    end
    for i=1:length(d1(1):dx(2):d1(2)) % 第二层防护服温度迭代
        if(i==1)
            deltau=-k1*(u2(i)-u1(end))-k2*(u2(i)-u2(i+1));
            u2(i)=u2(i)+deltau;
        elseif(i==length(d1(1):dx(2):d1(2)))
            deltau=-k2*(u2(i)-u2(i-1))-k3*(u2(i)-u3(1));
            u2(i)=u2(i)+deltau;
        else
            deltau=-k2*(u2(i)-u2(i-1))-k2*(u2(i)-u2(i+1));
            u2(i)=u2(i)+deltau;
        end
    end
    for i=1:length(d1(2):dx(3):d1(3)) % 第三层防护服温度迭代
        if(i==1)
            deltau=-k2*(u3(i)-u2(end))-k3*(u3(i)-u3(i+1));
            u3(i)=u3(i)+deltau;
        elseif(i==length(d1(2):dx(3):d1(3)))
            deltau=-k3*(u3(i)-u3(i-1))-k4*(u3(i)-u4(1));
            u3(i)=u3(i)+deltau;
        else
            deltau=-k3*(u3(i)-u3(i-1))-k3*(u3(i)-u3(i+1));
            u3(i)=u3(i)+deltau;
        end
    end
    for i=1:length(d1(3):dx(4):d1(4)) % 第四层防护服温度迭代
        if(i==1)
            deltau=-k3*(u4(i)-u3(end))-k4*(u4(i)-u4(i+1));
            u4(i)=u4(i)+deltau;
        elseif(i==length(d1(3):dx(4):d1(4)))
            deltau=-k4*(u4(i)-u4(i-1))-k5*(u4(i)-tn);
            u4(i)=u4(i)+deltau;
        else
            deltau=-k4*(u4(i)-u4(i-1))-k4*(u4(i)-u4(i+1));
            u4(i)=u4(i)+deltau;
        end
    end
    toc
    u(t)=u4(end);% 第四层靠近皮肤表面温度
end
此时需要做的工作是利用最小二乘编写优化算法,优化各层传热系数,使其逼近附件2值即可。
例(导弹追击问题):设为一导弹发射点A(0,0),(0,10)处有目标物向右以恒定速率v运动,0时刻处发射一枚导弹以2v的速度追赶目标物,试通过编程,输入b、v,输出导弹动态路线。
clear,clc
b=input('请输入初始时刻目标物的高度:');
v=input('请输入目标物运动速度大小:');
m=0;% 目标物横坐标
q=0;% 导弹横坐标
p=0;% 导弹纵坐标
dt=1;% 时间微元
for i=1:1000
    plot(m,b,'b*-',p,q,'rx--');% 绘制当前导弹和目标物位置
    axis([0,p+10,0,b+1]);
    hold on
    pause(0.2);% 动态显示
    m=m+v*dt;% 目标物横坐标运动
    d=sqrt((m-p)^2+(b-q)^2);% 目标物与导弹间距
    p=p+(2*v/d)*(m-p)*dt;% 导弹横坐标运动
    q=q+(2*v/d)*(b-q)*dt;% 导弹纵坐标运动
    if d<0.05% 当两者距离小于某值时停止
        break
    end
end
hold off
% 例(传染病SIR模型):假设人群分为健康者,病人,和病愈后具有免疫力不得病三类设任意时刻t,这三类人群占总人口比例分别为s(t),i(t),r(t),设病人的日接触率为lamda,日治愈率为miu,则有传染强度delta=lamda/miu,若人口总数为N,则有微分方程组:
s(t)+i(t)+r(t)=1;
dr/dt=miu*i;
ds/dt=-lamda*s*i;
di/dt=lamda*s*i-miu*i;
利用matlab仿真此类传染病的发病过程。
clear,clc
N=1e6;% 人口总数
s0=0.98;% 健康者初始比例
i0=0.02;% 病人初始比例
r0=0;% 免疫者初始比例
dt=1;% 单位时间/天
lamda=3;% 病人日接触率
miu=0.5;% 日治愈率
s=s0;i=i0;r=r0;
for k=1:1e2 % 总天数
    plot(k,N*s,'b*-',k,N*i,'rx--',k,N*r,'bx--');% 绘制当前各类人群人数
    axis([0,k+10,0,N]);
    legend('健康者','病人','免疫者');
    hold on
    pause(0.2);% 动态显示
    i=i+(lamda*s*i-miu*i)*dt;
    s=s+(-lamda*s*i)*dt;
    r=r+(miu*i)*dt;
end
hold off
利用该模型可以简单模拟传染病发病过程。
对于常微分方程,matlab提供了几个解常微分方程数值解的函数,如:ode45,ode23,ode113等,其中ode45采用四五阶龙格-库塔法(RK方法),是解非刚性微分方程的首选方法,ode23采用23阶RK方法,ode113采用多步法,效率一般比ode45高。
例(Logistic-弱肉强食模型):对于一种群,当存在自然资源,环境条件限制时,种群增长率将会随种群数量增加而减小,设t时刻种群数量为x(t),种群增长率r(x)为x的线性函数,即:r(x)=R-s*x;则不考虑其他条件时,该种群数量增长可表示为:
dx/dt=r(x)*x;
当自然界存在捕食者和被捕食者两种群时,设t时刻两种群数量分别为x1(t),x2(t),种群增长率分别为r1(x1)=(1e2-x1)/(50),r2(x2)=(1e3-x2)/(100);,捕杀量为0.2*x1*x2;利用matlab仿真两种群数量变化。
clear,clc
x10=5;% 捕食者初值
x20=500;% 被捕食者初值
dx=@(t,x)[(1e2-x(1))/(50)*x(1)+0.01*x(1)*x(2);
    (1e3-x(2))/(100)*x(2)-0.005*x(1)*x(2)];% 列写微分方程
[t,x]=ode45(dx,[0,100],[x10,x20]);
plot(t,x(:,1),'b*-',t,x(:,2),'rx--')
legend('捕食者','被捕食者')

1.2.指数平滑法

本节及之后讨论对时间序列数据的预测问题。假设预测时,离预测点越近的点作用越 大,随着时间变化权重以指数形式变化,则可以采用指数平滑法。一般情况下,时间 序列可以分解为水平部分(均值),趋势部分(上升或下降),季节性部分(周期重复), 随机波动(白噪声)。

当数据没有趋势和季节性时,可以采用一次指数平滑法;当数据有趋势无季节性时, 可以采用二次指数平滑法;当数据既有趋势又有季节性时,采用三次指数平滑法。 对于数据噪声,可以采用信号处理相关方法,如小波变换进行滤除。

现做以下规定:y(t)为t时刻数据,y'(t)为t时刻预测值,sn(t)为t时刻n次指数平滑值,alpha为记忆衰减因子。

一次指数平滑:y'(t+1)=alpha*y(t)+(1-alpha)*y'(t);s1(t)=y'(t);

二次指数平滑:y'(t+T)=a(t)+b(t)*T;其中,a(t)=2*s1(t)-s2(t);

b(t)=alpha/(1-alpha)*(s1(t)-s2(t));s2(t)=alpha*s1(t)+(1-alpha)*s2(t-1);

三次指数平滑:y'(t+T)=a(t)+b(t)*T+c(t)*T^2;其中,a(t)=3*s1(t)-3*s2(t)+s3(t);

b(t)=alpha/(2*(1-alpha)^2)*((6-5*alpha)*s1(t)-2*(5-4*alpha)*s2(t)+(4-3*alpha)*s3(t)); s3(t)=alpha*s2(t)+(1-alpha)*s3(t-1);

alpha取值标准:水平趋势,0.05~0.2;有波动,长期趋势变化不大,0.1~0.4;波动 大,长期趋势变化大,0.6~0.8;上升或下降的发展趋势,0.6~1。 本节利用指数平滑法对具有上升趋势和周期性的数据进行预测:

clc,clear
x=0:1e-1:20;y=2*exp(-0.16*x).*sin(2*x)+1/2*x;% 原始信号
alpha=0.6;% alpha值
n=length(y);% 原始信号长度
s10=mean(y(1:3));% 取y前三个数据均值作为各次指数平滑初值
s20=s10;
s30=s10;
s1=alpha*y(1)+(1-alpha)*s10;% 各次指数平滑值初始化
s2=alpha*s1(1)+(1-alpha)*s20;
s3=alpha*s2(1)+(1-alpha)*s30;
for i=2:n+20 % 预测原始信号后20位数据值
    if i>n
        y=[y,a+b+c];% T=1时,三次指数平滑预测值,此处进行滚动预测
    end
    s1=alpha*y(i)+(1-alpha)*s1;% 各次指数平滑值迭代
    s2=alpha*s1+(1-alpha)*s2;
    s3=alpha*s2+(1-alpha)*s3;
    a=3*s1-3*s2+s3;
    b=alpha/(2*(1-alpha)^2)*((6-5*alpha)*s1-2*(5-4*alpha)*s2+(4-3*alpha)*s3);
    c=alpha^2/(2*(1-alpha)^2)*(s1-2*s2+s3);
end
plot(1:n+20,y)% 绘制预测图

1.3.灰色预测

灰色预测适合处理小样本,中短期,指数变化的数据。设有时间序列x(n),n=1,2,...,N ,计算序列级比lamda(k)=x(k-1)/x(k),k=2,3,...,N,若对于任意 lamda(k) 属于 (exp(-2/(N+1)),exp(2/(N+1))),则可进行GM(1,1)预测。

现令:x1(n)=cumsum(x(n)),n=1,2,...,N为x(n)一次累加序列;

x0(n)=diff(x(n)),n=1,2,...,N-1为x(n)一次累减序列:

z(n)=1/2*(x1(2:end)+x1(1:end-1))为x1(n)均值序列。

GM(1,1)预测: 求解:[a,b]'=(B'*B)\B'*x(1:N)';

其中,B=[-z(1:N)',ones(N,1)];

则可预测得到: x1'(k+1)=(x(1)-b/a)*exp(-a*k)+b/a;x'(k+1)=x1'(k+1)-x1'(k);

进行GM(1,1)相对误差检验:delta(k)=abs(x'(k)-x(k))/x(k),k=1,2,...,N,当 delta(k)<0.1时,达到较高精度;delta(k)<0.2时,达到一般精度。

进行GM(1,1) 级比偏差值检验:rho(k)=1-((1-1/2*a)/(1+1/2*a))*lamda(k);rho(k)<0.1时,达到 较高精度;rho(k)<0.2时,达到一般精度.

GM(2,1)预测: 对于非单调摆动发展序列或有饱和s型序列,可以考虑建立GM(2,1),DGM,Verhulst模型:

称d2(x1(t))/xt2+a1*dx1(t)/dt+a2*x1(t)=b为GM(2,1)的白化方程。

令: B=[-x(2:end)',-z(1:end)',ones(1:N-1)'];Y=[x0(1:end)']; 则有GM(2,1)模型最小二乘估计:[a1,a2,b]'=(B'*B)\B'*Y; 得到参数后利用边界条件x1(1),x1(end)求解白化方程即可进行预测。

DGM(2,1)预测: 令B=[-x(2:end)',ones(1:N-1)'];Y=[x0(1:end)'];得到:[a,b]'=(B'*B)\B'*Y; 则可预测得到:x1'(k+1)=(b/a-x(1)/a)*exp(-a*k)+b/a*k+(1+a)/a*x(1)-b/a^2; Verhulst模型: 令B=[-z(1:end)',[z(1:end).^2]'];Y=[x(2:end)'];得到:[a,b]'=(B'*B)\B'*Y; x1'(k+1)=a*x(1)/(b*x(1)+(a-b*x(1))*exp(a*k)); 本节给出matlab简要计算各灰色预测模型的代码: 

clear,clc
% 生成各类型序列
x=0:0.1:10;
f1=diff(exp(-0.02*x)+1);% % 指数下降类序列,适合GM(1,1)预测
f2=diff(1e2*exp(x)-1e1*exp(0.001*x)+2);% 双指数综合,可选用GM(2,1)预测
f3=diff(2.7*exp(-0.4*x)+4*x+1);% 指数与一次综合,可选用DGM(2,1)预测
f4=diff(1./(1+50*exp(-0.01*x)));% s型曲线,可选用Verhulst预测
f=[f1;f2;f3;f4];
% 预测数据,与原始数据作比较
for i=1:4
    data=f(i,:);
    n=length(data);% 原始数据长度
    data0=diff(data);
    data1=cumsum(data);
    z_data=1/2*(data1(2:end)+data1(1:end-1));
    if(i==1)% GM(1,1)
        B=[-z_data(1:end)',ones(n-1,1)];
        Y=[data(2:end)]';
        u=(B'*B)\B'*Y;
        pre=@(k)(data(1)-u(2)/u(1))*exp(-u(1)*(k-1))+u(2)/u(1);
        fp=pre(1:n);fp=diff(fp);
    end
    if(i==2)% GM(2,1)
        B=[-data(2:end)',-z_data(1:end)',ones(n-1,1)];
        Y=data0';
        u=(B'*B)\B'*Y;
        syms x_un(t)
        x_un=dsolve(diff(x_un,2)+u(1)*diff(x_un)+u(2)*x_un==u(3), ...
            x_un(0)==data1(1),x_un(n-1)==data1(end));% 求解白化方程
        xt=vpa(x_un,6);
        fp=double(subs(x_un,t,0:n-1));
    end
    if(i==3)% DGM(2,1)
        B=[-data(2:end)',ones(n-1,1)];
        Y=data0';
        u=(B'*B)\B'*Y;
        pre=@(k)(u(2)/u(1)^2-data(1)/u(1))*exp(-u(1)*(k-1))+...
        u(2)/u(1)*(k-1)+(1+u(1))/u(1)*data(1)-u(2)/u(1)^2;
        fp=pre(1:n);fp=diff(fp);
    end
    if(i==4)% Verhulst
        B=[-z_data',z_data'.^2];
        Y=data(2:end)';
        u=(B'*B)\B'*Y;
        pre=@(k)u(1)*data(1)./(u(2)*data(1)+(u(1)-u(2)*data(1))*exp(u(1)*(k-1)));
        fp=pre(1:n);fp=diff(fp);
    end
    figure(i)
    plot(data,'b*-')
    hold on
    plot(fp,'r*-')
    hold off
    legend('实际值','预测值')
end
从预测结果中可以清楚的观察到部分结果预测并不准确,灰色预测适合小样本数据,当样本数据足够多时可选择其他时间序列预测方法,此外,样本变化幅度不能过大,在预测前需要对数据进行检验,灰色预测只能进行中短期预测,长期预测结果将出现较大误差。

2.马尔可夫预测Markov predicted

某一系统在已知现在情况的条件下,系统未来时刻的情况只与现在的情况有关,而与 历史无直接关系,则该系统可称其位马尔可夫系统。 有马氏链模型:某序列x(n)满足: p{x(n+1)=j|x(n)=i(n),x(n-1)=i(n-1),...,x(1)=i(1)}=p{x(n+1)=j|x(n)=i(n)}; 即未来时刻值只与现在值有关的序列。 时齐的马氏链:由一个状态转移到另一个状态的概率只与时间间隔有关,与起始时刻 无关。称以m步转移概率Pij(m)位元素的矩阵P(m)为马氏链的m步转移矩阵。 例:设某系统状态空间为:

E=[1,2,3,4];
% 有观察数据:
a=num2str([4,3,2,1,4,3,1,1,2,3,2,1,2,3,4,4,3,3,1,1,1,3,3,2,1,2,2,2,4,4,2,3,2,3,1,1,2,4,3,1]);
% 则有该马氏链一步转移矩阵P(1):
f=zeros(length(E));
for i=E
    for j=E
        f(i,j)=length(strfind(a,num2str([i j])));% 寻找a中出现[i,j]的次数
    end
end
p=sum(f,2);
for i=1:4
    f(i,:)=f(i,:)/p(i);
end% 一步转移矩阵

3.ARIMA预测ARIMA Forecast

3.1.平稳性Daniel检验

设有一组观察数据:(x(1),y(1)),(x(2),y(2)),...,(x(n),y(n)),设{x(1),x(2),...,x(n)} 的秩统计量为{R(1),R(2),...,R(n)}(秩统计量可理解为若将总数居从小到大排,该数 据排到第几),{y(1),y(2),...,y(n)}的秩统计量为{S(1),S(2),...,S(n)}。令: d={R(1)-S(1),R(2)-S(2),...,R(n)-S(n)}则观察数据有Spearman相关系数: Qxy=1-6/(n*(n^2-1))*sum(d.^2); 做假设检验:H0:rho(xy)=0,H1:rho(xy)~=0,其中rho(xy)为x,y总体相关系数。 令T=Qxy*sqrt(n-2)/sqrt(1-Qxy^2);,当abs(T)<=t(alpha/2)(n-2)时接受H0。 当H0检验通过时,说明观察数据序列平稳,不通过时说明存在上升或下降趋势。

clear,clc
% 设有数据列:
n=1000;
x=cumsum(ones(1,n));
y=randn(1,n);% 高斯噪声列
ry=sort(y);Ry=zeros(1,n);
for i=1:n
    [~,r]=find(y==ry(i));
    Ry(i)=r;% 获得y数据列秩统计量
end
d=x-Ry;
Qxy=1-6/(n*(n^2-1))*sum(d.^2);
T=Qxy*sqrt(n-2)/sqrt(1-Qxy^2);
t_0=tinv(0.975,n-2);% 计算上alpha/2分位数
% 由于abs(T)<t_0,则通过假设检验,数据列平稳。

3.2.自相关系数与偏自相关系数

clear,clc
% 设有数据列:
n=1000;
x=cumsum(ones(1,n));
y=randn(1,n);% 高斯噪声列
y_mean=mean(y);% 数据y平均值
y_0=y-y_mean;% 原始序列减去其均值
% 则自相关系数可表示为:
rho=zeros(1,n);
for k=0:n-1
    rho(k+1)=sum(y_0(1:n-k).*y_0(k+1:n))/sum(y_0.^2);
end
% 当序列为噪声序列时,自相关系数rho(k+1)~N(1,1/n)。
% 偏自相关系数计算:
phi=zeros(1,n-1);
for k=0:n-2
    rho1=rho(2:k+2);
    kp=length(rho1);
    rho2=zeros(kp);
    for k0=1:kp
        rho2(k0,k0:end)=rho(1:kp-k0+1);
        if(k0>=2)
            rho2(k0,1:k0-1)=fliplr(rho1(1:k0-1));
        end
    end
    phi1=rho2\rho1';phi(k+1)=phi1(end);
end

3.3.ARIMA

AR(p)模型:p阶自回归,y(t)=miu+sum(r(i)*y(t-i))(i=1:p)+et,其中,et为残差, miu,r为待定系数。 MA(q)模型:q阶滑动回归:y(t)=miu+sum(r(i)*e(t-i))(i=1:q)+et,其中,et为残差, miu为待定系数。 ARIMA中I指差分,ARIMA模型处理的数据为平稳数据,当数据有上升下降趋势时,需 要通过差分获得平稳数据列,若令需要差分的次数为d,则有完整的ARIMA(p,q,d)模 型:y(t)=miu+sum(r(i)*y(t-i))(i=1:p)+sum(r(i)*e(t-i))(i=1:q)+et。 本节给出matlab ARIMA模型预测的简单方法: (程序可见文件:ARIMA_model_prediction.m) ARIMA模型预测ARIMA model prediction

clear,clc
% 原始数据列:
n=1000;
x=cumsum(ones(1,n));
y=1/n*x+sin(x/(4*pi))/10+randn(1,n);% 具有上升趋势,小幅周期波动,噪声的数据列
% 判断数据列平稳性并做差分处理
dm=4;% 约束最大进行4次差分
yd=y;% 差分数列初始化
for dp=1:dm+2
    n1=length(yd);
    ry=sort(yd);Ry=zeros(1,n1);
    for i=1:n1
        [~,r]=find(yd==ry(i));
        Ry(i)=r;% 获得y数据列秩统计量
    end
    xp=cumsum(ones(1,n1));
    d=xp-Ry;
    Qxy=1-6/(n1*(n1^2-1))*sum(d.^2);
    T=Qxy*sqrt(n1-2)/sqrt(1-Qxy^2);
    t_0=tinv(0.975,n1-2);% 计算上alpha/2分位数
    if(abs(T)<=t_0)
        break
    elseif(dp==dm+2)
        disp('原始数据列无法差分平稳化,可能不适合ARIMA模型')
        break
    else
        yd=diff(yd);
    end
end
d_p=dp-1;% 最终差分次数
% 自相关系数,偏自相关系数计算
y_mean=mean(yd);% 数据yd平均值
y_0=yd-y_mean;% 序列减去其均值
rho=zeros(1,n1);
for k=0:n1-1
    rho(k+1)=sum(y_0(1:n1-k).*y_0(k+1:n1))/sum(y_0.^2);
end
phi=zeros(1,n1-1);
for k=0:n1-2
    tic
    rho1=rho(2:k+2);
    kp=length(rho1);
    rho2=zeros(kp);
    for k0=1:kp
        rho2(k0,k0:end)=rho(1:kp-k0+1);
        if(k0>=2)
            rho2(k0,1:k0-1)=fliplr(rho1(1:k0-1));
        end
    end
    phi1=rho2\rho1';phi(k+1)=phi1(end);
    toc
end
定阶预测
ARMA(p,q)预测一般需满足:数据自相关系数(ACF)q阶后衰减趋于0,偏自相关系数(PACF)p阶后衰减趋于0.
如果样本ACF和样本PACF在最初k阶明显大于2倍标准差,而后几乎95%的系数都落在2倍标准差的范围内,且非零系数衰减为小值波动的过程非常突然,通常视为k阶截尾;如果有超过5%的样本相关系数大于2倍标准差,或者非零系数衰减为小值波动的过程比较缓慢或连续,通常视为拖尾。
rho_std=std(rho);% 自相关系数标准差
phi_std=std(phi);% 偏自相关系数标准差
p=0;q=0;
for i=rho
    if(abs(i)>2*rho_std)
        p=p+1;
    else
        break
    end
end
for i=phi
    if(abs(i)>2*phi_std)
        q=q+1;
    else
        break
    end
end
以上程序简单进行了定阶,通常利用贝叶斯信息等可以更准确的进行定阶,限于篇幅本文不再赘述。
利用matlab内置函数进行预测;
AR_Order=p;
MA_Order=q;
AR=arima(AR_Order, d_p, MA_Order);
EstMdl=estimate(AR,y');
step=100;% 预测步数
fore_y=forecast(EstMdl,step,'Y0',y');
figure(1)
plot([y,fore_y'],'r*-')
hold on
plot(y,'b*-')
hold off
legend('预测值','实际值')
% 由结果分析可知,ARIMA模型适合短期的预测,且阶数的选择会影响预测结果。

4.插值与拟合Interpolation and fitting

4.1.插值

若数据有多项式特征,则可以用多项式函数进行插值。可利用待定系数法确定插值多 项式。

clear,clc
% 设有待插值原始数据列:
x0=(1:6)';y0=[16,18,21,17,15,12]';
A=vander(x0);% 返回x0范德蒙德矩阵
p=A\y0;% 求出多项式系数
x=0:0.1:6;
yh=polyval(p,x);plot(x,yh)% 绘制函数插值图像
% 用多项式做插值函数,随着插值节点的增加,插值多项式的次数越多,高次插值不但
% 计算复杂而且会产生龙格振荡现象。因此实际情况中一般会采用其他插值方法。
% 分段线性插值:将相邻节点用直线连接起来;
% 三次样条插值:插值函数在每个子区间均为三次多项式,且满足一定边界条件。
% 本节介绍matlab中一维插值函数:
clear,clc
% 设有待插值原始数据列:
x0=(1:6)';y0=[16,18,21,17,15,12]';
method='spline';% 设置插值方法,可以为'nearest'最近邻插值;'linear'线性插值;
% 'spline'三次样条插值;'cubic'立方插值
xq=0:0.1:6;% 待求点位置
vq=interp1(x0,y0,xq,method);plot(xq,vq)
% 二维插值函数interp2与interp1用法基本相同。

4.2.拟合

已知一组二维数据(x(i),y(i)),i=1,2,...,n,要寻找一个函数y=f(x)使f(x)在某准 则下与所有数据点最为接近。此时可以采用偏差平方和最小准则。即: min J=sum((f(1:n)-y(1:n)).^2);这一原则称为最小二乘元则。 线性最小二乘: 给定一个线性无关的函数系{phi(k)|,k=1,2,...,m},若拟合函数可以以其线性组合 形式出现:f(x)=sum(a(1:m).*phi(1:m)),则f(x)为关于系数{a(k)|,k=1,2,...,m} 的线性函数。 对于拟合函数:f(x)=a(1)*phi(1)+a(2)*phi(2)+,...,+a(m)*phi(m); 记:R=[phi(1:m)|x(1);phi(1:m)|x(2);...;phi(1:m)|x(n)];A=[a(1:m)];Y=[y(1:n)]; 有:A=(R'*R)\R'*Y。 matlab求约束线性最小二乘解:

clear,clc
x=[3,5,6,7,4,8,5,9]';y=[4,9,5,3,8,5,8,5]';% 观测数据
c=[exp(x),log(x)];% 线性无关函数系
a=lsqlin(c,y);% 拟合参数
% matlab求多项式拟合:
clear,clc
x=[3,5,6,7,4,8,5,9]';y=[4,9,5,3,8,5,8,5]';% 观测数据
p=polyfit(x,y,3);% 拟合3次多项式,返回值p为多项式系数,从高次幂到低次幂
% matlab fit和fittype函数
clear,clc
x=[3,5,6,7,4,8,5,9]';y=[4,9,5,3,8,5,8,5]';% 观测数据
g=fittype('a*exp(x)+b*log(x)','dependent',{'y'},'independent',{'x'});
f=fit(x,y,g,'StartPoint',rand(1,2));


Published with MATLAB® R2022b

04-26 20:47