-
Notifications
You must be signed in to change notification settings - Fork 188
/
metropolis.m
34 lines (30 loc) · 1.03 KB
/
metropolis.m
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
function [samples, naccept] = metropolis(target, proposal, xinit, Nsamples, targetArgs, proposalArgs)
% Metropolis.m
% Metropolis algorithm (assumes a symmetric proposal distribution)
% function samples = metropolis(target, proposal, xinit, Nsamples, targetArgs, proposalArgs)
%
% Inputs
% target is a function which will be called as 'p = target(x, targetArgs{:})'
% proposal is a function which will be called as 'xprime = proposal(x, proposalArgs{:})' where x is a 1xd vector
% xinit is a 1xd vector specifying the initial state
% Nsamples
%
% Outputs
% samples(t,:) is the t'th sample (of size d)
% naccept = number of accepted moves
if nargin < 5, targetArgs = {}; end
if nargin < 6, proposalArgs = {}; end
d = length(xinit);
samples = zeros(Nsamples, d);
x = xinit(:)';
naccept = 0;
for t=1:Nsamples
xprime = feval(proposal, x, proposalArgs{:});
alpha = min(1, feval(target, xprime, targetArgs{:})/feval(target, x, targetArgs{:}));
u = rand(1,1);
if u < alpha
x = xprime;
naccept = naccept + 1;
end
samples(t,:) = x;
end