Matlab-灰色预测

This is my blog.
灰色预测貌似是数模中的很重要的一部分
然后我知道就一个模版之后
嗯，敲了一遍，大概就是之后改数据就行
突然觉得编程很水
主要还是靠思路取胜的
果然我还是适合打酱油
不过，干这事估计我的英语阅读能提高吧！
哎，为什么我一个女孩子
英语可以那么弱啊
唔，不开心

灰色不是黑也不是白。。。

基本理论

灰色关联度的本质是依据曲线态势相近程度来分辨数列的相关度。
灰色模型是基于常微分理论的，其记号为\(GM(M,N)\)。其中，N表示变量的个数，M表示常微分方程的阶数。

适用范围

灰色模型只适用于连续、平滑的事件趋势中。且只有符合一阶线性常微分方程的数列或矩阵才能用灰色模型来预测，否则会有很大的误差。

经典灰色模型 GM(1,1)

求解步骤

数据的累加（AGO）
不仅可以降噪，还可以加强数据规律的显露
\(x^{(1)}(k)=\sum_{i=1}^kx^{(0)}(i),k=1,2,3,\cdot\cdot\cdot,n\)
求出数列\(x^{(1)}\)的导数
\(\left. \frac{dC}{dt} \right| _k=x^{(1)}(k)- x^{(1)}(k-1)=x^{(0)}(k)=d(k)\)
定义数列\(x^{(1)}(k) \)的紧邻均值（背景值）
\(z^{(1)}(k)=\frac{x^{(1)}(k)+x^{(1)}(k-1)}{2},k=2,3,4\cdot\cdot\cdot\)
称新数列\(z^{(1)}=(z^{(1)}(2),z^{(1)}(3),z^{(1)}(4),\cdot\cdot\cdot,z^{(1)}(n)）\)为\(x^{(1)}\)的紧邻均值数列。
定义\(GM(1,1)\)灰微分方程
\(x^{(0)}(k)+az^{(0)}(k)=b\)
a称为“发展系数”，b称为“灰作用量”。
令\(Y=(x^{(1)}(2),x^{(1)}(3),x^{(1)}(4),\cdot\cdot\cdot ,x^{(1)}(n))^T,u=(a,b)^T,B=\)

即可表示为\(Y=Bu\)
\(\breve{u}=(\breve{a},\breve{b})^T=(B^TB)^{-1}B^TY\)
这里所得到的\(\breve{u}\)是近似值

白化GM(1,1)模型
离散到连续的过程
白化模型
\(\frac{dx^{(1)}}{dk}+ax^{(1)}=b\)
求解模型
将5中的数据变形可得
连续函数表达式 \(X=(x^{(0)}(1)-\frac ba)e^{-ak},k=0,1,2,3,\cdot\cdot\cdot\)
具体数学表达式 \(\breve{x}^{(1)}(k+1)=x^{(0)}(1)-\frac ba)e^{-ak}+\frac ba\)
将数据累减(IAGO)从而还原为
\(\breve{x}^{(0)}(k+1)=\breve{x}^{(1)}(k+1)-\breve{x}^{(1)}(k),k=1,2,3,\cdot\cdot\cdot,n\)
模型检验
1）相对残差Q检验（越小越好）
\(\varepsilon_k=x^{(k)}-\breve{x}^{(k)}\)
\(Q=\frac 1n\sum_{k=1}^n\frac {\varepsilon_k}{x^{(0)}(k)}\)
2)方差比C检验（越小越好）
\(\varepsilon^{(0)}=(\varepsilon_1,\varepsilon_2,\varepsilon_3,\cdot\cdot\cdot,\varepsilon_n)\)
\(\overline{x}=\frac 1n\sum_{k=1}^nx^{(0)}(k)\)
\(\overline{\varepsilon}=\frac 1n\sum_{k=1}^n\varepsilon^{(0)}(k)\)
\(S1^2=\frac 1n\sum_{(k=1)}^n(x^{(0)}(k)-\overline{x})\)
\(S2^2=\frac 1n\sum_{(k=1)}^n(\varepsilon^{(0)}(k)-\overline{\varepsilon})\)
3)小误差概率p的检验（置信水平取0.5时）（p越大越好）
\(p=P(\left|\varepsilon_k-\overline{\varepsilon} \right| <0.6745S1)\)
模型的预测
通过计算数列的“级比”来判断是都可用GM(1,1)来预测
\(\lambda(k)=\frac {x^{(0)}(k-1)}{x^{(0)}(k)},k=2,3,4,\cdot\cdot\cdot\)
如果\(\lambda(k)\)能够落入到\(\left[e^{-\frac{2}{n+1}},e^{\frac{2}{n+1}}\right]\)区间内，则可用GM(1,1)来预测。

Matlab

x0 = [92.810 97.660 98.800 99.281 99.537 99.537 99.817 100.00];
n = length(x0);
% 做级比判断看是否适合用GM(1,1)建模
lamda = x0(1:n-1) ./ x0(2:n);
range = minmax(lamda);% minmax 用于获取数组中每一行的最小值和最大值
if range(1,1) < exp(-e/(n+2))) | range(1,2) > exp(2/(n+2))
    error('级比没有落入灰色模型的范围内');
