灰色马尔科夫预测 matlab代码（数据量太少，有局限性）

首先声明，本人只是个刚学matlab不到一周的纯小白，写灰色马尔科夫是因为数学建模培训练题的时候要用到，但是在网上找不到现成的能用的代码（啊没错，我就是那种白嫖党），而且找到的基本都是“付费观看”。我们组练的那道题主要用的模型并不是灰色马尔科夫，灰马在我们的模型里就相当于一个数据处理的环节，最后权重占得也不大（而且那题的优秀论文证明是我们思路偏了，原本根本用不上灰马），所以具体代码会有局限性，这里只是放出来分享一下。

我们当时的题是06年的A题，出版社资源配置。训练的时候不允许直接搜赛题内容，这也是我们第一次模拟练题，真的很菜，最后的思路也有点天马行空，反正和优秀例文的思路差的有点大，不过这里主要分享的是灰马，后面会提一下赛题内容，主要是为了讲解一下代码（对不起我就是有点啰里啰唆的）。

灰马的理论参考有csdn上的文章（灰色马尔科夫模型matlab实现 ），以及知网上找的论文，主要参考的有两篇

and

大家可以自己上知网找（ps：我们当时练题其实还是比较随便的，搞得差不多了之后就是队长在写论文，我比较闲，虽然前面没学过马尔科夫预测，但因为专业课程之前有信息论，所以有点马尔科夫链的基础，就想着随便搞搞灰马当作改进。啰嗦这些就是想说看灰马的除了要大概懂一懂灰色预测，没看过马尔可夫的最好把马尔可夫也看看）

马尔可夫部分的代码有参考大佬的文章马尔科夫预测MATLAB

概率转移矩阵那里我用的最传统的算法，有点笨，但我太菜了，当时也是只想随便整整能用就行，没搞明白常用的算法。

灰色预测的matlab代码参考了05年的长江水质问题的代码（我们用的是前期老师培训时给的matlab代码，好像是某个优秀论文的代码，csdn上也有，和我们参考的那份是一样的）

灰马在本菜比的理解就是先灰色预测，再用马尔可夫处理，给灰色预测出来的值加上波动。灰色预测用的累加或者累减处理，适合指数递增类型，不适用于具有波动性的数据。如果数据量大好像可以用时间序列之类的方法，但用灰色就是因为数据量少，所以用马尔可夫，分状态，算转移概率矩阵，最后得到改进的波动值，用波动值加上原本灰色预测出来的值，这样最后的值就具有波动性了。

我们当时就是大部分数据用了灰色预测，但有一份数据波动性大，我们这题因为数据量太少了有点不太会搞，灰马的话马尔可夫部分的状态最好三个以上，但是这题本身的数据就只有5个（72个科目，一个科目算一份，一份里面用前五年预测第六年的），灰色预测的话第一个预测值好像就是实际值（论文里也看到有第一个预测值和实际值不一样的，但我不知道怎么弄的，老师给的那份灰色预测代码从原理上出来的第一个预测值就是实际值，用spsspro弄也是这样）。因为后面马尔可夫用的数据是实际值和预测值的差值或相对变化率，所以第一个值没法用，这样能用的就只有4年的数据，感觉应该是能用的最低极限了。最后我是分了3个状态（我了解到的应该是最少三个，不然后面状态转移矩阵没法搞），参考了上面提到的第二篇论文。但论文里面数据分在哪个状态的判断我用matlab不会处理（果然还是太菜了），所以只能改了一下，但改完之后又会出来其他情况，于是就有了下面代码里面的一堆判断。

改完之后出来的波动值可能会跟按论文的方法出来的结果有点差别，但因为我们这题的数据值不是很大，影响比较小，加上真的是数据量太小了，压着能用的边儿用的（论文里的是5个数据，我们的只有4个数据，但要分3个状态，真的要死了），所以最后就这样了。如果数据量再大一点，比如有六个这样，改一改代码，原本论文里的方法就能继续用了，最后结果也会正常不少，我这里这么搞纯粹是因为我们的数据太少了。

clc,clear

format rat;%使用分数来表示值

load('ma.mat')

for n=1:72 %因为我的数据有72组，所以重复处理72次

a=ma(n,:)

b=sort(a);f=zeros(3,3);

d=find(b==a(4))%2005年所处的状态区间

if(d==1)

    b2=b(1:2);b1=b(3);b3=b(4);

     d=2;p=[b(3),b(1),b(2),b(4)];

elseif(d==4)

    b1=b(1);b3=b(2);b2=b(3:4);

    d=2;p=[b(1),b(3),b(4),b(2)];

else

    b1=b(1);b2=b(2:3);b3=b(4);

    p=b;

end



for i=1:3   

        if(ismember(a(i),b1))

            c1=1;

        elseif(ismember(a(i),b2))

            c1=2;

        else

        c1=find(p==a(i));

        end

        

         if(ismember(a(i+1),b1))

            c2=1;

        elseif(ismember(a(i+1),b2))

            c2=2;

        else

        c2=find(p==a(i+1));

        end

        

        if(ismember(a(i),b2))

            if(ismember(a(i+1),b2))

                f(2,2)=f(2,2)+1;

            else

                if(c1==4)c1=3;end

                if(c2==4)c2=3;end

                f(2,c2)=f(2,c2)+1;

            end

        else

            if(ismember(a(i+1),b2))

                if(c1==4)c1=3;end

                f(c1,2)=f(c1,2)+1;

            else

                if(c1==4)c1=3;end

                if(c2==4)c2=3;end

                f(c1,c2)=f(c1,c2)+1;

            end

        end

end



ni=sum(f,2);%按照行求和

phat=f./repmat(ni,1,size(f,2));%矩阵除法，repamt生成了一个矩阵，其维数为一维，即数组形式，元素分别为行和的副本

disp(phat);%size(f,2)求列数

format %恢复到短小数的表示形式



if(d==3)%找一步转移矩阵中2005年所处那行的最大概率是哪个

    m=max(phat(2,:));%此行中概率最大值

    s=find(phat(2,:)==m)%06年将会处于的状态

elseif (d==4)

    m=max(phat(3,:));

    s=find(phat(3,:)==m)

else

    m=max(phat(d,:));

    s=find(phat(d,:)==m)

end    

if(s==3)

    n=find(b==p(s+1));%状态3在p中对应的值是第4个

else

    n=find(b==p(s));

end

if(n==4)

    med=(b(n)+b(n-1))/2

else

med=(b(n)+b(n+1))/2

end

end

 

有问题球球大家轻喷，我真的只是个小白，放出来也只是单纯交流分享，不是教程更不是参考什么的。

原理以及代码说明后续有空再更