########################################################################
##
## Copyright (C) 1993-2021 The Octave Project Developers
##
## See the file COPYRIGHT.md in the top-level directory of this
## distribution or .
##
## This file is part of Octave.
##
## Octave is free software: you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING. If not, see
## .
##
########################################################################
## -*- texinfo -*-
## @deftypefn {} {[@var{smpl}, @var{neval}] =} slicesample (@var{start}, @var{nsamples}, "pdf", @var{pdf}, @var{property}, @var{value})
## Draws @var{nsamples} samples from a target stationary distribution @var{pdf}
## using slice sampling of Radford M. Neal.
##
## Input:
## @itemize
## @item
## @var{start} is a 1 by @var{dim} vector of the starting point of the
## Markov chain. Each column corresponds to a different dimension.
##
## @item
## @var{nsamples} is the number of samples, the length of the Markov chain.
## @end itemize
##
## Next, several property-value pairs can or must be specified, they are:
##
## (Required properties) one of:
##
## @itemize
## @item
## @var{"pdf"}: the value is a function handle of the target stationary
## distribution to be sampled. The function should accept different locations
## in each row and each column corresponds to a different dimension.
##
## or
##
## @item
## @var{logpdf}: the value is a function handle of the log of the target
## stationary distribution to be sampled. The function should accept different
## locations in each row and each column corresponds to a different dimension.
## @end itemize
##
## The following input property/pair values may be needed depending on the
## desired outut:
##
## @itemize
## @item
## "burnin" @var{burnin} the number of points to discard at the beginning, the default
## is 0.
##
## @item
## "thin" @var{thin} omitts @var{m}-1 of every @var{m} points in the generated
## Markov chain. The default is 1.
##
## @item
## "width" @var{width} the maximum Manhattan distance between two samples.
## The default is 10.
## @end itemize
##
## Outputs:
## @itemize
##
## @item
## @var{smpl} is a @var{nsamples} by @var{dim} matrix of random
## values drawn from @var{pdf} where the rows are different random values, the
## columns correspond to the dimensions of @var{pdf}.
##
## @item
## @var{neval} is the number of function evaluations per sample.
## @end itemize
## Example : Sampling from a normal distribution
##
## @example
## @group
## start = 1;
## nsamples = 1e3;
## pdf = @@(x) exp (-.5 * x .^ 2) / (pi ^ .5 * 2 ^ .5);
## [smpl, accept] = slicesample (start, nsamples, "pdf", pdf, "thin", 4);
## histfit (smpl);
## @end group
## @end example
##
## @seealso{rand, mhsample, randsample}
## @end deftypefn
function [smpl, neval] = slicesample (start, nsamples, varargin)
if (nargin < 4)
print_usage ();
endif
sizestart = size (start);
pdf = [];
logpdf = [];
width = 10;
burnin = 0;
thin = 1;
for k = 1:2:length (varargin)
if (ischar (varargin{k}))
switch lower (varargin{k})
case "pdf"
if (isa (varargin{k+1}, "function_handle"))
pdf = varargin{k+1};
else
error ("slicesample: pdf must be a function handle");
endif
case "logpdf"
if (isa (varargin{k+1}, "function_handle"))
pdf = varargin{k+1};
else
error ("slicesample: logpdf must be a function handle");
endif
case "width"
if (numel (varargin{k+1}) == 1 || numel (varargin{k+1}) == sizestart(2))
width = varargin{k+1}(:).';
else
error ("slicesample: width must be a scalar or 1 by dim vector");
endif
case "burnin"
if (varargin{k+1}>=0)
burnin = varargin{k+1};
else
error ("slicesample: burnin must be greater than or equal to 0");
endif
case "thin"
if (varargin{k+1}>=1)
thin = varargin{k+1};
else
error ("slicesample: thin must be greater than or equal to 1");
endif
otherwise
warning (["slicesample: Ignoring unknown option " varargin{k}]);
endswitch
else
error (["slicesample: " varargin{k} " is not a valid property."]);
endif
endfor
if (! isempty (pdf) && isempty (logpdf))
logpdf = @(x) rloge (pdf (x));
elseif (isempty (pdf) && isempty (logpdf))
error ("slicesample: pdf or logpdf must be input.");
endif
dim = sizestart(2);
smpl = zeros (nsamples, dim);
if (all (sizestart == [1 dim]))
smpl(1, :) = start;
else
error ("slicesample: start must be a 1 by dim vector.");
endif
maxit = 100;
neval = 0;
fgraterthan = @(x, fxc) logpdf (x) >= fxc;
ti = burnin + nsamples * thin;
rndexp = rande (ti, 1);
crand = rand (ti, dim);
prand = rand (ti, dim);
xc = smpl(1, :);
for i = 1:ti
neval++;
sliceheight = logpdf (xc) - rndexp(i);
c = width .* crand(i, :);
lb = xc - c;
ub = xc + width - c;
#Only for single variable as bounds can not be found with point when dim > 1
if (dim == 1)
for k=1:maxit
neval++;
if (! fgraterthan (lb, sliceheight))
break
endif
lb -= width;
end
if (k == maxit)
warning ("slicesample: Step out exceeded maximum iterations");
endif
for k = 1:maxit
neval++;
if (! fgraterthan (ub, sliceheight))
break
endif
ub += width;
end
if (k == maxit)
warning ("slicesample: Step out exceeded maximum iterations");
endif
end
xp = (ub - lb) .* prand(i, :) + lb;
for k=1:maxit
neval++;
isgt = fgraterthan (xp,sliceheight);
if (all (isgt))
break
endif
lc = ! isgt & xp < xc;
uc = ! isgt & xp > xc;
lb(lc) = xp(lc);
ub(uc) = xp(uc);
xp = (ub - lb) .* rand (1, dim) + lb;
end
if (k == maxit)
warning ("slicesample: Step in exceeded maximum iterations");
endif
xc = xp;
if (i > burnin)
indx = (i - burnin) / thin;
if rem (indx, 1) == 0
smpl(indx, :) = xc;
end
end
end
neval = neval / (nsamples * thin + burnin);
endfunction
function y = rloge (x)
y = -inf (size (x));
xg0 = x > 0;
y(xg0) = log (x(xg0));
endfunction
%!demo
%! ## Define function to sample
%! d = 2;
%! mu = [-1; 2];
%! Sigma = rand (d);
%! Sigma = (Sigma + Sigma');
%! Sigma += eye (d)*abs (eigs (Sigma, 1, "sa")) * 1.1;
%! pdf = @(x)(2*pi)^(-d/2)*det(Sigma)^-.5*exp(-.5*sum((x.'-mu).*(Sigma\(x.'-mu)),1));
%! ##Inputs
%! start = ones (1,2);
%! nsamples = 500;
%! K = 500;
%! m = 10;
%! [smpl, accept]=slicesample (start, nsamples, "pdf", pdf, "burnin", K, "thin", m, "width", [20, 30]);
%! figure;
%! hold on;
%! plot (smpl(:,1), smpl(:,2), 'x');
%! [x, y] = meshgrid (linspace (-6,4), linspace(-3,7));
%! z = reshape (pdf ([x(:), y(:)]), size(x));
%! mesh (x, y, z, "facecolor", "None");
%! ## Using sample points to find the volume of half a sphere with radius of .5
%! f = @(x) ((.25-(x(:,1)+1).^2-(x(:,2)-2).^2).^.5.*(((x(:,1)+1).^2+(x(:,2)-2).^2)<.25)).';
%! int = mean (f (smpl) ./ pdf (smpl));
%! errest = std (f (smpl) ./ pdf (smpl)) / nsamples^.5;
%! trueerr = abs (2/3*pi*.25^(3/2)-int);
%! fprintf("Monte Carlo integral estimate int f(x) dx = %f\n", int);
%! fprintf("Monte Carlo integral error estimate %f\n", errest);
%! fprintf("The actual error %f\n", trueerr);
%! mesh (x,y,reshape (f([x(:), y(:)]), size(x)), "facecolor", "None");
%!demo
%! ##Integrate truncated normal distribution to find normilization constant
%! pdf = @(x) exp (-.5*x.^2)/(pi^.5*2^.5);
%! nsamples = 1e3;
%! [smpl,accept] = slicesample (1, nsamples, "pdf", pdf, "thin", 4);
%! f = @(x) exp (-.5 * x .^ 2) .* (x >= -2 & x <= 2);
%! x=linspace(-3,3,1000);
%! area(x,f(x));
%! xlabel ('x');
%! ylabel ('f(x)');
%! int = mean (f (smpl)./pdf(smpl));
%! errest = std (f (smpl)./pdf(smpl))/nsamples^.5;
%! trueerr = abs (erf (2^.5)*2^.5*pi^.5-int);
%! fprintf("Monte Carlo integral estimate int f(x) dx = %f\n", int);
%! fprintf("Monte Carlo integral error estimate %f\n", errest);
%! fprintf("The actual error %f\n", trueerr);
%!test
%! start = rand (1, 1);
%! nsamples = 1e3;
%! pdf = @(x) exp (-.5*(x-1).^2)/(2*pi)^.5;
%! [smpl, accept] = slicesample (start, nsamples, "pdf", pdf, "thin", 2, "burnin", 0, "width", 5);
%! assert (mean (smpl, 1), 1, .1);
%! assert (var (smpl, 1), 1, .1);
%!error slicesample ();
%!error slicesample (1);
%!error slicesample (1, 1);