else
    disp('          ');
    disp('0');
end
% 做AGO累加处理
x1 = cumsum(x0);% cumsum 求各行的累加值
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 = dslove('Dx + a * x = b','x(0) = x0');
% subs函数的作用是替换元素，这里是把常量a,b,x0换成具体u(1),u(2),x1(1)数值
x = subs(x,{'a','b','x0'},{u(1),u(2),x1(1)});
forecast1 = subs(x,'t',[0:n-1]);
digits(6);
% R = vpa(A, d) 用d个位数代替当前设置的位数。
% vpa是在运算时就控制了精度，但是digits是对结果控制精度
y = vpa(x);
% 累减
% diff用于对符号表达式进行求导
% diff用于对向量时，是计算原向量相邻元素的差
exchange = diff(forecast1);
% 输出灰色预测的值
forcast = [x0(1),exchange]
% 计算残差
epsilon = x0 - forecast;
% 计算相对误差
delta = abs(epsilon ./ x0);
% 相对误差Q的检验
Q = mean(delta)
% 方差比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);
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('原始数列','模型数列');

也可以直接用6中的具体数学公式：\(\breve{x}^{(1)}(k+1)=x^{(0)}(1)-\frac ba)e^{-ak}+\frac ba\)来编写程序；

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 = inv(B' * B) * B' * Y;
% 利用GM(1,1)具体表达式计算原始数列的AGO
forecast1 = (x1(1) - u(2) ./ u(1)) .* exp(-u(1) .* ([0:n-1])) + u(2) ./ u(1);
% 累减
% diff用于对符号表达式进行求导
% diff用于对向量时，是计算原向量相邻元素的差
exchange = diff(forecast1);
% 输出灰色预测的值
forcast = [x0(1),exchange]
% 计算残差
epsilon = x0 - forecast;
% 计算相对误差
delta = abs(epsilon ./ x0);
% 相对误差Q的检验
Q = mean(delta)
% 方差比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);
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('原始数列','模型数列');

灰色Verhulst模型

曲线为S型

步骤

IAGO
\(x^{(0)}(k)=x^{(1)}(k)-x^{(1)}(k-1),k=2,3,4，\cdot\cdot\cdot,n\)
对原始数列\(x^{(1)}\)作紧邻均值
\(z^{(1)}(k)=\frac{x^{(1)}(k)+x^{(1)}(k-1)}{2},k=2,3,4\cdot\cdot\cdot\)
定义数据矩阵
B =
Y =
计算参数
\(\breve{u}=(\breve{a},\breve{b})^T=(B^TB)^{-1}B^TY\)
定义Verhulst模型白化矩阵
\(\frac{x^{(1)}}{dk}+ax^{(1)}=b(x^{(1)})^2\)
求解Verhulst模型白化矩阵
\(\breve{x}^{(1)}(k+1)=\frac{\breve{a}x^{(0)}(1)}{\breve{b}x^{(0)}(1)+\breve{a}-\breve{b}x^{(0)}(1))e^{\breve{a}k}}\)
检验模型

Matlab

x1 = [0.025 0.023 0.029 0.044 0.084 0.164 0.332 0.521 0.97 1.6 2.45 3.11 3.57 3.76 3.96 4 4.46 4.4 4.49 4.76 5.01];
n = length(x1);
x0 = diff(x1);
x0 = [x1(1),x0];
for i = 2:n
    z1(i) = 0.5 * (x1(i) + x1(i-1));
end
z1;
B = [-z1(2:end)',z1(2:end)' .^ 2];
Y = x0(2:end)';
abvalue =B \ Y;
x = dsolve('Dx + a * x = b * x ^ 2','x(0) = x0');
x = subs(x,{'a','b','x0'},{abvalue(1),abvalue(2),x1(1)});
forecast1 = subs(x,'t',0:n-1);
digits(6);
y = vpa(x);
forecast

结束语

一般用灰色模型的时候，都用经典型的，且可直接用上述代码，只需把数据更换一下即可。（解法中一般较多使用公式法的那种）
灰色模型也是一种拟合。

后记

早上踩着白雪
估计是没有人踩出路来
也估计是自己走路的时候总是爱想别的事的缘故
（突然想到，中午上课的时候，如果看到我们班的人，就很开心，
看着他鞋子走路，既不会迷路，也不会迟到，还可以乱想事情
hhhhh，感觉自己是个猥琐狂- -
直接单膝跪地摸到了软软的雪
看了琅琊榜2-36集，超级燃，大爱孙淳老师
重温了琅琊榜1，其实我觉得2拍的不烂
都超级好看呢！！！
感觉那才是大爱大情大义
bgm超级入境
找到了《月亮与六便士》
不过不是英版的，打算先看看
其实本来打算翻译这本的
偶然看到的片段觉得很美
嘻嘻嘻
打算上午出门，下午就窝宿舍了
最近偷懒一会咯

转载请注明出处，谢谢。

愿我是你的小太阳