【Matlab编程】马氏链随机模拟

本文是利用蒙特卡罗算法对马氏链过程的模拟。假设有10个状态，从每个状态到与之相邻状态的概率是相同的，仿真次数为1000，及进行了1000次状态转移。我们以动画的形式再现了状态转移的过程，并记录了到达每个状态的次数，具体实现如下：

close all;clc;clear;
figure;
s=1;
n=1000;
r=1; % 圆圈的半径
title('等概率情况的计算机模拟')
set(gcf,'doublebuffer','on'); % 设置图形渲染效果
xlabel('Please press "space" key and see the result!',...
   'fontsize',14,'color','r'); % 添加标注文字
hold on;axis equal; % 设置坐标轴属性
axis([-16,16,-16,16]); % 设置坐标轴范围
fill(r*sin(0:.1:2*pi)-7,r*cos(0:.1:2*pi)+13,'w'); % 画出固定点P1
hold on
fill(r*sin(0:.1:2*pi)-11,r*cos(0:.1:2*pi)+9,'w');hold on % 画出固定点P2
fill(r*sin(0:.1:2*pi)-3,r*cos(0:.1:2*pi)+9,'w'); hold on% 画出固定点P3
fill(r*sin(0:.1:2*pi)+5,r*cos(0:.1:2*pi)+9,'w');hold on % 画出固定点P4
fill(r*sin(0:.1:2*pi)+9,r*cos(0:.1:2*pi)+5,'w');hold on % 画出固定点P5
fill(r*sin(0:.1:2*pi)-15,r*cos(0:.1:2*pi)-3,'w');hold on % 画出固定点P6
fill(r*sin(0:.1:2*pi)+1,r*cos(0:.1:2*pi)-3,'w'); hold on% 画出固定点P7
fill(r*sin(0:.1:2*pi)+13,r*cos(0:.1:2*pi)-3,'w');hold on % 画出固定点P8
fill(r*sin(0:.1:2*pi)-7,r*cos(0:.1:2*pi)-11,'w');hold on % 画出固定点P9
fill(r*sin(0:.1:2*pi)+5,r*cos(0:.1:2*pi)-15,'w');hold on % 画出固定点P10
text(-15.4,-3,'6','FontSize',18);hold on
text(-11.4,9,'2','FontSize',18);hold on
text(-7.4,13,'1','FontSize',18);hold on
text(-7.4,-11,'9','FontSize',18);hold on
text(-3.4,9,'3','FontSize',18);hold on
text(0.6,-3,'7','FontSize',18);hold on
text(4.6,9,'4','FontSize',18);hold on
text(4.1,-15,'10','FontSize',18);hold on
text(8.6,5,'5','FontSize',18);hold on
text(12.6,-3,'8','FontSize',18);hold on
hold on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plot([8.5,6],[6,8.3],'r-')
hold on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x45=fliplr(8.3:-0.1:5.8);
y45=fliplr(linspace(5.9,8.2,length(x45)));
x54=8.3:-0.1:5.8;
y54=linspace(5.9,8.2,length(x54));
x85=fliplr(9.4:0.05:12.4);
y85=fliplr(linspace(4.1,-2.2,length(x85)));
x58=9.4:0.05:12.4;
y58=linspace(4.1,-2.2,length(x58));
x80=fliplr(5.8:0.1:12.2);
y80=fliplr(linspace(-14.4,-3.8,length(x80)));
x08=5.8:0.1:12.2;
y08=linspace(-14.4,-3.8,length(x08));
x87=fliplr(2:0.1:12);
y87=-3*ones(1,length(x87));
x78=2:0.1:12;
y78=-3*ones(1,length(x78));
x79=fliplr(-6.2:0.1:0.4);
y79=fliplr(linspace(-10.4,-3.8,length(x79)));
x97=-6.2:0.1:0.4;
y97=linspace(-10.4,-3.8,length(x97));
x73=fliplr(-2.6:0.06:0.9);
y73=fliplr(linspace(7.9,-2,length(x73)));
x37=-2.6:0.06:0.9;
y37=linspace(7.9,-2,length(x37));
x13=-6.4:0.1:-3.8;
y13=linspace(12.2,9.6,length(x13));
x31=fliplr(-6.4:0.1:-3.8);
y31=linspace(9.6,12.2,length(x31));
x67=-14:.11:0;
y67=-3*ones(1,length(x67));
x76=0:-0.11:-14;
y76=-3*ones(1,length(x76));
x21=-10.1:.1:-7.8;
y21=linspace(9.5,12.4,length(x21));
x12=-7.8:-0.1:-10.1;
y12=fliplr(linspace(9.5,12.4,length(x12)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t6=text(-15.3,-5,'0','FontSize',12,'Color',[1 0 0]);
t2=text(-11.4,7,'0','FontSize',12,'Color',[1 0 0]);
t1=text(-7.2,11,'0','FontSize',12,'Color',[1 0 0]);
t9=text(-7.3,-13,'0','FontSize',12,'Color',[1 0 0]);
t3=text(-4,7,'0','FontSize',12,'Color',[1 0 0]);
t7=text(0.6,-5,'0','FontSize',12,'Color',[1 0 0]);
t4=text(4.7,7,'0','FontSize',12,'Color',[1 0 0]);
t10=text(4.3,-13,'0','FontSize',12,'Color',[1 0 0]);
t5=text(8.3,3,'0','FontSize',12,'Color',[1 0 0]);
t8=text(12.6,-5,'0','FontSize',12,'Color',[1 0 0]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:length(x45)
    plot(x45(i),y45(i),'*')

end

for i=1:length(x54)
    plot(x54(i),y54(i),'*')

end

for i=1:length(x85)
    plot(x85(i),y85(i),'*')

end

for i=1:length(x58)
    plot(x58(i),y58(i),'*')

end

for i=1:length(x80)
    plot(x80(i),y80(i),'.')

end

for i=1:length(x08)
    plot(x08(i),y08(i),'*')
end

for i=1:length(x87)
    plot(x87(i),y87(i),'*')

end

for i=1:length(x78)
    plot(x78(i),y78(i),'*')
end

for i=1:length(x79)
    plot(x79(i),y79(i),'.')

end

for i=1:length(x97)
    plot(x97(i),y97(i),'*')

end

for i=1:length(x73)
    plot(x73(i),y73(i),'*')

end

for i=1:length(x37)
    plot(x37(i),y37(i),'.')

end

for i=1:length(x31)
    plot(x31(i),y31(i),'*')
end

for i=1:length(x13)
    plot(x13(i),y13(i),'*')

end

for i=1:length(x67)
    plot(x67(i),y67(i),'*')

end

for i=1:length(x76)
    plot(x76(i),y76(i),'*')

end

for i=1:length(x21)
    plot(x21(i),y21(i),'*')

end

for i=1:length(x12)
    plot(x12(i),y12(i),'*')

end
plot(x45,y45,'w.')
plot(x85,y85,'w.')
plot(x80,y80,'w.')
plot(x87,y87,'w.')
plot(x79,y79,'w.')
plot(x73,y73,'w.')
plot(x31,y31,'w.')
plot(x21,y21,'w.')
plot(x67,y67,'w.')
  plot(x54,y54,'w.')
plot(x58,y58,'w.')
plot(x08,y08,'w.')
plot(x78,y78,'w.')
plot(x97,y97,'w.')
plot(x37,y37,'w.')
plot(x13,y13,'w.')
plot(x12,y12,'w.')
plot(x76,y76,'w.')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   s=1;
    p1=0;%pn为到达n村庄的次数
 p2=0;p3=0;p4=0;p5=0;p6=0;p7=0;p8=0;p9=0;p10=0;

%plot([-14,0],[-3,-3],'b','linewidth',3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:n
    m=get(gcf,'currentkey'); % 获取键入按键的名称
   if strcmp(m,'space'); % 检查按下的按键是否为空格键
       break;
   end
    if s==1
        possible=round(rand(1));
        if possible==1
        s=3;
        p3=p3+1;
        for i=1:length(x13)
            plot(x13(i),y13(i),'.')
            pause(0.00000000001)
        end
        delete(t3);
        t3=text(-4,7,num2str(p3),'FontSize',12,'Color',[1 0 0]);
        else
            s=2;
            p2=p2+1;
         for i=1:length(x12)
            plot(x12(i),y12(i),'.')
            pause(0.00000000001)
         end
          delete(t2)
          t2=text(-11.4,7,num2str(p2),'FontSize',12,'Color',[1 0 0]);

        end

    elseif s==2
           s=1;
           p1=p1+1;
           for i=1:length(x21)
                plot(x21(i),y21(i),'.')
                pause(0.00000000001)
           end
            delete(t1)
            t1=text(-7.2,11,num2str(p1),'FontSize',12,'Color',[1 0 0]);
    elseif s==3
        possible=round(rand(1));
        if possible==1
        s=7;
        p7=p7+1;
        for i=1:length(x37)
        plot(x37(i),y37(i),'.')
        pause(0.00000000001)
        end
        delete(t7)
        t7=text(0.6,-5,num2str(p7),'FontSize',12,'Color',[1 0 0]);

        else
            s=1;
            p1=p1+1;
            for i=1:length(x31)
                plot(x31(i),y31(i),'.')
                pause(0.00000000001)
            end
            delete(t1)
             t1=text(-7.2,11,num2str(p1),'FontSize',12,'Color',[1 0 0]);
        end

    elseif s==4
           s=5;
           p5=p5+1;
        for i=1:length(x45)
            plot(x45(i),y45(i),'.')
            pause(0.00000000001)
        end
        delete(t5)
        t5=text(8.3,3,num2str(p5),'FontSize',12,'Color',[1 0 0]);

    elseif s==5
        possible=round(rand(1));
        if possible==1
        s=8;
        p8=p8+1;
        for i=1:length(x58)
            plot(x58(i),y58(i),'.')
            pause(0.00000000001)
        end
        delete(t8)
        t8=text(12.6,-5,num2str(p8),'FontSize',12,'Color',[1 0 0]);

        else
            s=4;
            p4=p4+1;
            for i=1:length(x54)
                plot(x54(i),y54(i),'.')
                pause(0.00000000001)
            end
            delete(t4)
            t4=text(4.7,7,num2str(p4),'FontSize',12,'Color',[1 0 0]);

            end

    elseif s==6
         s=7;
           p7=p7+1;
           for i=1:length(x67)
                plot(x67(i),y67(i),'.')
                pause(0.00000000001)
           end
           delete(t7)
           t7=text(0.6,-5,num2str(p7),'FontSize',12,'Color',[1 0 0]);
    elseif s==7
        possible=floor(rand(1)*4);
        if possible==0
            s=6;
            p6=p6+1;
            for i=1:length(x76)
                plot(x76(i),y76(i),'.')
                pause(0.00000000001)
            end
            delete(t6)
            t6=text(-15.3,-5,num2str(p6),'FontSize',12,'Color',[1 0 0]);

        elseif possible==1
            s=3;
            p3=p3+1;
            for i=1:length(x73)
                plot(x73(i),y73(i),'.')
                pause(0.00000000001)
            end
            delete(t3)
            t3=text(-4,7,num2str(p3),'FontSize',12,'Color',[1 0 0]);

        elseif possible==2
            s=9;
            p9=p9+1;
            for i=1:length(x79)
                plot(x79(i),y79(i),'.')
                pause(0.00000000001)
            end
            delete(t9)
            t9=text(-7.3,-13,num2str(p9),'FontSize',12,'Color',[1 0 0]);
        else
            s=8;
            p8=p8+1;
            for i=1:length(x78)
                plot(x78(i),y78(i),'.')
                pause(0.00000000001)
            end
            delete(t8)
            t8=text(12.6,-5,num2str(p8),'FontSize',12,'Color',[1 0 0]);

        end

    elseif s==8
         possible=floor(rand(1)*3);
        if possible==0
            s=7;
            p7=p7+1;
            for i=1:length(x87)
                plot(x87(i),y87(i),'.')
                pause(0.00000000001)
            end
            delete(t7)
            t7=text(0.6,-5,num2str(p7),'FontSize',12,'Color',[1 0 0]);
        elseif possible==1
            s=5;
            p5=p5+1;
            for i=1:length(x85)
                plot(x85(i),y85(i),'.')
                pause(0.00000000001)
            end
            delete(t5)
            t5=text(8.3,3,num2str(p5),'FontSize',12,'Color',[1 0 0]);

        else
            s=10;
            p10=p10+1;
            for i=1:length(x80)
                plot(x80(i),y80(i),'.')
                pause(0.00000000001)
            end
            delete(t10)
            t10=text(4.3,-13,num2str(p10),'FontSize',12,'Color',[1 0 0]);
        end

    elseif s==9
        s=7;
        p7=p7+1;
        for i=1:length(x97)
            plot(x97(i),y97(i),'.')
            pause(0.00000000001)
        end
        delete(t7)
        t7=text(0.6,-5,num2str(p7),'FontSize',12,'Color',[1 0 0]);

    else
        s=8;
        p8=p8+1;
        for i=1:length(x08)
            plot(x08(i),y08(i),'.')
            pause(0.00000000001)
        end
        delete(t8)
        t8=text(12.6,-5,num2str(p8),'FontSize',12,'Color',[1 0 0]);

    end
plot(x45,y45,'w.')
plot(x85,y85,'w.')
plot(x80,y80,'w.')
plot(x87,y87,'w.')
plot(x79,y79,'w.')
plot(x73,y73,'w.')
plot(x31,y31,'w.')
plot(x21,y21,'w.')
plot(x67,y67,'w.')
plot(x54,y54,'w.')
plot(x58,y58,'w.')
plot(x08,y08,'w.')
plot(x78,y78,'w.')
plot(x97,y97,'w.')
plot(x37,y37,'w.')
plot(x13,y13,'w.')
plot(x12,y12,'w.')
plot(x76,y76,'w.')
end

for j=i:n
    if s==1
        possible=round(rand(1));
        if possible==1
        s=3;
        p3=p3+1;
        else
            s=2;
            p2=p2+1;
        end

    elseif s==2
           s=1;
           p1=p1+1;

    elseif s==3
        possible=round(rand(1));
        if possible==1
        s=7;
        p7=p7+1;
        else
            s=1;
            p1=p1+1;
        end

    elseif s==4
           s=5;
           p5=p5+1;

    elseif s==5
        possible=round(rand(1));
        if possible==1
        s=8;
        p8=p8+1;
        else
            s=4;
            p4=p4+1;
        end

    elseif s==6
         s=7;
           p7=p7+1;

    elseif s==7
        possible=floor(rand(1)*4);
        if possible==0
            s=6;
            p6=p6+1;
        elseif possible==1
            s=3;
            p3=p3+1;
        elseif possible==2
            s=9;
            p9=p9+1;
        else
            s=8;
            p8=p8+1;
        end

    elseif s==8
         possible=floor(rand(1)*3);
        if possible==0
            s=7;
            p7=p7+1;
        elseif possible==1
            s=5;
            p5=p5+1;
        else
            s=10;
            p10=p10+1;
        end

    elseif s==9
        s=7;
        p7=p7+1;

    else
        s=8;
        p8=p8+1;
    end

end
delete(t1)
delete(t2)
delete(t3)
delete(t4)
delete(t5)
delete(t6)
delete(t7)
delete(t8)
delete(t9)
delete(t10)
t6=text(-15.3,-5,num2str(p6),'FontSize',12,'Color',[1 0 0]);
t2=text(-11.4,7,num2str(p2),'FontSize',12,'Color',[1 0 0]);
t1=text(-7.2,11,num2str(p1),'FontSize',12,'Color',[1 0 0]);
t9=text(-7.3,-13,num2str(p9),'FontSize',12,'Color',[1 0 0]);
t3=text(-4,7,num2str(p3),'FontSize',12,'Color',[1 0 0]);
t7=text(0.6,-5,num2str(p7),'FontSize',12,'Color',[1 0 0]);
t4=text(4.7,7,num2str(p4),'FontSize',12,'Color',[1 0 0]);
t10=text(4.3,-13,num2str(p10),'FontSize',12,'Color',[1 0 0]);
t5=text(8.3,3,num2str(p5),'FontSize',12,'Color',[1 0 0]);
t8=text(12.6,-5,num2str(p8),'FontSize',12,'Color',[1 0 0]);

仿真过程如下：

最终的结果如下：

【Matlab编程】马氏链随机模拟

时间： 2024-10-24 18:25:23

【Matlab编程】马氏链随机模拟的相关文章

【随机过程】马氏链的理论与仿真

在2014年终总结中,我提到要对这学期学过的数学课中的部分算法进行仿真实现.<数值分析>和<工程优化>这两门数学课里面还有些专门讲算法的,可以用来仿真.在<随机过程>这门课中,几乎全都是公式推导,定理证明,实在难以仿真实现.最后发现,马尔科夫链这一章比较适合仿真,况且先前也写过类似的程序,更重要的是之前有人也问过关于马氏链的Matlab实现问题.关于马氏链的理论原理在这就不作描述,下面直接用程序来实现具体问题的求解. 假设有9个状态,其状态转移图如下所示: 根据状态转移

[数学建模]马氏链模型(一)基本介绍

马氏链是一个随机过程 由第n部的状态i 第m部 转移到 状态 j 的概率 与 n部之前的概率是无关的 ,  我们称为马氏链 如果n部概率与时间无关, 那么成为齐时性的马氏链 给出一个简单的马氏链应用,  a是一条 E属于 0,1 空间的马氏链, 我们的一步转移矩阵可以求出: a1='1110010011111110011110111111001111111110001101101'; a2='111011011010111101110111101111110011011111100111'; a

MATLAB求马氏距离(Mahalanobis distance)

MATLAB求马氏距离(Mahalanobis distance) 作者:凯鲁嘎吉 - 博客园 http://www.cnblogs.com/kailugaji/ 1.马氏距离计算公式 d2(xi, xj)=(xi-xj)TS-1(xi-xj) 其中,S是总体的协方差矩阵,而不是样本的协方差矩阵. 2.matlab中现有的函数 >> x=[155 66;180 71;190 73;160 60;190 68;150 58;170 75] x = 155 66 180 71 190 73 160

随机模拟（MCMC）

http://cos.name/2013/01/lda-math-mcmc-and-gibbs-sampling/ http://blog.csdn.net/lin360580306/article/details/51240398 随机模拟(或者统计模拟)方法有一个很酷的别名是蒙特卡罗方法(Monte Carlo Simulation).这个方法的发展始于20世纪40年代,和原子弹制造的曼哈顿计划密切相关,当时的几个大牛,包括乌拉姆.冯.诺依曼.费米.费曼.Nicholas Metropoli

随机采样和随机模拟：吉布斯采样Gibbs Sampling

http://blog.csdn.net/pipisorry/article/details/51373090 马氏链收敛定理 马氏链定理: 如果一个非周期马氏链具有转移概率矩阵P,且它的任何两个状态是连通的,那么limn→∞Pnij 存在且与i无关,记limn→∞Pnij=π(j), 我们有 limn→∞Pn=???????π(1)π(1)?π(1)?π(2)π(2)?π(2)??????π(j)π(j)?π(j)????????????? π(j)=∑i=0∞π(i)Pij π 是方程 πP

随机模拟与采样方法

随机模拟方法简介 随机模拟方法又称为蒙特卡罗方法(Monte Carlo Method).蒙特卡洛模拟方法的原理是当问题或对象本身具有概率特征时,可以用计算机模拟的方法产生抽样结果,根据抽样计算统计量或者参数的值:随着模拟次数的增多,可以通过对各次统计量或参数的估计值求平均的方法得到稳定结论.由于涉及到时间序列的反复生成,蒙特卡洛模拟法是以高容量和高速度的计算机为前提条件的,因此只是在近些年才得到广泛推广. 从上面的描述我们不难看出随机模拟方法主要是针对那些确定算法不好解或者解不出来的情况.因此

各种距离 欧式距离、曼哈顿距离、切比雪夫距离、闵可夫斯基距离、标准欧氏距离、马氏距离、余弦距离、汉明距离、杰拉德距离、相关距离、信息熵

1. 欧氏距离(Euclidean Distance) 欧氏距离是最容易直观理解的距离度量方法,我们小学.初中和高中接触到的两个点在空间中的距离一般都是指欧氏距离. 二维平面上点a(x1,y1)与b(x2,y2)间的欧氏距离: 三维空间点a(x1,y1,z1)与b(x2,y2,z2)间的欧氏距离: n维空间点a(x11,x12,…,x1n)与b(x21,x22,…,x2n)间的欧氏距离(两个n维向量): Matlab计算欧氏距离: Matlab计算距离使用pdist函数.若X是一个m×n的矩阵,

基于欧氏距离和马氏距离的异常点检测—matlab实现

前几天接的一个小项目,基于欧氏距离和马氏距离的异常点检测,已经交接完毕,现在把代码公开. 基于欧式距离的: load data1.txt %导入数据,行为样本,列为特征 X=data1; %赋值给X u=mean(X); %求均值 [m,n]=size(X); for i=1:m dist(i)=sqrt(sum(X(i,:)-u).^2); end [a,b]=sort(dist);%对欧氏距离进行排序 T=ceil(m*0.02)%设置阀值 Threshold=a(m-T);%定为阀值 le

蒙特卡罗随机模拟

Monte Carlo方法也叫随机模拟.随机抽样或者统计实验方法.其主要用途用于模拟一些无法用数值产生的随机系统.比如当系统的各个单元的特征量已知,但系统过于复杂导致无法预测其准确数学模型,这个时候可以用随机模拟法计算系统的相关参数. 蒙特卡罗的基本思想,为了求解数学.物理.生产管理等方面的问题,首先需要建立一个概率模型或者一个抽样试验来计算所求参数的统计特征,最后通过计算机模拟求出近似解. 我们说蒙特卡罗模拟的是随机现象,那么不同的问题的随机现象是不同的.比如掷骰子,各个点数等概率出现,模拟的