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