MCMC方法获取指定概率分布的模拟样本实验的MATLAB实现（附代码）

一、实验目的

了解马尔科夫抽样与M-H抽样的原理，阅读已有程序代码实现通过MCMC方法获取指定概率分布模拟样本的目的。以下介绍来自网络搜索：

MCMC方法就是构造合适的马尔科夫链进行抽样而使用蒙特卡洛方法进行积分计算，既然马尔科夫链可以收敛到平稳分布。我们可以建立一个以π为平稳分布的马尔科夫链，对这个链运行足够长时间之后，可以达到平稳状态。此时马尔科夫链的值就相当于在分布π(x)中抽取样本。利用马尔科夫链进行随机模拟的方法就是MCMC。

第一个MC: Monte Carlo(蒙特卡洛)。这个简单来说是让我们使用随机数（随机抽样）来解决计算问题。在MCMC中意味着：后验分布作为一个随机样本生成器，我们利用它来生成样本(simulation)，然后通过这些样本对一些感兴趣的计算问题（特征数，预测）进行估计。

第二个MC：Markov Chain(马尔科夫链)。第二个MC是这个方法的关键，因为我们在第一个MC中看到，我们需要利用后验分布生成随机样本，但后验分布太复杂，当这些样本独立时，利用大数定律样本均值会收敛到期望值。如果得到的样本是不独立的，那么就要借助于马尔科夫链进行抽样，利用Markov Chain的平稳分布这个概念实现对复杂后验分布的抽样。

二、实验内容：

1. 阅读已有程序代码。
2. 理解程序中接受-拒绝抽样的原理。

     3.了解马尔科夫抽样与M-H抽样的原理。

     4. 实现MCMC方法获取指定概率分布的模拟样本。

     5. 输出获取样本的概率分布图，验证结果。

三、实验程序及结果

clear; close all; home; format long g;
 
%% Initial data
% Distributions 
nu           = 5;                                       % Target PDF parameter
p            = @(t) (t.^(-nu/2-1)).*(exp(-1./(2*t)));   % Target "PDF", Ref. 1 - Ex. 2
proposal_PDF = @(mu) unifrnd(0,3);                      % Proposal PDF
% Parameters
N        = 1000;              % Number of samples (iterations)
nn       = 0.1*(N);           % Number of samples for examine the AutoCorr
theta    = zeros(1,N);        % Samples drawn form the target "PDF"
theta(1) = 0.3;               % Initial state of the chain
 
%% Procedure
for i = 1:N
   theta_ast = proposal_PDF([]);          % Sampling from the proposal PDF
   alpha     = p(theta_ast)/p(theta(i));  % Ratio of the density at theta_ast and theta points
   if rand <= min(alpha,1)
      % Accept the sample with prob = min(alpha,1)
      theta(i+1) = theta_ast;
   else
      % Reject the sample with prob = 1 - min(alpha,1)
      theta(i+1) = theta(i);
   end
end
 
 
 
% Samples and distributions
xx = eps:0.01:10;   % x-axis (Graphs)
figure;
% Histogram and target dist
% subplot(2,1,1); 
[n1, x1] = hist(theta, ceil(sqrt(N)));
bar(x1,n1/(N*(x1(2)-x1(1))));   colormap summer;   hold on;   % Normalized histogram
plot(xx, p(xx)/trapz(xx,p(xx)), 'r-', 'LineWidth', 3);        % Normalized function
grid on;   xlim([0 max(theta)]);   
title('Distribution of samples', 'FontSize', 14);
ylabel('Probability density function', 'FontSize', 12);

四、心得体会

掌握了马尔科夫抽样与M-H抽样的原理，了解了通过MCMC方法获取指定概率分布模拟样本的目的。