1########################################################################
2##
3## Copyright (C) 1993-2021 The Octave Project Developers
4##
5## See the file COPYRIGHT.md in the top-level directory of this
6## distribution or <https://octave.org/copyright/>.
7##
8## This file is part of Octave.
9##
10## Octave is free software: you can redistribute it and/or modify it
11## under the terms of the GNU General Public License as published by
12## the Free Software Foundation, either version 3 of the License, or
13## (at your option) any later version.
14##
15## Octave is distributed in the hope that it will be useful, but
16## WITHOUT ANY WARRANTY; without even the implied warranty of
17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18## GNU General Public License for more details.
19##
20## You should have received a copy of the GNU General Public License
21## along with Octave; see the file COPYING.  If not, see
22## <https://www.gnu.org/licenses/>.
23##
24########################################################################
25
26## -*- texinfo -*-
27## @deftypefn  {} {[@var{smpl}, @var{neval}] =} slicesample (@var{start}, @var{nsamples}, @var{property}, @var{value}, @dots{})
28## Draws @var{nsamples} samples from a target stationary distribution @var{pdf}
29## using slice sampling of Radford M. Neal.
30##
31## Input:
32## @itemize
33## @item
34## @var{start} is a 1 by @var{dim} vector of the starting point of the
35## Markov chain. Each column corresponds to a different dimension.
36##
37## @item
38## @var{nsamples} is the number of samples, the length of the Markov chain.
39## @end itemize
40##
41## Next, several property-value pairs can or must be specified, they are:
42##
43## (Required properties) One of:
44##
45## @itemize
46## @item
47## @var{"pdf"}: the value is a function handle of the target stationary
48## distribution to be sampled.  The function should accept different locations
49## in each row and each column corresponds to a different dimension.
50##
51## or
52##
53## @item
54## @var{logpdf}: the value is a function handle of the log of the target
55## stationary distribution to be sampled. The function should accept different
56## locations in each row and each column corresponds to a different dimension.
57## @end itemize
58##
59## The following input property/pair values may be needed depending on the
60## desired outut:
61##
62## @itemize
63## @item
64## "burnin" @var{burnin} the number of points to discard at the beginning, the default
65## is 0.
66##
67## @item
68## "thin" @var{thin} omitts @var{m}-1 of every @var{m} points in the generated
69## Markov chain. The default is 1.
70##
71## @item
72## "width" @var{width} the maximum Manhattan distance between two samples.
73## The default is 10.
74## @end itemize
75##
76## Outputs:
77## @itemize
78##
79## @item
80## @var{smpl} is a @var{nsamples} by @var{dim} matrix of random
81## values drawn from @var{pdf} where the rows are different random values, the
82## columns correspond to the dimensions of @var{pdf}.
83##
84## @item
85## @var{neval} is the number of function evaluations per sample.
86## @end itemize
87## Example : Sampling from a normal distribution
88##
89## @example
90## @group
91## start = 1;
92## nsamples = 1e3;
93## pdf = @@(x) exp (-.5 * x .^ 2) / (pi ^ .5 * 2 ^ .5);
94## [smpl, accept] = slicesample (start, nsamples, "pdf", pdf, "thin", 4);
95## histfit (smpl);
96## @end group
97## @end example
98##
99## @seealso{rand, mhsample, randsample}
100## @end deftypefn
101
102function [smpl, neval] = slicesample (start, nsamples, varargin)
103
104  if (nargin < 4)
105    print_usage ();
106  endif
107
108  sizestart = size (start);
109  pdf = [];
110  logpdf = [];
111  width = 10;
112  burnin = 0;
113  thin = 1;
114  for k = 1:2:length (varargin)
115    if (ischar (varargin{k}))
116      switch lower (varargin{k})
117        case "pdf"
118          if (isa (varargin{k+1}, "function_handle"))
119            pdf = varargin{k+1};
120          else
121            error ("slicesample: pdf must be a function handle");
122          endif
123        case "logpdf"
124          if (isa (varargin{k+1}, "function_handle"))
125            pdf = varargin{k+1};
126          else
127            error ("slicesample: logpdf must be a function handle");
128          endif
129        case "width"
130          if (numel (varargin{k+1}) == 1 || numel (varargin{k+1}) == sizestart(2))
131            width = varargin{k+1}(:).';
132          else
133            error ("slicesample: width must be a scalar or 1 by dim vector");
134          endif
135        case "burnin"
136          if (varargin{k+1}>=0)
137            burnin = varargin{k+1};
138          else
139            error ("slicesample: burnin must be greater than or equal to 0");
140          endif
141        case "thin"
142          if (varargin{k+1}>=1)
143            thin = varargin{k+1};
144          else
145            error ("slicesample: thin must be greater than or equal to 1");
146          endif
147        otherwise
148          warning (["slicesample: Ignoring unknown option " varargin{k}]);
149      endswitch
150    else
151      error (["slicesample: " varargin{k} " is not a valid property."]);
152    endif
153  endfor
154
155  if (! isempty (pdf) && isempty (logpdf))
156    logpdf = @(x) rloge (pdf (x));
157  elseif (isempty (pdf) && isempty (logpdf))
158    error ("slicesample: pdf or logpdf must be input.");
159  endif
160  dim = sizestart(2);
161  smpl = zeros (nsamples, dim);
162
163  if (all (sizestart == [1 dim]))
164    smpl(1, :) = start;
165  else
166    error ("slicesample: start must be a 1 by dim vector.");
167  endif
168
169  maxit = 100;
170  neval = 0;
171
172  fgraterthan = @(x, fxc) logpdf (x) >= fxc;
173
174  ti = burnin + nsamples * thin;
175
176  rndexp = rande (ti, 1);
177  crand = rand (ti, dim);
178  prand = rand (ti, dim);
179
180  xc = smpl(1, :);
181  for i = 1:ti
182    neval++;
183    sliceheight = logpdf (xc) - rndexp(i);
184    c = width .* crand(i, :);
185    lb = xc - c;
186    ub = xc + width - c;
187    #Only for single variable as bounds can not be found with point when dim > 1
188    if (dim == 1)
189      for k=1:maxit
190        neval++;
191        if (! fgraterthan (lb, sliceheight))
192          break
193        endif
194        lb -= width;
195      end
196      if (k == maxit)
197        warning ("slicesample: Step out exceeded maximum iterations");
198      endif
199      for k = 1:maxit
200        neval++;
201        if (! fgraterthan (ub, sliceheight))
202          break
203        endif
204        ub += width;
205      end
206      if (k == maxit)
207        warning ("slicesample: Step out exceeded maximum iterations");
208      endif
209    end
210    xp = (ub - lb) .* prand(i, :) + lb;
211    for k=1:maxit
212      neval++;
213      isgt = fgraterthan (xp,sliceheight);
214      if (all (isgt))
215        break
216      endif
217      lc = ! isgt & xp < xc;
218      uc = ! isgt & xp > xc;
219      lb(lc) = xp(lc);
220      ub(uc) = xp(uc);
221      xp = (ub - lb) .* rand (1, dim) + lb;
222    end
223    if (k == maxit)
224      warning ("slicesample: Step in exceeded maximum iterations");
225    endif
226    xc = xp;
227    if (i > burnin)
228      indx = (i - burnin) / thin;
229      if rem (indx, 1) == 0
230        smpl(indx, :) = xc;
231      end
232    end
233  end
234  neval = neval / (nsamples * thin + burnin);
235endfunction
236
237function y = rloge (x)
238
239  y = -inf (size (x));
240  xg0 = x > 0;
241  y(xg0) = log (x(xg0));
242
243endfunction
244
245
246%!demo
247%! ## Define function to sample
248%! d = 2;
249%! mu = [-1; 2];
250%! Sigma = rand (d);
251%! Sigma = (Sigma + Sigma');
252%! Sigma += eye (d)*abs (eigs (Sigma, 1, "sa")) * 1.1;
253%! pdf = @(x)(2*pi)^(-d/2)*det(Sigma)^-.5*exp(-.5*sum((x.'-mu).*(Sigma\(x.'-mu)),1));
254%! ##Inputs
255%! start = ones (1,2);
256%! nsamples = 500;
257%! K = 500;
258%! m = 10;
259%! [smpl, accept]=slicesample (start, nsamples, "pdf", pdf, "burnin", K, "thin", m, "width", [20, 30]);
260%! figure;
261%! hold on;
262%! plot (smpl(:,1), smpl(:,2), 'x');
263%! [x, y] = meshgrid (linspace (-6,4), linspace(-3,7));
264%! z = reshape (pdf ([x(:), y(:)]), size(x));
265%! mesh (x, y, z, "facecolor", "None");
266%! ## Using sample points to find the volume of half a sphere with radius of .5
267%! f = @(x) ((.25-(x(:,1)+1).^2-(x(:,2)-2).^2).^.5.*(((x(:,1)+1).^2+(x(:,2)-2).^2)<.25)).';
268%! int = mean (f (smpl) ./ pdf (smpl));
269%! errest = std (f (smpl) ./ pdf (smpl)) / nsamples^.5;
270%! trueerr = abs (2/3*pi*.25^(3/2)-int);
271%! fprintf("Monte Carlo integral estimate int f(x) dx = %f\n", int);
272%! fprintf("Monte Carlo integral error estimate %f\n", errest);
273%! fprintf("The actual error %f\n", trueerr);
274%! mesh (x,y,reshape (f([x(:), y(:)]), size(x)), "facecolor", "None");
275
276%!demo
277%! ##Integrate truncated normal distribution to find normilization constant
278%! pdf = @(x) exp (-.5*x.^2)/(pi^.5*2^.5);
279%! nsamples = 1e3;
280%! [smpl,accept] = slicesample (1, nsamples, "pdf", pdf, "thin", 4);
281%! f = @(x) exp (-.5 * x .^ 2) .* (x >= -2 & x <= 2);
282%! x=linspace(-3,3,1000);
283%! area(x,f(x));
284%! xlabel ('x');
285%! ylabel ('f(x)');
286%! int = mean (f (smpl)./pdf(smpl));
287%! errest = std (f (smpl)./pdf(smpl))/nsamples^.5;
288%! trueerr = abs (erf (2^.5)*2^.5*pi^.5-int);
289%! fprintf("Monte Carlo integral estimate int f(x) dx = %f\n", int);
290%! fprintf("Monte Carlo integral error estimate %f\n", errest);
291%! fprintf("The actual error %f\n", trueerr);
292
293
294%!test
295%! start = 0.5;
296%! nsamples = 1e3;
297%! pdf = @(x) exp (-.5*(x-1).^2)/(2*pi)^.5;
298%! [smpl, accept] = slicesample (start, nsamples, "pdf", pdf, "thin", 2, "burnin", 0, "width", 5);
299%! assert (mean (smpl, 1), 1, .1);
300%! assert (var (smpl, 1), 1, .1);
301
302%!error slicesample ();
303%!error slicesample (1);
304%!error slicesample (1, 1);
305
306