diff --git a/matlab/prior_draw.m b/matlab/prior_draw.m index 9e49f8d3685d140c73d57ff6444c188ef6e92fa9..8cebcdfe479ab59dbdb08dc2806808820362237e 100644 --- a/matlab/prior_draw.m +++ b/matlab/prior_draw.m @@ -175,368 +175,100 @@ if weibull_draws end %@test:1 -%$ %% Initialize required structures -%$ options_.prior_trunc=0; -%$ options_.initialize_estimated_parameters_with_the_prior_mode=0; -%$ -%$ M_.dname='test'; -%$ M_.param_names = 'alp'; -%$ ndraws=100000; -%$ global estim_params_ -%$ estim_params_.var_exo = []; -%$ estim_params_.var_endo = []; -%$ estim_params_.corrx = []; -%$ estim_params_.corrn = []; -%$ estim_params_.param_vals = []; -%$ estim_params_.param_vals = [1, NaN, (-Inf), Inf, 1, 0.356, 0.02, NaN, NaN, NaN ]; -%$ -%$ %% beta -%$ estim_params_.param_vals(1,3)= -Inf; % LB -%$ estim_params_.param_vals(1,4)= +Inf; % UB -%$ estim_params_.param_vals(1,5)= 1; % Shape -%$ estim_params_.param_vals(1,6)=0.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=NaN; -%$ estim_params_.param_vals(1,9)=NaN; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_) -%$ -%$ pdraw = prior_draw(1,0); pdraw -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec<0) || any(pdraw_vec>1) -%$ error('Beta prior wrong') -%$ end -%$ -%$ -%$ %% Gamma -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 2;%Shape -%$ estim_params_.param_vals(1,6)=0.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=NaN; -%$ estim_params_.param_vals(1,9)=NaN; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec<0) -%$ error('Gamma prior wrong') -%$ end -%$ -%$ %% Normal -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 3;%Shape -%$ estim_params_.param_vals(1,6)=0.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=NaN; -%$ estim_params_.param_vals(1,9)=NaN; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 -%$ error('Normal prior wrong') -%$ end -%$ -%$ %% inverse gamma distribution (type 1) -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 4;%Shape -%$ estim_params_.param_vals(1,6)=0.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=NaN; -%$ estim_params_.param_vals(1,9)=NaN; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec<0) -%$ error('inverse gamma distribution (type 1) prior wrong') -%$ end -%$ -%$ %% Uniform -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 5;%Shape -%$ estim_params_.param_vals(1,6)=0.5; -%$ estim_params_.param_vals(1,7)=sqrt(12)^(-1)*(1-0); -%$ estim_params_.param_vals(1,8)=NaN; -%$ estim_params_.param_vals(1,9)=NaN; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-2 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-3 || any(pdraw_vec<0) || any(pdraw_vec>1) -%$ error('Uniform prior wrong') -%$ end -%$ -%$ %% inverse gamma distribution (type 2) -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 6;%Shape -%$ estim_params_.param_vals(1,6)=0.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=NaN; -%$ estim_params_.param_vals(1,9)=NaN; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec<0) -%$ error('inverse gamma distribution (type 2) prior wrong') -%$ end -%$ -%$ -%$ %%%%%%%%%%%%%%%%%%%%%% Generalized distributions %%%%%%%%%%%%%%%%%%%%% -%$ -%$ %% beta -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 1;%Shape -%$ estim_params_.param_vals(1,6)=1.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=1; -%$ estim_params_.param_vals(1,9)=3; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec<estim_params_.param_vals(1,3)) || any(pdraw_vec>estim_params_.param_vals(1,4)) -%$ error('Beta prior wrong') -%$ end -%$ -%$ %% Gamma -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 2;%Shape -%$ estim_params_.param_vals(1,6)=1.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=1; -%$ estim_params_.param_vals(1,9)=NaN; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec<estim_params_.param_vals(1,8)) -%$ error('Gamma prior wrong') -%$ end -%$ -%$ %% inverse gamma distribution (type 1) -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 4;%Shape -%$ estim_params_.param_vals(1,6)=1.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=1; -%$ estim_params_.param_vals(1,9)=NaN; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec<estim_params_.param_vals(1,8)) -%$ error('inverse gamma distribution (type 1) prior wrong') -%$ end -%$ -%$ %% Uniform -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 5;%Shape -%$ estim_params_.param_vals(1,6)=1.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=NaN; -%$ estim_params_.param_vals(1,9)=NaN; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec<estim_params_.param_vals(1,3)) || any(pdraw_vec>estim_params_.param_vals(1,4)) -%$ error('Uniform prior wrong') -%$ end -%$ -%$ %% inverse gamma distribution (type 2) -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 6;%Shape -%$ estim_params_.param_vals(1,6)=1.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=1; -%$ estim_params_.param_vals(1,9)=NaN; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec<estim_params_.param_vals(1,8)) -%$ error('inverse gamma distribution (type 2) prior wrong') -%$ end -%$ -%$ %%%%%%%%%%%% With prior truncation -%$ options_.prior_trunc=.4; -%$ -%$ %% beta -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 1;%Shape -%$ estim_params_.param_vals(1,6)=1.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=1; -%$ estim_params_.param_vals(1,9)=3; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ bounds = prior_bounds(bayestopt_,options_.prior_trunc)'; -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>5e-3 || any(pdraw_vec<bounds.lb) || any(pdraw_vec>bounds.ub) -%$ error('Beta prior wrong') -%$ end -%$ -%$ %% Gamma -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 2;%Shape -%$ estim_params_.param_vals(1,6)=1.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=1; -%$ estim_params_.param_vals(1,9)=NaN; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ bounds = prior_bounds(bayestopt_,options_.prior_trunc)'; -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>5e-3 || any(pdraw_vec<bounds.lb) || any(pdraw_vec>bounds.ub) -%$ error('Gamma prior wrong') -%$ end -%$ -%$ %% inverse gamma distribution (type 1) -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 4;%Shape -%$ estim_params_.param_vals(1,6)=1.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=1; -%$ estim_params_.param_vals(1,9)=NaN; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ bounds = prior_bounds(bayestopt_,options_.prior_trunc)'; -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>5e-3 || any(pdraw_vec<bounds.lb) || any(pdraw_vec>bounds.ub) -%$ error('inverse gamma distribution (type 1) prior wrong') -%$ end -%$ -%$ %% Uniform -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 5;%Shape -%$ estim_params_.param_vals(1,6)=1.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=NaN; -%$ estim_params_.param_vals(1,9)=NaN; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ bounds = prior_bounds(bayestopt_,options_.prior_trunc)'; -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>5e-3 || any(pdraw_vec<bounds.lb) || any(pdraw_vec>bounds.ub) -%$ error('Uniform prior wrong') -%$ end -%$ -%$ -%$ %% inverse gamma distribution (type 2) -%$ estim_params_.param_vals(1,3)= -Inf;%LB -%$ estim_params_.param_vals(1,4)= +Inf;%UB -%$ estim_params_.param_vals(1,5)= 6;%Shape -%$ estim_params_.param_vals(1,6)=1.5; -%$ estim_params_.param_vals(1,7)=0.01; -%$ estim_params_.param_vals(1,8)=1; -%$ estim_params_.param_vals(1,9)=NaN; -%$ -%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); -%$ bounds = prior_bounds(bayestopt_,options_.prior_trunc)'; -%$ -%$ pdraw = prior_draw(1,0); -%$ pdraw_vec=NaN(ndraws,1); -%$ for ii=1:ndraws -%$ pdraw_vec(ii)=prior_draw(0,0); -%$ end -%$ -%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>5e-3 || any(pdraw_vec<bounds.lb) || any(pdraw_vec>bounds.ub) -%$ error('inverse gamma distribution (type 2) prior wrong') -%$ end -%$ -%@eof:1 +%$ % Fill global structures with required fields... +%$ prior_trunc = 1e-10; +%$ p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape +%$ p1 = .4*ones(14,1); % Prior mean +%$ p2 = .2*ones(14,1); % Prior std. +%$ p3 = NaN(14,1); +%$ p4 = NaN(14,1); +%$ p5 = NaN(14,1); +%$ p6 = NaN(14,1); +%$ p7 = NaN(14,1); +%$ +%$ for i=1:14 +%$ switch p0(i) +%$ case 1 +%$ % Beta distribution +%$ p3(i) = 0; +%$ p4(i) = 1; +%$ [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i)); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 1); +%$ case 2 +%$ % Gamma distribution +%$ p3(i) = 0; +%$ p4(i) = Inf; +%$ [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i)); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 2); +%$ case 3 +%$ % Normal distribution +%$ p3(i) = -Inf; +%$ p4(i) = Inf; +%$ p6(i) = p1(i); +%$ p7(i) = p2(i); +%$ p5(i) = p1(i); +%$ case 4 +%$ % Inverse Gamma (type I) distribution +%$ p3(i) = 0; +%$ p4(i) = Inf; +%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 4); +%$ case 5 +%$ % Uniform distribution +%$ [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i)); +%$ p3(i) = p6(i); +%$ p4(i) = p7(i); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 5); +%$ case 6 +%$ % Inverse Gamma (type II) distribution +%$ p3(i) = 0; +%$ p4(i) = Inf; +%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 6); +%$ case 8 +%$ % Weibull distribution +%$ p3(i) = 0; +%$ p4(i) = Inf; +%$ [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i)); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 8); +%$ otherwise +%$ error('This density is not implemented!') +%$ end +%$ end +%$ +%$ BayesInfo.pshape = p0; +%$ BayesInfo.p1 = p1; +%$ BayesInfo.p2 = p2; +%$ BayesInfo.p3 = p3; +%$ BayesInfo.p4 = p4; +%$ BayesInfo.p5 = p5; +%$ BayesInfo.p6 = p6; +%$ BayesInfo.p7 = p7; +%$ +%$ ndraws = 1e5; +%$ m0 = BayesInfo.p1; %zeros(14,1); +%$ v0 = diag(BayesInfo.p2.^2); %zeros(14); +%$ +%$ % Call the tested routine +%$ try +%$ % Initialization of prior_draws. +%$ prior_draw(BayesInfo, prior_trunc, false); +%$ % Do simulations in a loop and estimate recursively the mean and teh variance. +%$ for i = 1:ndraws +%$ draw = transpose(prior_draw()); +%$ m1 = m0 + (draw-m0)/i; +%$ m2 = m1*m1'; +%$ v0 = v0 + ((draw*draw'-m2-v0) + (i-1)*(m0*m0'-m2'))/i; +%$ m0 = m1; +%$ end +%$ t(1) = true; +%$ catch +%$ t(1) = false; +%$ end +%$ +%$ if t(1) +%$ t(2) = all(abs(m0-BayesInfo.p1)<3e-3); +%$ t(3) = all(all(abs(v0-diag(BayesInfo.p2.^2))<2e-3)); +%$ end +%$ T = all(t); +%@eof:1 \ No newline at end of file