Draws a vector of candidate deep parameters in the joint prior density. B&K conditions have to be tested on each draw...
git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@999 ac1d8469-bf42-47a9-8791-bf33cf982152time-shift
parent
e66f3884db
commit
22e0f21220
|
@ -2,13 +2,16 @@ function pdraw = prior_draw(init,cc)
|
||||||
% Build one draw from the prior distribution.
|
% Build one draw from the prior distribution.
|
||||||
%
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
% o SampleSize [integer] Size of the sample to build
|
% o init [integer] scalar equal to 1 (first call) or 0.
|
||||||
|
% o cc [double] two columns matrix (same as in
|
||||||
|
% metropolis.m), constraints over the
|
||||||
|
% parameter space (upper and lower bounds).
|
||||||
%
|
%
|
||||||
% OUTPUTS
|
% OUTPUTS
|
||||||
% None.
|
% o pdraw [double] draw from the joint prior density.
|
||||||
%
|
%
|
||||||
% ALGORITHM
|
% ALGORITHM
|
||||||
% None.
|
% ...
|
||||||
%
|
%
|
||||||
% SPECIAL REQUIREMENTS
|
% SPECIAL REQUIREMENTS
|
||||||
% None.
|
% None.
|
||||||
|
@ -17,7 +20,7 @@ function pdraw = prior_draw(init,cc)
|
||||||
% part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
|
% part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
|
||||||
% Gnu Public License.
|
% Gnu Public License.
|
||||||
global M_ options_ estim_params_
|
global M_ options_ estim_params_
|
||||||
persistent fname npar bounds pshape pmean pstd a b p3 p4
|
persistent fname npar bounds pshape pmean pstd a b p3 p4 condition
|
||||||
|
|
||||||
if init
|
if init
|
||||||
nvx = estim_params_.nvx;
|
nvx = estim_params_.nvx;
|
||||||
|
@ -43,19 +46,37 @@ function pdraw = prior_draw(init,cc)
|
||||||
bounds = [-Inf Inf];
|
bounds = [-Inf Inf];
|
||||||
end
|
end
|
||||||
for i = 1:npar
|
for i = 1:npar
|
||||||
if pshape(i) == 3
|
switch pshape(i)
|
||||||
|
case 3% Gaussian prior
|
||||||
b(i) = pstd(i)^2/(pmean(i)-p3(i));
|
b(i) = pstd(i)^2/(pmean(i)-p3(i));
|
||||||
a(i) = (pmean(i)-p3(i))/b(i);
|
a(i) = (pmean(i)-p3(i))/b(i);
|
||||||
elseif pshape(i) == 1
|
case 1% Beta prior
|
||||||
mu = (p1(i)-p3(i))/(p4(i)-p3(i));
|
mu = (p1(i)-p3(i))/(p4(i)-p3(i));
|
||||||
stdd = p2(i)/(p4(i)-p3(i));
|
stdd = p2(i)/(p4(i)-p3(i));
|
||||||
a(i) = (1-mu)*mu^2/stdd^2 - mu;
|
a(i) = (1-mu)*mu^2/stdd^2 - mu;
|
||||||
b(i) = a*(1/mu - 1);
|
b(i) = a*(1/mu - 1);
|
||||||
elseif pshape(i) ==
|
case 2;%Gamma prior
|
||||||
|
mu = p1(i)-p3(i);
|
||||||
|
b(i) = p2(i)^2/mu;
|
||||||
|
a(i) = mu/b;
|
||||||
|
case {5,4,6}
|
||||||
|
% Nothing to do here
|
||||||
|
%
|
||||||
|
% 4: Inverse gamma, type 1, prior
|
||||||
|
% p2(i) = nu
|
||||||
|
% p1(i) = s
|
||||||
|
% 6: Inverse gamma, type 2, prior
|
||||||
|
% p2(i) = nu
|
||||||
|
% p1(i) = s
|
||||||
|
% 5: Uniform prior
|
||||||
|
% p3(i) and p4(i) are used.
|
||||||
|
otherwise
|
||||||
|
disp('prior_draw :: Error!')
|
||||||
|
disp('Unknown prior shape.')
|
||||||
|
return
|
||||||
end
|
end
|
||||||
pdraw = zeros(npar,1);
|
pdraw = zeros(npar,1);
|
||||||
end
|
end
|
||||||
|
|
||||||
condition = 1;
|
condition = 1;
|
||||||
pdraw = zeros(npar,1);
|
pdraw = zeros(npar,1);
|
||||||
return
|
return
|
||||||
|
@ -88,22 +109,29 @@ function pdraw = prior_draw(init,cc)
|
||||||
y2 = gamma_draw(b(i),1,0);
|
y2 = gamma_draw(b(i),1,0);
|
||||||
tmp = y1/(y1+y2);
|
tmp = y1/(y1+y2);
|
||||||
if tmp >= bounds(i,1) && tmp <= bounds(i,2)
|
if tmp >= bounds(i,1) && tmp <= bounds(i,2)
|
||||||
pdraw(i) = tmp;
|
pdraw(i) = pmean(i)+tmp*pstd(i);
|
||||||
break
|
break
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
case 4% INV-GAMMA1 distribution
|
case 4% INV-GAMMA1 distribution
|
||||||
|
while condition
|
||||||
case 6% INV-GAMMA2 distribution
|
tmp = sqrt(1/gamma_draw(p2(i)/2,1/p1(i),0));
|
||||||
|
if tmp >= bounds(i,1) && tmp <= bounds(i,2)
|
||||||
|
pdraw(i) = tmp;
|
||||||
|
break
|
||||||
otherwise
|
end
|
||||||
disp('prior_draw:: Error!')
|
end
|
||||||
disp('Unknown prior distribution.')
|
case 6% INV-GAMMA2 distribution
|
||||||
pdraw(i) = NaN;
|
while condition
|
||||||
|
tmp = 1/gamma_draw(p2(i)/2,1/p1(i),0);
|
||||||
|
if tmp >= bounds(i,1) && tmp <= bounds(i,2)
|
||||||
|
pdraw(i) = tmp;
|
||||||
|
break
|
||||||
|
end
|
||||||
|
end
|
||||||
|
otherwise
|
||||||
|
% Nothing to do here.
|
||||||
end
|
end
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
|
@ -116,13 +144,14 @@ function gamma_draw(a,b,c)
|
||||||
z = randn;
|
z = randn;
|
||||||
g = b*(z+sqrt(4*a-1))^2/4 + c;
|
g = b*(z+sqrt(4*a-1))^2/4 + c;
|
||||||
else
|
else
|
||||||
|
condi = 1
|
||||||
|
while condi
|
||||||
x = -1;
|
x = -1;
|
||||||
while x<0
|
while x<0
|
||||||
u1 = rand;
|
u1 = rand;
|
||||||
y = tan(pi*u1);
|
y = tan(pi*u1);
|
||||||
x = y*sqrt(2*a-1)+a-1;
|
x = y*sqrt(2*a-1)+a-1;
|
||||||
end
|
end
|
||||||
while condition
|
|
||||||
u2 = rand;
|
u2 = rand;
|
||||||
if log(u2) <= log(1+y^2)+(a-1)*log(x/(a-1))-y*sqrt(2*a-1);
|
if log(u2) <= log(1+y^2)+(a-1)*log(x/(a-1))-y*sqrt(2*a-1);
|
||||||
break
|
break
|
||||||
|
|
Loading…
Reference in New Issue