本文是利用蒙特卡罗算法对马氏链过程的模拟。假设有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]);仿真过程如下:
最终的结果如下:
原文:http://blog.csdn.net/tengweitw/article/details/34063833
作者:nineheadedbird