-
Notifications
You must be signed in to change notification settings - Fork 8
/
RJMCMCStep.m
40 lines (35 loc) · 1.25 KB
/
RJMCMCStep.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
35
36
37
38
39
function [newState, draw] = RJMCMCStep(oldState, settings)
newState = oldState;
%Get draw
draw = getDraw(oldState, settings);
%Evaluate joint log posterior
[draw.logPosterior] = evaluatePosterior(draw, settings);
%Compute log acceptance probability
alpha = draw.logPosterior - oldState.logPosterior + evaluateProposalRatio(oldState, draw, settings);
%Accept draw with probability alpha
if rand <= exp(min(alpha,0))
% Accept the candidate
newState.logPosterior = draw.logPosterior;
if draw.ps > 0
newState.arParameters = draw.arParameters;
newState.arPacs = draw.arPacs;
else
newState.arParameters = [];
newState.arPacs = [];
end;
if draw.qs > 0
newState.maParameters = draw.maParameters;
newState.maPacs = draw.maPacs;
else
newState.maParameters = [];
newState.maPacs = [];
end;
newState.ps = draw.ps;
newState.qs = draw.qs;
draw.accepted = 1; % Note the acceptance
newState.sigmaEs = draw.sigmaEs;
else
%reject the candidate
draw.accepted = 0; % Note the rejection
end;
end