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