########################################################################
##
## Copyright (C) 1993-2020 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{accept}] =} mhsample (@var{start}, @var{nsamples}, @var{property}, @var{value})
## Draws @var{nsamples} samples from a target stationary distribution @var{pdf}
## using Metropolis-Hastings algorithm.
##
## @var{smpl} is a @var{nsamples} by @var{dim} by @var{nchain} tensor of random
## values drawn from @var{pdf} where the rows are different random values, the
## columns correspond to the dimensions of @var{pdf} and the third dimension
## corresponds to different Markov chains.
##
## @var{accept} is a vector of the acceptance rate for each chain.
##
## @var{start} is a @var{nchain} by @var{dim} matrix of starting points for each
## Markov chain. Each row is the starting point of a different chain and each
## column corresponds to a different dimension.
##
## @var{nsamples} is the number of samples, the length of each Markov chain.
##
## Several property-value pairs must be specified, in pairs, They are:
##
## "pdf", @var{pdf} 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.
##
## "logpdf", @var{logpdf} 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. Only one of
## @var{pdf} or @var{logpdf} need to be specified.
##
## "proppdf" @var{proppdf} a function handle of the proposal distribution that
## is sampled from with @var{proprnd} to give the next point in the chain. The
## function should accept two inputs, the random variable and the current
## location each input should accept different locations in each row and each
## column corresponds to a different dimension.
##
## "logproppdf" @var{logproppdf} the log of "proppdf". Only one of
## @var{logproppdf} or @var{logproppdf} need to be specified.
##
## "proprnd" @var{proprnd} a function handle which generates random numbers from
## @var{proppdf} The function should accept different locations in each row and
## each column corresponds to a different dimension corresponding with the current
## location.
##
## "symmetric" @var{sym} true or false based on whether @var{proppdf} is a
## symmetric distribution. If true @var{proppdf} does not need to be given.
##
## "burnin" @var{k} the number of points to discard at the beginning, the default
## is 0.
##
## "thin" @var{m} omitts @var{m}-1 of every @var{m} points in the generated
## Markov chain. The default is 1.
##
## "nchain" @var{nchain} the number of Markov chains to generate. The default
## is 1.
##
## Example : Sampling from a normal distribution
##
## @example
## @group
## start = 1;
## nsamples = 1e3;
## pdf = @@(x) exp (-.5 * x .^ 2) / (pi ^ .5 * 2 ^ .5);
## proppdf = @@(x,y) 1 / 6;
## proprnd = @@(x) 6 * (rand (size (x)) - .5) + x;
## [smpl, accept] = mhsample (start, nsamples, "pdf", pdf, "proppdf", ...
## proppdf, "proprnd", proprnd, "thin", 4);
## histfit (smpl);
## @end group
## @end example
##
## @seealso{rand, slicesample}
## @end deftypefn
function [smpl, accept]=mhsample (start, nsamples, varargin)
if (nargin < 6)
print_usage ();
endif
sizestart = size (start);
pdf = [];
proppdf = [];
logpdf = [];
logproppdf = [];
proprnd = [];
sym = false;
K = 0;#burnin
m = 1;#thin
nchain = 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 ("mhsample: pdf must be a function handle");
endif
case "proppdf"
if (isa (varargin{k+1}, "function_handle"))
proppdf = varargin{k+1};
else
error ("mhsample: proppdf must be a function handle");
endif
case "logpdf"
if (isa (varargin{k+1}, "function_handle"))
pdf = varargin{k+1};
else
error ("mhsample: logpdf must be a function handle");
endif
case "logproppdf"
if (isa (varargin{k+1}, "function_handle"))
proppdf = varargin{k+1};
else
error ("mhsample: logproppdf must be a function handle");
endif
case "proprnd"
if (isa (varargin{k+1}, "function_handle"))
proprnd = varargin{k+1};
else
error ("mhsample: proprnd must be a function handle");
endif
case "symmetric"
if (isa (varargin{k+1}, "logical"))
sym = varargin{k+1};
else
error ("mhsample: sym must be true or false");
endif
case "burnin"
if (varargin{k+1}>=0)
K = varargin{k+1};
else
error ("mhsample: K must be greater than or equal to 0");
endif
case "thin"
if (varargin{k+1}>=1)
m = varargin{k+1};
else
error ("mhsample: m must be greater than or equal to 1");
endif
case "nchain"
if (varargin{k+1}>=1)
nchain = varargin{k+1};
else
error ("mhsample: nchain must be greater than or equal to 1");
endif
otherwise
warning (["mhsample: Ignoring unknown option " varargin{k}]);
endswitch
else
error (["mhsample: " 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 ("mhsample: pdf or logpdf must be input.");
endif
if (! isempty (proppdf) && isempty (logproppdf))
logproppdf = @(x,y) rloge (proppdf (x,y));
elseif (isempty(proppdf) && isempty(logproppdf) && ! sym)
error ("mhsample: proppdf or logproppdf must be input unless 'symetrical' is true.");
endif
if (! isa (proprnd, "function_handle"))
error ("mhsample: proprnd must be a function handle.");
endif
if length (sizestart)==2
sizestart = [sizestart 0];
end
smpl = zeros (nsamples, sizestart(2), nchain);
if (all (sizestart([1 3]) == [1 nchain]))#Could remove, not Matlab compatable but allows continuing chains
smpl(1,:,:) = start;
elseif (all (sizestart([1 3]) == [nchain 0]))
smpl(1,:,:) = permute (start, [3, 2, 1]);
elseif (all (sizestart([1 3]) == [1 0]))#Could remove, not Matlab compatable but allows all chains to start at the same location
smpl(1,:,:) = repmat (start,[1, 1, nchain]);
else
error ("mhsample: start must be a nchain by dim matrix.");
endif
cx = permute (smpl(1,:,:),[3, 2, 1]);
accept = zeros (nchain, 1);
i = 1;
rnd = log (rand (nchain, nsamples*m+K));
for k = 1:nsamples*m+K
canacc = rem (k-K,m) == 0;
px = proprnd (cx);
if (sym)
A = logpdf (px) - logpdf(cx);
else
A = (logpdf (px) + logproppdf (cx, px)) - (logpdf (cx) + logproppdf (px, cx));
endif
ac = rnd(:,k) < min (A, 0);
cx(ac,:) = px(ac,:);
accept(ac)++;
if (canacc)
smpl(i,:,:) = permute (cx,[3, 2, 1]);
end
if k > K && canacc
i++;
endif
endfor
accept ./= (nsamples * m + K);
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;
%! sym = true;
%! K = 500;
%! m = 10;
%! proprnd = @(x) (rand (size (x)) - .5) * 3 + x;
%! [smpl, accept]=mhsample (start, nsamples, "pdf", pdf, "proprnd", proprnd, "symmetric", sym, "burnin", K, "thin", m);
%! 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;
%! proprnd = @(x) (rand (size (x)) - .5)*3+x;
%! [smpl,accept] = mhsample(1,nsamples,"pdf",pdf,"proprnd", proprnd,"symmetric",true,"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
%! nchain=1e4;
%! start = rand (nchain, 1);
%! nsamples = 1e3;
%! pdf = @(x) exp (-.5*(x-1).^2)/(2*pi)^.5;
%! proppdf = @(x,y) 1/3;
%! proprnd = @(x) 3 * (rand (size (x)) - .5) + x;
%! [smpl, accept] = mhsample (start, nsamples, "pdf", pdf, "proppdf", ...
%! proppdf, "proprnd", proprnd, "thin", 2, "nchain", nchain, "burnin", 0);
%! assert (mean (mean (smpl, 1), 3), 1, .01);
%! assert (mean (var (smpl, 1), 3), 1, .01)
%!error mhsample ();
%!error mhsample (1);
%!error mhsample (1, 1);
%!error mhsample (1, 1, "pdf", @(x)x);
%!error mhsample (1, 1, "pdf", @(x)x, "proprnd", @(x)x+rand(size(x)));