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

Tianyi_Li 发表于 2021/12/06 09:20:03 2021/12/06
【摘要】 一、实验目的了解马尔科夫抽样与M-H抽样的原理,阅读已有程序代码实现通过MCMC方法获取指定概率分布模拟样本的目的。二、实验内容: 阅读已有程序代码。理解程序中接受-拒绝抽样的原理。     3.了解马尔科夫抽样与M-H抽样的原理。     4. 实现MCMC方法获取指定概率分布的模拟样本。     5. 输出获取样本的概率分布图,验证结果。三、实验程序及结果clear; close all...

一、实验目的

了解马尔科夫抽样与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方法获取指定概率分布模拟样本的目的。

 

【版权声明】本文为华为云社区用户原创内容,转载时必须标注文章的来源(华为云社区),文章链接,文章作者等基本信息,否则作者和本社区有权追究责任。如果您发现本社区中有涉嫌抄袭的内容,欢迎发送邮件至:cloudbbs@huaweicloud.com进行举报,并提供相关证据,一经查实,本社区将立刻删除涉嫌侵权内容。
  • 点赞
  • 收藏
  • 关注作者

评论(0

0/1000
抱歉,系统识别当前为高风险访问,暂不支持该操作

全部回复

上滑加载中

设置昵称

在此一键设置昵称,即可参与社区互动!

*长度不超过10个汉字或20个英文字符,设置后3个月内不可修改。

*长度不超过10个汉字或20个英文字符,设置后3个月内不可修改。