Markovs Chains采样

1 引言

蒙特卡洛积分通过生成符合特定分布的随机变量来近似计算积分值。例如:

上面的式子可以理解成求函数 f(x)f(x) 的期望因此根据大数定理我们生成符合 p(x)p(x) 分布的 N 个随机变量然后用下式来近似计算函数的期望等效于上式积分值。

但是如过分布函数 p(x)过于复杂我们是很难生成符合条件的随机变量的计算机只能生成符合诸如均匀分布正态分布等简单分布的随机变量对于复杂分布我们需要用一些特殊的技巧才能生成其分布这里就用到马尔可夫链了。

2 有限维空间马尔可夫链

Example: Predicting the weather with a finite state-space Markov chain

假如加州大学伯克利分校的天气有三种晴天雾天和雨天这用来模拟状态空间的三个离散的值。这里的天气状态十分稳定所以伯克利的天气预报员可以轻松地根据当前的天气来预测下周的天气按照下面转移概率矩阵进行计算

2.1 直接采用MC公式进行求解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

clear all
close all
clc

P = [0.8  0.15 0.05;   % SUNNY

   0.4  0.5  0.1;    % FOGGY

   0.1  0.3  0.6];    % RAINY

nWeeks = 25;

   

% 初始化状态为RAINY概率为1

X(1,:) = [0 0 1];

   
% RUN MARKOV CHAIN

for iB = 2:nWeeks

    X(iB,:) = X(iB-1,:)*P; % TRANSITION
end

   
% DISPLAY
figure; hold on

h(1) = plot(1:nWeeks,X(:,1),'r','Linewidth',1);   % SUNNY概率

h(2) = plot(1:nWeeks,X(:,2),'k','Linewidth',1);   % FOGGY概率

h(3) = plot(1:nWeeks,X(:,3),'b','Linewidth',1);   % RAINY概率

h(4) = plot([15 15],[0 1],'g--','Linewidth',1);  %暂态阶段热机阶段
hold off

legend(h, {'Sunny','Foggy','Rainy','Burn In'});

xlabel('Week')

ylabel('p(Weather)')

xlim([1,nWeeks]);

ylim([0 1]);

   
% PREDICTIONS

fprintf('\np(weather) in 1 week -->'), disp(X(2,:))

fprintf('\np(weather) in 24 weeks -->'), disp(X(24,:))

fprintf('\np(weather) in 25 weeks -->'), disp(X(25,:))

输出结果：

p(weather) in 1 week --> 0.1000 0.3000 0.6000

p(weather) in 24 weeks --> 0.5965 0.2632 0.1404

p(weather) in 25 weeks --> 0.5965 0.2632 0.1404

2.1 采用MC抽样进行求解

上面的代码展示了经过多次迭代后马尔可夫链的分布趋于一个稳态分布。我们将利用马尔可夫链的这一性质进行抽样，从而得到符合这个稳态分布的样例。

之所以讨论MCMC方法，就是为了得到符合某种分布的样例，通常，我们不容易对这种分布进行直接抽样（知道了分布函数也不行），我们计算机只能对均匀分布、正态分布等一系列简单分布进行抽样，对复杂分布是无法抽样的。所以，我们的思路是：构造这样的马尔可夫链，使其稳态分布符合我们要抽样的函数分布，然后沿着马尔可夫链进行抽样，从而得到符合稳态分布的样例，也就是间接地得到了符合目标分布函数的样例。下面我们继续利用上例的转移矩阵，实地抽样，看看得到的样例分布是否符合其稳态分布。

方法一：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

clear all
close all
clc

P = [0.8 0.15 0.05;   % SUNNY

   0.4 0.5  0.1;   % FOGGY

   0.1 0.3  0.6];  % RAINY

nWeeks = 25000; % 为了精确，我们进行25000次抽样

X = zeros(nWeeks,3); % x 用来存储该周天气

p = zeros(nWeeks,3); % p 用来存储到当前周为止，每种天气出现的频率

X(1,:) = [0,0,1]; % 第一周天气为雨天

