######################################################################## ## ## 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) %! 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) heaviside(x+2).*exp(-.5*x.^2).*(1-heaviside(x-2)); %! int = mean (f (smpl)./pdf(smpl)) %! errest = std (f (smpl)./pdf(smpl))/nsamples^.5 %! trueerr = abs (erf (2^.5)*2^.5*pi^.5-int) %!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)));