基础的统计学预测方法分析，内容参考国防工业出版社-司守奎，孙玺菁主编-《数学建模算法与应用（第三版）》。本文结合实际应用对文章内容进行了提取，结合matlab算法进行程序编写。
本文所涉及的所有代码内容可通过百度网盘进行下载：下载地址
统计学预测方法Statistical data prediction methods
Contents
- 1.典型微分方程、指数平滑法和灰色预测
- 1.1.典型微分方程
- 1.2.指数平滑法
- 1.3.灰色预测
- 2.马尔可夫预测Markov predicted
- 3.ARIMA预测ARIMA Forecast
- 3.1.平稳性Daniel检验
- 3.2.自相关系数与偏自相关系数
- 3.3.ARIMA
- 4.插值与拟合Interpolation and fitting
- 4.1.插值
- 4.2.拟合
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:tmticfor 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;elsedeltau=-k1*(u1(i)-u1(i-1))-k1*(u1(i)-u1(i+1));u1(i)=u1(i)+deltau;endendfor 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;elsedeltau=-k2*(u2(i)-u2(i-1))-k2*(u2(i)-u2(i+1));u2(i)=u2(i)+deltau;endendfor 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;elsedeltau=-k3*(u3(i)-u3(i-1))-k3*(u3(i)-u3(i+1));u3(i)=u3(i)+deltau;endendfor 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;elsedeltau=-k4*(u4(i)-u4(i-1))-k4*(u4(i)-u4(i+1));u4(i)=u4(i)+deltau;endendtocu(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:1000plot(m,b,'b*-',p,q,'rx--');% 绘制当前导弹和目标物位置axis([0,p+10,0,b+1]);hold onpause(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% 当两者距离小于某值时停止breakend 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 onpause(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>ny=[y,a+b+c];% T=1时，三次指数平滑预测值，此处进行滚动预测ends1=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:4data=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);endif(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));endif(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);endif(i==4)% VerhulstB=[-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);endfigure(i)plot(data,'b*-')hold onplot(fp,'r*-')hold offlegend('实际值','预测值')
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=Efor j=Ef(i,j)=length(strfind(a,num2str([i j])));% 寻找a中出现[i,j]的次数end end p=sum(f,2); for i=1:4f(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-1rho(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-2rho1=rho(2:k+2);kp=length(rho1);rho2=zeros(kp);for k0=1:kprho2(k0,k0:end)=rho(1:kp-k0+1);if(k0>=2)rho2(k0,1:k0-1)=fliplr(rho1(1:k0-1));endendphi1=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+2n1=length(yd);ry=sort(yd);Ry=zeros(1,n1);for i=1:n1[~,r]=find(yd==ry(i));Ry(i)=r;% 获得y数据列秩统计量endxp=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)breakelseif(dp==dm+2)disp('原始数据列无法差分平稳化，可能不适合ARIMA模型')breakelseyd=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-1rho(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-2ticrho1=rho(2:k+2);kp=length(rho1);rho2=zeros(kp);for k0=1:kprho2(k0,k0:end)=rho(1:kp-k0+1);if(k0>=2)rho2(k0,1:k0-1)=fliplr(rho1(1:k0-1));endendphi1=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=rhoif(abs(i)>2*rho_std)p=p+1;elsebreakend
end
for i=phiif(abs(i)>2*phi_std)q=q+1;elsebreakend
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