for i=2:nWeeks

    px = X(i-1,:)*P; % 根据前一周的天气来计算这一周三种天气出现的概率

    % 我们把三种天气的累计频率计算出来，比如px=[0.3,0.5,0.2]

    % 那么 cump = [0.3,0.8,1]，然后生成0到1之间符合均匀分布的随机数，

    % 这个随机数小于0.3的概率是0.3，大于0.3小于0.8的概率是0.5，大于0.8

    % 小于1的概率是0.2，这样就可以模拟随机抽样的过程了

    cump = cumsum(px);

    % 生成一个在0到1之间均匀分布的随机数

    rnd = rand;

    % 看这个随机数落到哪个区间：如果小于第一个数，那么就是晴天，如果在第一个数

    % 到第二数个之间，就是雾天，如果在第二个数到第三个数之间，就是雨天

    if(rnd<=cump(1))

        X(i,:) = [1,0,0];

    elseif(rnd>cump(1)&&rnd<=cump(2))

        X(i,:) = [0,1,0];

    else

        X(i,:) = [0,0,1];

    end
end

cumx = cumsum(X,1); % 每一行的元素代表到这一周为止，出现的晴天、雾天、雨天的次数
% 每一行各种天气的频数除以当前行对应周数

p = cumx./repmat(transpose(1:nWeeks),1,3);
hold on;

h(1) = plot(1:nWeeks,p(:,1),'r','Linewidth',1);

h(2) = plot(1:nWeeks,p(:,2),'k','Linewidth',1);

h(3) = plot(1:nWeeks,p(:,3),'b','Linewidth',1);
hold off

legend(h, {'Sunny','Foggy','Rainy'});

xlabel('Week')

ylabel('p(Weather)')

xlim([1,nWeeks]);

ylim([0 1]);
% PREDICTIONS

fprintf('\np(weather) in 1 week -->'), disp(p(2,:))

fprintf('\np(weather) in 24999 weeks -->'), disp(p(nWeeks-1,:))

fprintf('\np(weather) in 25000 weeks -->'), disp(p(nWeeks,:))

输出结果：

p(weather) in 1 week --> 0 0.5000 0.5000

p(weather) in 24999 weeks --> 0.6002 0.2569 0.1429

p(weather) in 25000 weeks --> 0.6002 0.2569 0.1429

方法二：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

clear all
close all
clc

P = [.8 .15 .05;  % SUNNY

     .4 .5  .1;   % FOGGY

     .1 .3  .6];  % RAINY

nWeeks = 25000; % 为了精确，我们进行25000次抽样

X = ones(1,nWeeks); % x 用来存储该周天气

p = zeros(nWeeks,3); % p 用来存储到当前周为止，每种天气出现的频率

X(1) = 3; % 第一周天气为雨天

C=cumsum(P,2);

for i=2:nWeeks    

    % 生成一个在0到1之间均匀分布的随机数

    rnd = rand;

    % 看这个随机数落到哪个区间：如果小于第一个数，那么就是晴天，如果在第一个数

    % 到第二数个之间，就是雾天，如果在第二个数到第三个数之间，就是雨天

    if(rnd<=C(X(i-1),1))

        X(i) = 1;

    elseif(rnd>C(X(i-1),1)&&rnd<=C(X(i-1),2))

        X(i) = 2;

    else

        X(i) = 3;

    end
end
X=ind2vec(X);
X=X';

cumx = cumsum(X,1); % 每一行的元素代表到这一周为止，出现的晴天、雾天、雨天的次数
% 每一行各种天气的频数除以当前行对应周数

p = cumx./repmat(transpose(1:nWeeks),1,3);
hold on;

h(1) = plot(1:nWeeks,p(:,1),'r','Linewidth',1);

h(2) = plot(1:nWeeks,p(:,2),'k','Linewidth',1);

h(3) = plot(1:nWeeks,p(:,3),'b','Linewidth',1);
hold off

legend(h, {'Sunny','Foggy','Rainy'});

xlabel('Week')

ylabel('p(Weather)')

xlim([1,nWeeks]);

ylim([0 1]);
% PREDICTIONS

fprintf('\np(weather) in 1 week -->'), disp(p(2,:))

fprintf('\np(weather) in 24999 weeks -->'), disp(p(nWeeks-1,:))

fprintf('\np(weather) in 25000 weeks -->'), disp(p(nWeeks,:))

输出结果：

p(weather) in 1 week --> (1,3) 1

p(weather) in 24999 weeks --> (1,1) 0.6071

(1,2) 0.2536

(1,3) 0.1393

p(weather) in 25000 weeks --> (1,1) 0.6071

(1,2) 0.2536

(1,3) 0.1393

3 连续状态空间马尔可夫链

