y空间兑换代码_4行关键代码实现灰色模型GM(1, 1)

最新推荐文章于 2023-02-11 07:49:32 发布

weixin_39869733

最新推荐文章于 2023-02-11 07:49:32 发布

阅读量260

点赞数

文章标签： y空间兑换代码投入产出模型matlab代码指数增长模型matlab代码矩阵求逆c语言实现

4行关键代码实现灰色模型GM(1, 1)

1、灰色模型GM(1, 1)

先抄书（嫌烦可略过，建议略过）：

灰色GM(1, 1)模型：

其白化方程通常写为:

这里

在邓和刘的书中通常称为紧邻均值序列，有时也称为背景值，它的公式为:

这里

即是原始数据的累加值，表示为

参数估计式：

其中：

离散响应式（通过求解白化方程得到，取初值为

最后计算还原式：

2、重新梳理计算步骤（重点）

许多人通常在看书的时候就会觉得，一会又是矩阵，一会又是方程，感觉有点麻烦。在实操过程中稍遇到些问题，一想直背后太多知识点就放弃，十分可惜。

下面，我们单纯以实现这个模型的计算为目的，重新梳理其整个建模步骤。

数据准备：假设我们有原始数据

第一步：计算其累加值：

第二步：再算背景值：

第三步：构造矩阵

(在具体实操中，我们并不需要专门构造矩阵

)：

第四步：算矩阵的转置、乘积、求逆再算矩阵乘向量（后面马上就来具体实操，并不难）：

第五步：算时间响应式和还原值（取

）：

3、MATLAB代码手把手实现（以下高能）

在上一节里我们将模型的主要步骤分成了5个步骤，实际上真正的计算过程就是按这个顺序逐步执行的。接下来，我们就来一步步写出各步的代码。

在以下内容中，你将看到熟悉MATLAB编程特性的情况下，实现这一模型是何其简单。

数据准备：

x0 = exp(0.1*[1:10])' + rand(10,1);

注：这里是生成了一组10*1的矩阵，以指数函数为真实值，加入0到1之间的噪声。这里我们采用列向量方便后面计算。

另外，我们设n=8,p=2

n = 8;
p = 2;

第一步：计算其累加值：

x1 = cumsum(x0);

注： cumsum()是MATLAB内置函数，直接求出某一组列向量的累加值，返回一个和原数据同样的矩阵。那么也就是说，这段代码直接实现了累加的所有计算，不需要任何循环。

第二步：计算背景值：

z1 = 0.5 * ( x1(1:n-1) + x1(2:n) );

注：如果上述操作是小高能，那么这一步就是高能操作。即是利用MATLAB中对矩阵元素的取出方式，简单实现了背景值的计算。同样，不需要任何循环。

第三步：构造矩阵:

B = [ -z1 ones(n-1,1)];
Y = x0(2:n);

注：如果上面是高能，那么这里就算是小超高能。利用MATLAB分块矩阵的记法，轻松实现两个矩阵的构造，一句循环都不需要，同时代码简洁易懂。

第四步：计算矩阵乘积、求逆、、、、（不想打了，直接看代码）：

val_a_b = pinv(B)*Y;
a = val_a_b(1);
b = val_a_b(2);

注：不明真相的群众需要补充一些知识，即广义逆矩阵（也称为伪逆矩阵）。见百科：

https://zh.wikipedia.org/wiki/%E5%B9%BF%E4%B9%89%E9%80%86%E9%98%B5zh.wikipedia.org

这里：函数pinv()即是矩阵广义逆的函数，当矩阵行大于列时，这个函数计算出的值即是：

所以我们所看到的这一大坨矩阵计算，一个函数就可以搞定。

第五步: 计算响应式和还原式

这里你以为我要用循环了？（图样图桑破！）看代码：

k = [1:10]';
hat_x1 = ( x0(1) - b/a ) * exp(-a*k) + b/a;
hat_x0 = [ x0(1);hat_x1(2:end)-hat_x1(1:end-1) ];

