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