12 3 / 3 页下一页

查看: 17202|回复: 26

[问题求助] 灰色预测Matlab 程序

[复制链接]

字体大小: 正常放大

kelimasa

7 主题	5 听众	188 积分

升级 44%

TA的每日心情

	开心 2012-9-10 21:57

签到天数: 53 天

[LV.5]常住居民I

电梯直达

1^#

发表于 2011-12-15 09:26 |只看该作者 |倒序浏览

|招呼Ta 关注Ta

标签：灰色模型 gm(1 1) 二次拟合 matlab 分类：技术点滴

%by allen @ 红嘴海鸥
%灰色模型预测是在数据不呈现一定规律下可以采取的一种建模和预测方法，其预测数据与原始数据存在一定的规律相似性

%下面程序是灰色模型GM(1,1)程序二次拟合和等维新陈代谢改进预测程序,matlab6.5 ,使用本程序请注明，程序存储为gm1.m

%x = [5999,5903,5848,5700,7884];gm1(x);  测试数据

%二次拟合预测GM(1,1)模型
function  gmcal=gm1(x)
sizexd2 = size(x,2);
%求数组长度

k=0;
for y1=x
k=k+1;
if k>1
      x1(k)=x1(k-1)+x(k);
      %累加生成
      z1(k-1)=-0.5*(x1(k)+x1(k-1));
      %z1维数减1，用于计算B
      yn1(k-1)=x(k);
else
      x1(k)=x(k);
end
end
%x1,z1,k,yn1

sizez1=size(z1,2);
%size(yn1);
z2 = z1';
z3 = ones(1,sizez1)';

YN = yn1'; %转置
%YN