注：这里属于核能。

这段代码有几个技巧需要说明：

(1) 在计算响应式的时候，我们预先知道需要算出 $k = 1, 2,...,10 $的值，通常的思路是写一个FOR循环让$k$取遍1到10. 不过MATLAB提供了十分灵活的矩阵计算方式。用第二行代码算出的值直接就是响应式所有值构成的列向量，因此，不需要循环。

(2) 在计算还原式的时候，采用了2个技巧。第1个仍然是分块矩阵的运用，这一点不用多说。第2个技巧是MATLAB对下标的操作。这里 hat_x1(2:end) 指的是这一向量中第2个到最后一个元素，同理hat_x1(1:end-1)即是第1到倒数第二个元素。这两组相减，刚好每个元素的值就是

，

所以这里，我们仍然不需要任何循环。

4、福利完整代码：

x0 = exp(0.1*[1:10])' + rand(10,1);
n = 8;
p = 2;
x1 = cumsum(x0);
z1 = 0.5 * ( x1(1:n-1) + x1(2:n) ); 
B = [ -z1 ones(n-1,1)];
Y = x0(2:n);
val_a_b = pinv(B)*Y;
a = val_a_b(1);
b = val_a_b(2);
k = [1:10]';
hat_x1 = ( x0(1) - b/a ) * exp(-a*k) + b/a;
hat_x0 = [ x0(1);hat_x1(2:end)-hat_x1(1:end-1) ];

至于模型的评估，计算误差这些问题，自己简单写一写即可。

提示：仍然全部使用MATLAB自带函数，如mean(),abs() 以及向量的表达式。

5、进一步讨论（核能继续）

上述代码其实仍有简化的空间，有几个地方其实是完全没必要的。

（1）矩阵B和Y并没有必要显式构造出来，同时背景值也只用来做了一下矩阵构造，所以，这三行可以不要。

（2）上述三行代码只是用来算参数，因此把这三行合并到这一行来即可。

即：

z1 = 0.5 * ( x1(1:n-1) + x1(2:n) );  
B = [ -z1 ones(n-1,1)]; 
Y = x0(2:n); 
val_a_b = pinv(B)*Y;

缩减到一行：

val_a_b = pinv([ -0.5 * ( x1(1:n-1) + x1(2:n) ) ones(n-1,1)] )*x0(2:n);

（3）单独用变量a,b,k来存储一些值，只是为了让代码更好懂，其实可以全删了。直接用val_a_b(1), val_a_b(2) 分别代替a,b; 后面的k直接写成[1:n+p]'即可，这里就不拆了。

直接上最终简洁版：

x0 = exp(0.1*[1:10])' + rand(10,1);
n = 8;
p = 2;
x1 = cumsum(x0);
val_a_b = pinv([ -0.5 * ( x1(1:n-1) + x1(2:n) ) ones(n-1,1)] )*x0(2:n);
hat_x1 = ( x0(1) - val_a_b(2)/val_a_b(1) ) * exp(-val_a_b(1)*[1:n+p]') + val_a_b(2)/val_a_b(1);
hat_x0 = [ x0(1);hat_x1(2:end)-hat_x1(1:end-1) ];

仔细看看，前3行是必要的数据准备和参数设置。真正的计算过程4行代码即可。

6、小结及后记

读完本文，如果认真读的话，可能你花得不止10分钟。

读完本文，相信你已经能学会灰色模型的实现，换上你自己的数据试一试，或许会有更多发现。

通过本文，相信你已经能体会到一些MATLAB编程的乐趣，尤其是亲自动手实现过这个模型或者类似模型的朋友。

更多关于灰色预测模型的最新成果，可以参考以下链接：

https://www.researchgate.net/profile/Sifeng_Liuwww.researchgate.netBo Zeng | Chongqing Technology and Business University, Chongqingwww.researchgate.net