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

Posted 2021-07-25 Z.Q.Feng

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

• 马尔可夫链
• 生成马尔可夫链的蒙特卡罗方法(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文章合集

---

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