B=[z2 z3];
au0=inv(B'*B)*B'*YN;
au = au0';
%B,au0,au

afor = au(1);
ufor = au(2);
ua = au(2)./au(1);
%afor,ufor,ua
%输出预测的  a u 和 u/a的值

constant1 = x(1)-ua;
afor1 = -afor;
x1t1 = 'x1(t+1)';
estr = 'exp';
tstr = 't';
leftbra = '(';
rightbra = ')';
%constant1,afor1,x1t1,estr,tstr,leftbra,rightbra

strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+',leftbra,num2str(ua),rightbra)
%输出时间响应方程

%******************************************************
%二次拟合

k2 = 0;
for y2 = x1
k2 = k2 + 1;
if k2 > k
else
      ze1(k2) = exp(-(k2-1)*afor);
end
end
%ze1

sizeze1 = size(ze1,2);
z4 = ones(1,sizeze1)';
G=[ze1' z4];
X1 = x1';
au20=inv(G'*G)*G'*X1;
au2 = au20';
%z4,X1,G,au20

Aval = au2(1);
Bval = au2(2);
%Aval,Bval
%输出预测的  A,B的值

strcat(x1t1,'=',num2str(Aval),estr,leftbra,num2str(afor1),tstr,rightbra,'+',leftbra,num2str(Bval),rightbra)
%输出时间响应方程

nfinal = sizexd2-1 + 1;
%决定预测的步骤数5  这个步骤可以通过函数传入

%nfinal = sizexd2 - 1 + 1;
%预测的步骤数 1

for  k3=1:nfinal
x3fcast(k3) = constant1*exp(afor1*k3)+ua;
end
%x3fcast
%一次拟合累加值

for  k31=nfinal:-1:0
if k31>1
      x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1);
else
      if k31>0
         x31fcast(k31+1) = x3fcast(k31)-x(1);
      else
         x31fcast(k31+1) = x(1);
      end
end

end
x31fcast
%一次拟合预测值

for  k4=1:nfinal
x4fcast(k4) = Aval*exp(afor1*k4)+Bval;
end
%x4fcast

for  k41=nfinal:-1:0
if k41>1
      x41fcast(k41+1) = x4fcast(k41)-x4fcast(k41-1);
else
      if k41>0
         x41fcast(k41+1) = x4fcast(k41)-x(1);
      else
         x41fcast(k41+1) = x(1);
      end
end

end
x41fcast,x
%二次拟合预测值

%***精度检验p C************//////////////////////////////////
k5 = 0;
for y5 = x
k5 = k5 + 1;
if k5 > sizexd2
else
      err1(k5) = x(k5) - x41fcast(k5);
end
end
%err1
%绝对误差

xavg = mean(x);
%xavg
%x平均值

err1avg = mean(err1);
%err1avg
%err1平均值

k5 = 0;
s1total = 0 ;
for y5 = x
k5 = k5 + 1;
if k5 > sizexd2
else
      s1total = s1total + (x(k5) - xavg)^2;
end
end
s1suqare = s1total ./ sizexd2;
s1sqrt = sqrt(s1suqare);
%s1suqare,s1sqrt
%s1suqare  残差数列x的方差  s1sqrt 为x方差的平方根S1

k5 = 0;
s2total = 0 ;
for y5 = x
k5 = k5 + 1;
if k5 > sizexd2
else
      s2total = s2total + (err1(k5) - err1avg)^2;
end
end
s2suqare = s2total ./ sizexd2;
%s2suqare 残差数列err1的方差S2

Cval = sqrt(s2suqare ./ s1suqare);
Cval
%nnn = 0.6745 * s1sqrt
%Cval  C检验值

k5 = 0;
pnum = 0 ;
for y5 = x
k5 = k5 + 1;
if abs( err1(k5) - err1avg ) < 0.6745 * s1sqrt
      pnum = pnum + 1;
      %ppp = abs( err1(k5) - err1avg )
else
end
end
pval = pnum ./ sizexd2;
pval
%p检验值

%arr1 = x41fcast(1:6)

灰色预测MATLAB程序.txt (3.86 KB, 下载次数: 170)

zan

模型, 测试, function, matlab, 程序

转播0 淘帖0 分享0 收藏3 支持6 反对0 微信

点评

zgtjdxlhz matlab在数学建模中的应用书中有更简便的代码发表于 2013-1-27 23:09

我也说一句

发表

使用道具举报

kelimasa

7 主题	5 听众	188 积分

升级 44%

TA的每日心情

	开心 2012-9-10 21:57

签到天数: 53 天

[LV.5]常住居民I

3^#

发表于 2011-12-16 13:11 |只看该作者 |招呼Ta 关注Ta

厚积薄发发表于 2011-12-15 09:38
9 }- j2 j1 h9 n( \- X; V不错，好东西

嘿嘿，大家一起加油啊~

我也说一句

发表

使用道具举报

mesproc

0 主题	4 听众	23 积分

升级 18.95%

TA的每日心情

	开心 2012-1-28 10:38

签到天数: 1 天

[LV.1]初来乍到

自我介绍: 好孩子一枚

4^#

发表于 2012-1-19 18:13 |只看该作者 |招呼Ta 关注Ta

真心不错是好东西哇~~收下了

我也说一句

发表

使用道具举报

lqg0920

0 主题	5 听众	7 积分

升级 2.11%

该用户从未签到

5^#

发表于 2012-1-30 12:27 |只看该作者 |招呼Ta 关注Ta

不错，好东西

我也说一句

发表

使用道具举报

狗王张董

2 主题	8 听众	3694 积分

升级 56.47%

TA的每日心情

	郁闷 2016-3-28 12:46

签到天数: 1182 天

[LV.10]以坛为家III

6^#

发表于 2012-2-27 15:12 |只看该作者 |招呼Ta 关注Ta

就喜欢你这样的

我也说一句

发表

使用道具举报

梓爱

5 主题	5 听众	58 积分

升级 55.79%

TA的每日心情

	无聊 2012-8-28 02:10

签到天数: 8 天

[LV.3]偶尔看看II

自我介绍: 喜爱建模

7^#

发表于 2012-7-11 09:49 |只看该作者 |招呼Ta 关注Ta

请教下楼主，最后面的c检验值和p检验值是什么意思，标准是什么，什么值算好呀？？

我也说一句

发表

使用道具举报

梓爱

5 主题	5 听众	58 积分

升级 55.79%

TA的每日心情

	无聊 2012-8-28 02:10

签到天数: 8 天

[LV.3]偶尔看看II

自我介绍: 喜爱建模

8^#

发表于 2012-7-11 10:06 |只看该作者 |招呼Ta 关注Ta

梓爱发表于 2012-7-11 09:49
5 n6 M( H0 j, D% J% t, Z" k请教下楼主，最后面的c检验值和p检验值是什么意思，标准是什么，什么值算好呀？？

此问题已解决，详见如下：
if p>0.95 & c<0.35
disp('The model is good,and the forecast is:'),
disp(Hatx0(length(x0)+T))
elseif p>0.85 & c<0.5
disp('The model is eligibility,and the forecast is:'),
disp(Hatx0(length(x0)+T))
elseif p>0.70 & c<0.65
disp('The model is not good,and the forecast is:'),
disp(Hatx0(length(x0)+T))
else p<=0.70 & c>0.65
disp('The model is bad,and try again')