蒙特卡罗方法生成指定状态空间下对应长度的马尔可夫链--MATLAB源程序

Posted Z.Q.Feng

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了蒙特卡罗方法生成指定状态空间下对应长度的马尔可夫链--MATLAB源程序相关的知识,希望对你有一定的参考价值。

MCMC方法生成指定状态空间下对应长度的马尔可夫链


马尔可夫链

定义:


生成马尔可夫链的蒙特卡罗方法(MCMC方法)

详细原理可参考百度百科:马尔科夫链蒙特卡洛方法

在这里我们只介绍MCMC算法中的Metropolis-Hastings算法

Hastings-Metropolis Algorithm
        Step1:  Choose an irreducible Markov transition 
        		probability matrix Q with transition probabilities 	
        		Q(i,j), i, j = 1, 2, ...,m. Also, choose
                some integer value k between 1 and m.
        Step2:  Let n = 0 and X0 = k
        Step3:  Generate a random variable X such that 
        		Pr{X = j} = Q(Xn, j) and 
        		generate a random number U.
        Step4:  If U < [Pie(X)Q(X,Xn)]/[Pie(X)Q(Xn,X)], 
        		then NS = X; else NS = Xn.
        Step5:  n = n + 1, Xn = NS;
        Step6:  Go to Step3.

伪码如下:


代码实现

MATLAB源码:

function S = random_MarkovChain(Pie,N)
%  -----------------------------------------------------------------------
%  Function S = random_MarkovChain(Pie,N) generates a sequence   
%    of the Markov Chain with transition probability matrix Q  
%    and stationary distribution Pie by the Hastings-Metropolis 
%    algorithm.
%  -----------------------------------------------------------------------
%    Input 
%          Pie  a probability vector
%          N    the length of the sequence S, its default 
%				value is 4000
%    Output
%          S    a sequence of Markov Chain with transition
%				matrix Q 
%    Remark
%          state space: {1,2,...,m}
%  -----------------------------------------------------------------------
%  by Feng Zhiqiang

if nargin == 1
    N = 4000;
end

m = max(size(Pie));

% Q is the transition probability matrix, 
% here Q(Xn,:) is a uniform distribution
Q = ones(m)/m; 

% Set the initial state
k = floor(rand(1,1)*m) + 1;
Xn = k;

S_long = zeros(1,3*N);
S_long(1) = Xn;

for n = 1:3*N
   
    % for a general Q(Xn,:), 
    % X = random_AccRejMethod(Q(Xn,:)) is OK;
    X = floor(rand(1,1)*m) + 1; 
    
    U = rand(1,1);
    if U < Pie(X)*Q(X, Xn)/(Pie(Xn)*Q(Xn,X))
        NS = X;
    else
        NS = Xn;
    end
    
    S_long(n) = NS;
    Xn = NS;
end

%   N, N+1, ..., 2N -1, there are N elements in the sequence S
S = S_long(N:2*N-1);  

return

输出示例


调用 random_MarkovChain 函数,输出结果如下:

>> p = [ 0.1 0.2 0.1 0.3 0.3];
>> X = random_MarkovChain(p,14)
X =
	5    2    4    4    4    4    5    4    4    5    1    4    3    5

最优化相关算法设计数学原理:最优化/Optimization文章合集


有帮助可以点赞哦,谢谢大家的支持~

以上是关于蒙特卡罗方法生成指定状态空间下对应长度的马尔可夫链--MATLAB源程序的主要内容,如果未能解决你的问题,请参考以下文章

一份数学小白也能读懂的「马尔可夫链蒙特卡洛方法」入门指南

隐马尔可夫模型

R语言中实现马尔可夫链蒙特卡罗MCMC模型

马尔可夫链蒙特卡罗法

隐马尔可夫模型

在 Swift 中使用马尔可夫链生成文本