灰色预测模型程序
- 格式：doc
- 大小：25.00 KB
- 文档页数：4
下载文档原格式
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
朱龙
灰色预测模型程序
>> clear
syms a b;
c=[a b]';
A=[89677,99215,109655,120333,135823,159878,182321,209401,246619,300670]; B=cumsum(A);
n=length(A);
for i=1:(n-1)
C(i)=(B(i)+B(i+1))/2; %生成累加矩阵
end
D=A;D(1)=[];
D=D';
E=[-C;ones(1,n-1)];
c=inv(E*E')*E*D;
c=c';
a=c(1);b=c(2); %预测后续数据
F=[];F(1)=A(1);
for i=2:(n+10)
F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a;
end
G=[];G(1)=A(1);
for i=2:(n+10)
G(i)=F(i)-F(i-1); %得到预测出来的数据
end
t1=1999:2008;
t2=1999:2018;
G
plot(t1,A,'O',t2,G);
/tjsj/ndsj/2011/indexch.htm
/view/2d7ab038376baf1ffc4fadb2.html
谢阳12:04:12
function gg
x0=[];
GM(x0)
end
function GM(x0) %灰色系统GM（1，1）预测
T=input('请输入T：');
x1=zeros(1,length(x0));
B=zeros(length(x0)-1,2);
yn=zeros(length(x0)-1,1);
hatx0=zeros(1,length(x0)+T);
hatx00=zeros(1,length(x0));
hatx1=zeros(1,length(x0)+T);
epsilon=zeros(length(x0),1);
omega=zeros(length(x0),1);
for i=1:length(x0)
for j=1:i
x1(i)=x1(i)+x0(j);%累加生成
end
end
x1
for i=1:length(x0)-1
B(i,1)=(-1/2)*(x1(i)+x1(i+1));
B(i,2)=1;
yn(i)=x0(i+1);
end
hatA=(inv(B'*B))*B'*yn %求a，u
for k=1:length(x0)+T
hatx1(k)=(x0(1)-hatA(2)/hatA(1))*exp(-hatA(1)*(k-1))+hatA(2)/hatA(1); end
hatx0(1)=hatx1(1);
for k=2:length(x0)+T
hatx0(k)=hatx1(k)-hatx1(k-1);%累减还原
end
hatx0 %历史数据的模拟值
for i=1:length(x0)%开始检验
epsilon(i)=x0(i)-hatx0(i);
omega(i)=(epsilon(i)/x0(i))*100;
end
c=std(epsilon)/std(x0);
p=0;
for i=1:length(x0)
if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)
p=p+1;
end
end
p=p/length(x0)
if p>0.95&c<0.35
disp('预测精度好，预测值为:')
disp(hatx0(length(x0)+T))
elseif p>0.85&c<0.5
disp('预测精度合格，预测值为:')
disp(hatx0(length(x0)+T))
elseif p>0.70&c<0.65
disp('预测精度勉强合格，预测值为:')
disp(hatx0(length(x0)+T))
elseif p<=0.7&c>0.65
disp('预测精度不合格')
end
for i=1:length(x0)
hatx00(i)=hatx0(i);
end
z=1:length(x0);
plot(z,x0,'*',z,hatx00,'-')
text(2,x0(2),'历史数据：实线')
text(length(x0)/2,hatx00(length(x0))/2,'模拟数据：虚线') end
>> x2=[0 24 49 73 98 147 196 245 294 342];
y2=[33.46 32.47 36.06 37.96 41.04 40.09 41.26 42.17 40.36 42.73]; x1=[0 34 67 101 135 202 259 336 404 471];
y1=[15.18 21.36 25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75]; x3=[0 47 93 140 186 279 372 465 558 651];
y3=[18.98 27.35 34.86 38.52 38.44 37.73 38.43 43.87 42.77 46.22]; subplot(1,3,1);plot(x1,y1,'r*')
subplot(1,3,2);plot(x2,y2,'b*')
subplot(1,3,3);plot(x3,y3,'g*')
y1=-0.00034*(x1).^2+0.1972*(x1)+14.742;
y2=-0.000138*(x2).^2+0.07182*x2+32.9161;
y3=0.03035*x3+28.24622;
-
>> x1=[0 34 67 101 135 202 259 336 404 471];
y1=[15.18 21.36 25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75]; p=polyfit(x1,y1,2)
p =
-0.0003 0.1971 14.7416
>> x=0:20:471;
>> y=polyval(p,x);
>> plot(x1,y1,'*r',x,y,'-b')