-
Notifications
You must be signed in to change notification settings - Fork 0
/
slice_sample.m
64 lines (59 loc) · 1.61 KB
/
slice_sample.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
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
function [xnew, varargout] = slice_sample(x,f,b,varargin)
% b: a vector specifying an interval, a scalar specifying a width or a cell
% array of the form: {@sampling_function [l,r] {arguments}}
M = 5; % maximum width = M*b
% draw a level to define a slice. This assumes that f(x) = log(p(x)).
fx = feval(f, x, varargin{:});
y = fx + log(rand); % y = log(p(x)) - e, where e ~ exp(1)
if iscell(b)
%l = -1/eps;
%r = 1/eps;
l = b{2}(1);
r = b{2}(2);
g = b{1};
if length(b) > 2
gargs = b{3};
else
gargs = {};
end
else
if length(b) == 1
% step out to find a bound large enough to contain the slice
l = x - b*rand; % left bound
r = l + b; % right bound
j = floor(M*rand);
k = M - 1 - j;
% left bound
while j > 0 && y < feval(f, l, varargin{:})
l = l - b;
j = j - 1;
end
% right bound
while k > 0 && y < feval(f, r, varargin{:})
r = r + b;
k = k - 1;
end
else
% just use the bounds provided
l = b(1);
r = b(2);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% shrink the bound to sample from the slice
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ll = l; rr = r; vout = cell(nargout-1,1);
while true
if iscell(b)
xnew = feval(g,ll,rr,gargs{:});
else
xnew = ll + rand*(rr-ll);
end
[fxnew, vout{:}] = feval(f, xnew, varargin{:});
if y < fxnew
varargout = vout;
break
end
if xnew < x, ll = xnew; else rr = xnew; end
end
end