GM(1,1)程序设计
GM(1,1)是灰色模型中较为常见的模型,下面是程序,x0是数据,可更改。
之前编辑忘了说了,一般就是给定一组数据,自己根据这些数据拟合一个灰色模型,底下的代码可以得到该模型对应的公式。
代码:
% GM(1,1) % 程序有详尽注释 clc; clear all; x0=[92.810 97.660 98.800 99.281 99.537 99.537 99.817 100.000]; n=length(x0); % 做级比判断,看看是否适合用GM(1,1)建模 lamda=x0(1:n-1)./x0(2:n); range=minmax(lamda); % 判定是否适合用一阶灰色模型建模 if range(1,1)<exp(-(2/(n+2)))| range(1,2)>exp(2/(n+2)) error('级比没有落入灰色模型的范围内,不适合用GM(1,1)建模'); else % 空行输出 disp(' '); disp('可用GM(1,1)建模'); end %做AGO累加处理 x1=cumsum(x0); for i=2:n %计算紧邻均值,也就是白化背景值 z(i)=0.5*(x1(i)+x1(i-1)); end B=[-z(2:n)',ones(n-1,1)]; Y=x0(2:n)'; % 矩阵做除法,计算发展系数和灰色作用量 % 注意:这里是右除不是左除 u=B\Y; % 在MATLAB中,用大写字母D表示导数,dsolve函数用来求解符号常微分方程 x=dsolve('Dx+a*x=b','x(0)=x0'); % subs函数的作用是替换元素,这里是把常量a,b,x0换成具体u(1),u(2),x(1)数值 x=subs(x,{'a','b','x0'},{u(1),u(2),x1(1)}); disp('函数:'); vpa(x,6) forecast1=subs(x,'t',[0:n-1]); % digits和vpa函数用来控制计算的有效数字的位数 digits(6) % y值是AGO形式的(还是累加的) y=vpa(x); % 把AGO输出值进行累减 % diff用于对符号表达时进行求导 % 但是如果输入的是向量,则表示计算原向量相邻元素的差 exchange=diff(forecast1); % 输出灰色模型预测的值 disp('输出预测模型预测的值:'); forecast=[x0(1),exchange] % 计算残差 epsilon=x0-forecast; % 计算相对误差 delta=abs(epsilon./x0); % 检验模型的误差 % 检验方法一:相对误差Q检验法 disp('相对误差Q检验值:'); Q=mean(delta) % 检验方法二:方差比C检验法 % 计算标准差函数为std(x,a) % 如果后面一个参数a取0表示的是除以n-1,如果是1就是最后除以n disp('方差比C检验值:'); C=std(epsilon,1)/std(x0,1) % 检验方法三:小误差概率P检验法 S1=std(x0,1); S1_new=S1*0.6745; temp_P=find(abs(epsilon-mean(epsilon))<S1_new); disp('小误差概率P检验值:'); P=length(temp_P)/n %绘制原始数列和灰色模型预测得出的数列差异折线图 plot(1:n,x0,'ro','markersize',11); hold on plot(1:n,forecast,'k-','LineWidth',2.5); grid on; axis tight; xlabel('x'); ylabel('y'); title('保有量比例与时间序列的时间关系'); legend('原始数列','模型数列');
运行结果,程序输出的图形如下所示:
版权声明:本文为博主原创文章,未经博主允许不得转载。
posted on 2015-09-05 23:29 Tob__yuhong 阅读(528) 评论(0) 编辑 收藏 举报