matlab 马尔可夫链代码

写了一个简单的马尔可夫链的模型。

function chain = markov_mine(P,P0,chain_length)
%马尔可夫链模型
%初态由P0的概率随机产生。P0为各个状态的初始概率
%用于确定马尔可夫链的第一个状态是哪一个。
%之后的状态由P来产生。P为概率转移矩阵
%最终马尔可夫链的长度由chain_length来决定
%Author:JayimLiu 12/16/2021
chain=zeros(1,chain_length);%最后想输出的马尔可夫链
prob=cumsum(P0);%计算累积概率
choice=rand();%生成一个随机数
temp=find(prob>=choice);%找到累积概率大于这个随机数的位置
target=temp(1);%取出第一个位置
chain(1)=target;%保存初态
for i=1:chain_length%算之后的链
    if i~=1
        P_temp=P(chain(i-1),:);%用上一个态的位置，锁定需要使用的转移概率
        prob=cumsum(P_temp);
        choice=rand();
        temp=find(prob>=choice);
        target=temp(1);
        chain(i)=target;%保存这次的结果
    end
end
end