where supp(π) = {x : π(x) > 0} is the support of distribution π (i.e., the set of points with non-zero probability).
This condition just says that our proposal must have a non-zero probability of moving to the states that have non-zero
probability in the target. The resulting transition distribution of the algorithm is
p
M
(x
0
|x) = q(x
0
|x)r(x
0
|x) + I(x
0
= x)(1 −
Z
q(x
0
|x)r(x
0
|x)dx
0
) (20)
An easy way to implement step (c) is to generate U ∼ U(0, 1) and to set X
s+1
= x
0
if U < r, and to set
x
s+1
= X
s
otherwise.
Note that when evaluating α, we only need to knowthe target density p up to a normalization constant. In particular,
suppose π(x) =
1
Z
π
0
(x), where π
0
(x) is an unnormalized distribution and Z is the normalization constant. Then
α =
(π
0
(x
0
)/Z) q(x|x
0
)
(π
0
(x)/Z) q(x
0
|x)
(21)
so the Z’s cancel. Thus, using MH, we can sample from π even if we can only compute π
0
. Later we will see many
examples where it is hard to evaluate Z but easy to evaluate π
0
.
If we have a symmetric proposal distribution q(x
0
|x) = q(x|x
0
), then the acceptance ratio simplifies to
α =
π(x
0
)
π(x)
(22)
This is called the Metropolis algorithm. For example, it is common to use a Gaussian as a proposal distribution:
q(x
0
|x) = N(x
0
|x, σ
2
). This is called a random walk MH algorithm. It is crucial to pick the right σ
2
to ensure that
a reasonable number (say 50%) of the proposals are accepted: see Figure 3.
In the Metropolis algorithm, if the new state x
0
is more probable than the current state x, the proposal is always
accepted r(x
0
|x) = 1, otherwise it is accepted with probability π(x
0
)/π(x).
A special case of the Metropolis algorithm is when the proposal is independent of the current state: q(x
0
|x) =
q(x
0
). Then the acceptance probability is
α =
π(x
0
)/q(x
0
)
π(x)/q(x)
(23)
This is called the independence sampler, and is similar to importance sampling.
Below is some generic MH code. If the proposal distribution is symmetric, it is not necessary to compute the actual
probabilities q(x
0
|x) and q(x|x
0
), it is only necessary to be able to sample from q(x
0
|x).
function [samples, naccept] = MH(target, proposal, xinit, Nsamples, targetArgs, proposalArgs, proposalProb)
% Metropolis-Hastings algorithm
%
% Inputs
% target returns the unnormalized log posterior, called as ’p = exp(target(x, targetArgs{:}))’
% proposal is a fn, as ’xprime = proposal(x, proposalArgs{:})’ where x is a 1xd vector
% xinit is a 1xd vector specifying the initial state
% Nsamples - total number of samples to draw
% targetArgs - cell array passed to target
% proposalArgs - cell array passed to proposal
% proposalProb - optional fn, called as ’p = proposalProb(x,xprime, proposalArgs{:})’,
% computes q(xprime|x). Not needed for symmetric proposals (Metropolis algorithm)
%
% Outputs
% samples(s,:) is the s’th sample (of size d)
% naccept = number of accepted moves
if nargin < 5, targetArgs = {}; end
if nargin < 6, proposalArgs = {}; end
if nargin < 7, proposalProb = []; end
d = length(xinit);
samples = zeros(Nsamples, d);
x = xinit(:)’;
naccept = 0;
logpOld = feval(target, x, targetArgs{:});
for t=1:Nsamples
4
剩余18页未读，继续阅读
瑞脑销金兽
- 粉丝: 0
- 资源: 6
我的内容管理
展开
最新资源
- AirKiss技术详解：无线传递信息与智能家居连接
- Hibernate主键生成策略详解
- 操作系统实验：位示图法管理磁盘空闲空间
- JSON详解：数据交换的主流格式
- Win7安装Ubuntu双系统详细指南
- FPGA内部结构与工作原理探索
- 信用评分模型解析：WOE、IV与ROC
- 使用LVS+Keepalived构建高可用负载均衡集群
- 微信小程序驱动餐饮与服装业创新转型：便捷管理与低成本优势
- 机器学习入门指南：从基础到进阶
- 解决Win7 IIS配置错误500.22与0x80070032
- SQL-DFS：优化HDFS小文件存储的解决方案
- Hadoop、Hbase、Spark环境部署与主机配置详解
- Kisso：加密会话Cookie实现的单点登录SSO
- OpenCV读取与拼接多幅图像教程
- QT实战：轻松生成与解析JSON数据
资源上传下载、课程学习等过程中有任何疑问或建议，欢迎提出宝贵意见哦~我们会及时处理！
点击此处反馈
