1function pdraw = prior_draw(BayesInfo, prior_trunc, uniform) % --*-- Unitary tests --*-- 2 3% This function generate one draw from the joint prior distribution and 4% allows sampling uniformly from the prior support (uniform==1 when called with init==1) 5% 6% INPUTS 7% o init [integer] scalar equal to: 8% 1: first call to set up persistent variables 9% describing the prior 10% 0: subsequent call to get prior 11% draw 12% o uniform [integer] scalar used in initialization (init=1), equal to: 13% 1: sample uniformly from prior 14% support (overwrites prior shape used for sampling within this function) 15% 0: sample from joint prior distribution 16% 17% OUTPUTS 18% o pdraw [double] 1*npar vector, draws from the joint prior density. 19% 20% 21% SPECIAL REQUIREMENTS 22% none 23% 24% NOTE 1. Input arguments 1 an 2 are only needed for initialization. 25% NOTE 2. A given draw from the joint prior distribution does not satisfy BK conditions a priori. 26% NOTE 3. This code relies on bayestopt_ as created in the base workspace 27% by the preprocessor (or as updated in subsequent pieces of code and handed to the base workspace) 28% 29% Copyright (C) 2006-2017 Dynare Team 30% 31% This file is part of Dynare. 32% 33% Dynare is free software: you can redistribute it and/or modify 34% it under the terms of the GNU General Public License as published by 35% the Free Software Foundation, either version 3 of the License, or 36% (at your option) any later version. 37% 38% Dynare is distributed in the hope that it will be useful, 39% but WITHOUT ANY WARRANTY; without even the implied warranty of 40% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 41% GNU General Public License for more details. 42% 43% You should have received a copy of the GNU General Public License 44% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 45 46persistent p6 p7 p3 p4 lb ub 47persistent uniform_index gaussian_index gamma_index beta_index inverse_gamma_1_index inverse_gamma_2_index weibull_index 48persistent uniform_draws gaussian_draws gamma_draws beta_draws inverse_gamma_1_draws inverse_gamma_2_draws weibull_draws 49 50if nargin>0 51 p6 = BayesInfo.p6; 52 p7 = BayesInfo.p7; 53 p3 = BayesInfo.p3; 54 p4 = BayesInfo.p4; 55 bounds = prior_bounds(BayesInfo, prior_trunc); 56 lb = bounds.lb; 57 ub = bounds.ub; 58 number_of_estimated_parameters = length(p6); 59 if nargin>2 && uniform 60 prior_shape = repmat(5,number_of_estimated_parameters,1); 61 else 62 prior_shape = BayesInfo.pshape; 63 end 64 beta_index = find(prior_shape==1); 65 if isempty(beta_index) 66 beta_draws = false; 67 else 68 beta_draws = true; 69 end 70 gamma_index = find(prior_shape==2); 71 if isempty(gamma_index) 72 gamma_draws = false; 73 else 74 gamma_draws = true; 75 end 76 gaussian_index = find(prior_shape==3); 77 if isempty(gaussian_index) 78 gaussian_draws = false; 79 else 80 gaussian_draws = true; 81 end 82 inverse_gamma_1_index = find(prior_shape==4); 83 if isempty(inverse_gamma_1_index) 84 inverse_gamma_1_draws = false; 85 else 86 inverse_gamma_1_draws = true; 87 end 88 uniform_index = find(prior_shape==5); 89 if isempty(uniform_index) 90 uniform_draws = false; 91 else 92 uniform_draws = true; 93 end 94 inverse_gamma_2_index = find(prior_shape==6); 95 if isempty(inverse_gamma_2_index) 96 inverse_gamma_2_draws = false; 97 else 98 inverse_gamma_2_draws = true; 99 end 100 weibull_index = find(prior_shape==8); 101 if isempty(weibull_index) 102 weibull_draws = false; 103 else 104 weibull_draws = true; 105 end 106 pdraw = NaN(number_of_estimated_parameters,1); 107 return 108end 109 110if uniform_draws 111 pdraw(uniform_index) = rand(length(uniform_index),1).*(p4(uniform_index)-p3(uniform_index)) + p3(uniform_index); 112 out_of_bound = find( (pdraw(uniform_index)'>ub(uniform_index)) | (pdraw(uniform_index)'<lb(uniform_index))); 113 while ~isempty(out_of_bound) 114 pdraw(uniform_index) = rand(length(uniform_index),1).*(p4(uniform_index)-p3(uniform_index)) + p3(uniform_index); 115 out_of_bound = find( (pdraw(uniform_index)'>ub(uniform_index)) | (pdraw(uniform_index)'<lb(uniform_index))); 116 end 117end 118 119if gaussian_draws 120 pdraw(gaussian_index) = randn(length(gaussian_index),1).*p7(gaussian_index) + p6(gaussian_index); 121 out_of_bound = find( (pdraw(gaussian_index)'>ub(gaussian_index)) | (pdraw(gaussian_index)'<lb(gaussian_index))); 122 while ~isempty(out_of_bound) 123 pdraw(gaussian_index(out_of_bound)) = randn(length(gaussian_index(out_of_bound)),1).*p7(gaussian_index(out_of_bound)) + p6(gaussian_index(out_of_bound)); 124 out_of_bound = find( (pdraw(gaussian_index)'>ub(gaussian_index)) | (pdraw(gaussian_index)'<lb(gaussian_index))); 125 end 126end 127 128if gamma_draws 129 pdraw(gamma_index) = gamrnd(p6(gamma_index),p7(gamma_index))+p3(gamma_index); 130 out_of_bound = find( (pdraw(gamma_index)'>ub(gamma_index)) | (pdraw(gamma_index)'<lb(gamma_index))); 131 while ~isempty(out_of_bound) 132 pdraw(gamma_index(out_of_bound)) = gamrnd(p6(gamma_index(out_of_bound)),p7(gamma_index(out_of_bound)))+p3(gamma_index(out_of_bound)); 133 out_of_bound = find( (pdraw(gamma_index)'>ub(gamma_index)) | (pdraw(gamma_index)'<lb(gamma_index))); 134 end 135end 136 137if beta_draws 138 pdraw(beta_index) = (p4(beta_index)-p3(beta_index)).*betarnd(p6(beta_index),p7(beta_index))+p3(beta_index); 139 out_of_bound = find( (pdraw(beta_index)'>ub(beta_index)) | (pdraw(beta_index)'<lb(beta_index))); 140 while ~isempty(out_of_bound) 141 pdraw(beta_index(out_of_bound)) = (p4(beta_index(out_of_bound))-p3(beta_index(out_of_bound))).*betarnd(p6(beta_index(out_of_bound)),p7(beta_index(out_of_bound)))+p3(beta_index(out_of_bound)); 142 out_of_bound = find( (pdraw(beta_index)'>ub(beta_index)) | (pdraw(beta_index)'<lb(beta_index))); 143 end 144end 145 146if inverse_gamma_1_draws 147 pdraw(inverse_gamma_1_index) = ... 148 sqrt(1./gamrnd(p7(inverse_gamma_1_index)/2,2./p6(inverse_gamma_1_index)))+p3(inverse_gamma_1_index); 149 out_of_bound = find( (pdraw(inverse_gamma_1_index)'>ub(inverse_gamma_1_index)) | (pdraw(inverse_gamma_1_index)'<lb(inverse_gamma_1_index))); 150 while ~isempty(out_of_bound) 151 pdraw(inverse_gamma_1_index(out_of_bound)) = ... 152 sqrt(1./gamrnd(p7(inverse_gamma_1_index(out_of_bound))/2,2./p6(inverse_gamma_1_index(out_of_bound))))+p3(inverse_gamma_1_index(out_of_bound)); 153 out_of_bound = find( (pdraw(inverse_gamma_1_index)'>ub(inverse_gamma_1_index)) | (pdraw(inverse_gamma_1_index)'<lb(inverse_gamma_1_index))); 154 end 155end 156 157if inverse_gamma_2_draws 158 pdraw(inverse_gamma_2_index) = ... 159 1./gamrnd(p7(inverse_gamma_2_index)/2,2./p6(inverse_gamma_2_index))+p3(inverse_gamma_2_index); 160 out_of_bound = find( (pdraw(inverse_gamma_2_index)'>ub(inverse_gamma_2_index)) | (pdraw(inverse_gamma_2_index)'<lb(inverse_gamma_2_index))); 161 while ~isempty(out_of_bound) 162 pdraw(inverse_gamma_2_index(out_of_bound)) = ... 163 1./gamrnd(p7(inverse_gamma_2_index(out_of_bound))/2,2./p6(inverse_gamma_2_index(out_of_bound)))+p3(inverse_gamma_2_index(out_of_bound)); 164 out_of_bound = find( (pdraw(inverse_gamma_2_index)'>ub(inverse_gamma_2_index)) | (pdraw(inverse_gamma_2_index)'<lb(inverse_gamma_2_index))); 165 end 166end 167 168if weibull_draws 169 pdraw(weibull_index) = wblrnd(p7(weibull_index), p6(weibull_index)) + p3(weibull_index); 170 out_of_bound = find( (pdraw(weibull_index)'>ub(weibull_index)) | (pdraw(weibull_index)'<lb(weibull_index))); 171 while ~isempty(out_of_bound) 172 pdraw(weibull_index(out_of_bound)) = wblrnd(p7(weibull_index(out_of_bound)),p6(weibull_index(out_of_bound)))+p3(weibull_index(out_of_bound)); 173 out_of_bound = find( (pdraw(weibull_index)'>ub(weibull_index)) | (pdraw(weibull_index)'<lb(weibull_index))); 174 end 175end 176 177%@test:1 178%$ % Fill global structures with required fields... 179%$ prior_trunc = 1e-10; 180%$ p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape 181%$ p1 = .4*ones(14,1); % Prior mean 182%$ p2 = .2*ones(14,1); % Prior std. 183%$ p3 = NaN(14,1); 184%$ p4 = NaN(14,1); 185%$ p5 = NaN(14,1); 186%$ p6 = NaN(14,1); 187%$ p7 = NaN(14,1); 188%$ 189%$ for i=1:14 190%$ switch p0(i) 191%$ case 1 192%$ % Beta distribution 193%$ p3(i) = 0; 194%$ p4(i) = 1; 195%$ [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i)); 196%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 1); 197%$ case 2 198%$ % Gamma distribution 199%$ p3(i) = 0; 200%$ p4(i) = Inf; 201%$ [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i)); 202%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 2); 203%$ case 3 204%$ % Normal distribution 205%$ p3(i) = -Inf; 206%$ p4(i) = Inf; 207%$ p6(i) = p1(i); 208%$ p7(i) = p2(i); 209%$ p5(i) = p1(i); 210%$ case 4 211%$ % Inverse Gamma (type I) distribution 212%$ p3(i) = 0; 213%$ p4(i) = Inf; 214%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false); 215%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 4); 216%$ case 5 217%$ % Uniform distribution 218%$ [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i)); 219%$ p3(i) = p6(i); 220%$ p4(i) = p7(i); 221%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 5); 222%$ case 6 223%$ % Inverse Gamma (type II) distribution 224%$ p3(i) = 0; 225%$ p4(i) = Inf; 226%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false); 227%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 6); 228%$ case 8 229%$ % Weibull distribution 230%$ p3(i) = 0; 231%$ p4(i) = Inf; 232%$ [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i)); 233%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 8); 234%$ otherwise 235%$ error('This density is not implemented!') 236%$ end 237%$ end 238%$ 239%$ BayesInfo.pshape = p0; 240%$ BayesInfo.p1 = p1; 241%$ BayesInfo.p2 = p2; 242%$ BayesInfo.p3 = p3; 243%$ BayesInfo.p4 = p4; 244%$ BayesInfo.p5 = p5; 245%$ BayesInfo.p6 = p6; 246%$ BayesInfo.p7 = p7; 247%$ 248%$ ndraws = 1e5; 249%$ m0 = BayesInfo.p1; %zeros(14,1); 250%$ v0 = diag(BayesInfo.p2.^2); %zeros(14); 251%$ 252%$ % Call the tested routine 253%$ try 254%$ % Initialization of prior_draws. 255%$ prior_draw(BayesInfo, prior_trunc, false); 256%$ % Do simulations in a loop and estimate recursively the mean and teh variance. 257%$ for i = 1:ndraws 258%$ draw = transpose(prior_draw()); 259%$ m1 = m0 + (draw-m0)/i; 260%$ m2 = m1*m1'; 261%$ v0 = v0 + ((draw*draw'-m2-v0) + (i-1)*(m0*m0'-m2'))/i; 262%$ m0 = m1; 263%$ end 264%$ t(1) = true; 265%$ catch 266%$ t(1) = false; 267%$ end 268%$ 269%$ if t(1) 270%$ t(2) = all(abs(m0-BayesInfo.p1)<3e-3); 271%$ t(3) = all(all(abs(v0-diag(BayesInfo.p2.^2))<3e-3)); 272%$ end 273%$ T = all(t); 274%@eof:1 275