1function [theta, fxsim, neval] = rotated_slice_sampler(objective_function,theta,thetaprior,sampler_options,varargin) 2% ---------------------------------------------------------- 3% ROTATED SLICE SAMPLER - with stepping out (Neal, 2003) 4% extension of the orthogonal univarite sampler (slice_sampler.m) 5% copyright M. Ratto (European Commission) 6% 7% objective_function(theta,varargin): -log of any unnormalized pdf 8% with varargin (optional) a vector of auxiliaty parameters 9% to be passed to f( ). 10% ---------------------------------------------------------- 11% 12% INPUTS 13% objective_function: objective function (expressed as minus the log of a density) 14% theta: last value of theta 15% thetaprior: bounds of the theta space 16% sampler_options: posterior sampler options 17% varargin: optional input arguments to objective function 18% 19% OUTPUTS 20% theta: new theta sample 21% fxsim: value of the objective function for the new sample 22% neval: number of function evaluations 23% 24% SPECIAL REQUIREMENTS 25% none 26 27% Copyright (C) 2015-2017 Dynare Team 28% 29% This file is part of Dynare. 30% 31% Dynare is free software: you can redistribute it and/or modify 32% it under the terms of the GNU General Public License as published by 33% the Free Software Foundation, either version 3 of the License, or 34% (at your option) any later version. 35% 36% Dynare is distributed in the hope that it will be useful, 37% but WITHOUT ANY WARRANTY; without even the implied warranty of 38% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 39% GNU General Public License for more details. 40% 41% You should have received a copy of the GNU General Public License 42% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 43 44theta=theta(:); 45npar = length(theta); 46neval = zeros(npar,1); 47W1=[]; 48if isfield(sampler_options,'WR') 49 W1 = sampler_options.WR; 50end 51if ~isempty(sampler_options.mode) 52 mm = sampler_options.mode; 53 n = length(mm); 54 for j=1:n 55 distance(j)=sqrt(sum((theta-mm(j).m).^2)); 56 end 57 [m, im] = min(distance); 58 59 r=im; 60 V1 = mm(r).m; 61 jj=0; 62 for j=1:n 63 if j~=r 64 jj=jj+1; 65 tmp=mm(j).m-mm(r).m; 66 %tmp=mm(j).m-theta; 67 V1(:,jj)=tmp/norm(tmp); 68 end 69 end 70 resul=randperm(n-1,n-1); 71 V1 = V1(:,resul); 72 73 %V1 = V1(:, randperm(n-1)); 74 % %d = chol(mm(r).invhess); 75 % %V1 = transpose(feval(sampler_options.proposal_distribution, transpose(mm(r).m), d, npar)); 76 % 77 % V1=eye(npar); 78 % V1=V1(:,randperm(npar)); 79 % for j=1:2, 80 % V1(:,j)=mm(r(j)).m-theta; 81 % V1(:,j)=V1(:,j)/norm(V1(:,j)); 82 % end 83 % % Gram-Schmidt 84 % for j=2:npar, 85 % for k=1:j-1, 86 % V1(:,j)=V1(:,j)-V1(:,k)'*V1(:,j)*V1(:,k); 87 % end 88 % V1(:,j)=V1(:,j)/norm(V1(:,j)); 89 % end 90 % for j=1:n, 91 % distance(j)=sqrt(sum((theta-mm(j).m).^2)); 92 % end 93 % [m, im] = min(distance); 94 % if im==r, 95 % fxsim=[]; 96 % return, 97 % else 98 % theta1=theta; 99 % end 100else 101 V1 = sampler_options.V1; 102end 103npar=size(V1,2); 104 105for it=1:npar 106 theta0 = theta; 107 neval(it) = 0; 108 xold = 0; 109 % XLB = thetaprior(3); 110 % XUB = thetaprior(4); 111 tb=sort([(thetaprior(:,1)-theta)./V1(:,it) (thetaprior(:,2)-theta)./V1(:,it)],2); 112 XLB=max(tb(:,1)); 113 XUB=min(tb(:,2)); 114 if isempty(W1) 115 W = (XUB-XLB); %*0.8; 116 else 117 W = W1(it); 118 end 119 120 % ------------------------------------------------------- 121 % 1. DRAW Z = ln[f(X0)] - EXP(1) where EXP(1)=-ln(U(0,1)) 122 % THIS DEFINES THE SLICE S={x: z < ln(f(x))} 123 % ------------------------------------------------------- 124 125 fxold = -feval(objective_function,theta,varargin{:}); 126 %I have to be sure that the rotation is for L,R or for Fxold, theta(it) 127 neval(it) = neval(it) + 1; 128 Z = fxold + log(rand(1,1)); 129 % ------------------------------------------------------------- 130 % 2. FIND I=(L,R) AROUND X0 THAT CONTAINS S AS MUCH AS POSSIBLE 131 % STEPPING-OUT PROCEDURE 132 % ------------------------------------------------------------- 133 u = rand(1,1); 134 L = max(XLB,xold-W*u); 135 R = min(XUB,L+W); 136 137 %[L R]=slice_rotation(L, R, alpha); 138 while(L > XLB) 139 xsim = L; 140 theta = theta0+xsim*V1(:,it); 141 fxl = -feval(objective_function,theta,varargin{:}); 142 neval(it) = neval(it) + 1; 143 if (fxl <= Z) 144 break 145 end 146 L = max(XLB,L-W); 147 end 148 while(R < XUB) 149 xsim = R; 150 theta = theta0+xsim*V1(:,it); 151 fxr = -feval(objective_function,theta,varargin{:}); 152 neval(it) = neval(it) + 1; 153 if (fxr <= Z) 154 break 155 end 156 R = min(XUB,R+W); 157 end 158 % ------------------------------------------------------ 159 % 3. SAMPLING FROM THE SET A = (I INTERSECT S) = (LA,RA) 160 % ------------------------------------------------------ 161 fxsim = Z-1; 162 while (fxsim < Z) 163 u = rand(1,1); 164 xsim = L + u*(R - L); 165 theta = theta0+xsim*V1(:,it); 166 fxsim = -feval(objective_function,theta,varargin{:}); 167 neval(it) = neval(it) + 1; 168 if (xsim > xold) 169 R = xsim; 170 else 171 L = xsim; 172 end 173 end 174end 175 176% if ~isempty(sampler_options.mode), 177% dist1=sqrt(sum((theta-mm(r).m).^2)); 178% if dist1>distance(r), 179% theta=theta1; 180% fxsim=[]; 181% end 182% end 183end 184