0

我正在尝试模拟参数 theta 的分布f= theta ^(z_f+n+alpha-1)*(1-theta)^(n+1-z_f-k+ beta-1)，其中除 theta 之外的所有参数都是已知的。我正在使用 Metro Polish hasting 算法进行 MCMC 模拟。我的提案密度是带有参数 alpha 和 beta 的 beta 分布。我的模拟代码如下。为此，我正在使用名为 mhsample() 的 buitlin Matlab 代码，我如何知道我的代码是否正常工作？

clear
clc
alpha=2;
beta=2;
z_f=1;
n=6;
k=5;

nsamples = 3000;
pdf= @(x) x^(z_f+n+alpha-1)*(1-x)^(n+1-z_f-k+beta-1); % here x acts as theta
proppdf= @(x,y) betapdf(x, alpha, beta);
proprnd =@(x) betarnd(alpha,beta,1);

smpl = mhsample(0.1,nsamples,'pdf',pdf,'proprnd',proprnd,'proppdf',proppdf);
4

1 回答 1

0

当你说“我怎么知道我的代码是否正常工作”时，我不确定你在问什么——我假设它执行了？但是为了直观地比较您的函数与模拟，您可以绘制 PDF 和从 mhsample 获得的数据，如下所示：

% i'm assuming you ran the code above so that smpl and @pdf are both defined...

fplot(pdf,[0 1]); % fplot takes your function and plots it between x-limit [0,1]
figure % new figure
hist(smpl,30); % 30 here is bin size, change it to your preference

下图：

  • 左侧输出的直方图smpl，即您的模拟
  • pdf右侧 [0,1] 中的函数，用于与您的模拟进行比较

左边是 smpl 的输出，右边是 [0,1] 中的 pdf

这只是一个疯狂的猜测，因为这两个数字彼此相似，并且也是 beta-distribution-esque。

如果您想要比这更复杂的分析，恐怕我还不精通 MCMC :)

于 2014-08-26T23:07:36.723 回答