我们可以通过连续状态空间的马尔可夫链的稳态分布，来对一个连续概率密度函数取样（通常这个概率密度函数都不容易直接取样，原因是函数太过复杂）：我们运行充分次数的马尔可夫链从而让它达到稳态分布，然后保留马尔可夫链到达稳态分布后抽取的样例。

下面的例子里，我们定义一个连续状态空间马尔可夫链。转移算子是正态分布，其方差为1，均值为当前状态与0点的距离的一半，初始状态的分布是均值为0方差为1的标准正态分布。

为了保证链条已经移动得离初始状态足够远，达到了稳态分布，我们将扔掉马尔可夫链最初的50个抽样（也就是 burn in 状态）。我们同时运行多个链，从而更加致密地得到稳态分布抽样。这里我们选择同时运行5个马尔可夫链。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

% EXAMPLE OF CONTINUOUS STATE-SPACE MARKOV CHAIN
clear all
close all
clc  
% INITIALIZE

randn('seed',12345)

nBurnin = 50; % # BURNIN

nChains = 5;  % # MARKOV CHAINS

   
% DEFINE TRANSITION OPERATOR
% 定义内联函数,用来产生符合正态分布的随机变量normrnd(均值,方差,行,列)

P = inline('normrnd(.5*x,1,1,nChains)','x','nChains');

nTransitions = 1000;

% 将所有状态初始化为0，抽样时对这个矩阵进行更新
x = zeros(nTransitions,nChains);
%初始状态，由于是简单的正态分布，所以可以用函数得到

x(1,:) = randn(1,nChains);

   
% RUN THE CHAINS

for iT = 2:nTransitions

    %后一个状态由前一个状态决定，这里源代码等号右边是 P(x(iT-1),nChains)不准确，

    %因为5个马尔可夫链上，每个抽样都是由上一个抽样决定的，而不是由第一个链的上一个状态决定

    x(iT,:) = P(x(iT-1,:),nChains);
end

   
% DISPLAY BURNIN
figure

subplot(221); plot(x(1:100,:)); hold on;
minn = min(x(:));
maxx = max(x(:));

l = line([nBurnin nBurnin],[minn maxx],'color','k','Linewidth',2);
ylim([minn maxx])

legend(l,'~Burn-in','Location','SouthEast')

title('First 100 Samples'); hold off

   
% DISPLAY ENTIRE MARKOV CHAIN

subplot(223); plot(x);hold on;

l = line([nBurnin nBurnin],[minn maxx],'color','k','Linewidth',2);

legend(l,'~Burn-in','Location','SouthEast')

title('Entire Chain');

   
% DISPLAY SAMPLES FROM STATIONARY DISTRIBUTION

samples = x(nBurnin+1:end,:);

subplot(122);

% 这里 hist 将样例的区间分成100份，返回bins为各个区间的中点，counts是落到各个区间的数目

[counts,bins] = hist(samples(:),100); colormap hot
b = bar(bins,counts);

legend(b,sprintf('Markov Chain\nSamples'));

title(['\mu=',num2str(mean(samples(:))),' \sigma=',num2str(var(samples(:)))])

输出结果：

左上方的输出，我们可以看到1000次转移中的前100次，由5个马尔可夫链同时运行得到。burn in 的切割点用黑线显示。左下方可以看到马尔可夫链的整体转移过程。右边的状态分布直方图，表示 burn in 之后的抽样的分布情况，我们可以发现大部分抽样值都集中在-2到2之间，这个是一个正态分布，均值为0，方差为1.3。

4 结语

上面的例子中，我们能够根据 burn in 期之后的抽样推断出马尔可夫链的稳态分布。然而，为了能够利用马尔可夫链来对目标分布函数进行抽样，我们需要设计转移算子，使得经过多次迭代后，链最终达到与目标函数吻合的分布。这就是诸多 MCMC 方法，例如 Metropolis sampler，Metropolis-Hastings sampler，以及 Gibbs sampler 所要实现的目标。

参考：

https://theclevermachine.wordpress.com/2012/09/24/a-brief-introduction-to-markov-chains/

https://theclevermachine.wordpress.com/2012/10/05/mcmc-the-metropolis-sampler/

https://theclevermachine.wordpress.com/2012/10/20/mcmc-the-metropolis-hastings-sampler/

http://www.cnblogs.com/yinxiangnan-charles/p/5018876.html

本文转自stock0991 51CTO博客，原文链接：http://blog.51cto.com/qing0991/1793526，如需转载请自行联系原作者

相关推荐