From cfefa6951c0b43c4513e9417242252f9680b6bc2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Dec 2015 17:51:54 +0100 Subject: [PATCH] Changed inputs of inverse_gamma_specification routine and added unit tests. --- .../inverse_gamma_specification.m | 185 ++++++++++++------ matlab/set_prior.m | 14 +- 2 files changed, 123 insertions(+), 76 deletions(-) diff --git a/matlab/distributions/inverse_gamma_specification.m b/matlab/distributions/inverse_gamma_specification.m index 20b1b9b35..f6634c186 100644 --- a/matlab/distributions/inverse_gamma_specification.m +++ b/matlab/distributions/inverse_gamma_specification.m @@ -1,50 +1,25 @@ -function [s,nu] = inverse_gamma_specification(mu,sigma,type,use_fzero_flag) +function [s,nu] = inverse_gamma_specification(mu, sigma2, type, use_fzero_flag, name) % --*-- Unitary tests --*-- + % Computes the inverse Gamma hyperparameters from the prior mean and standard deviation. +% +% INPUTS +% - mu [double] scalar, prior mean. +% - sigma2 [double] positive scalar, prior variance. +% - type [integer] scalar equal to 1 or 2, type of the inverse gamma distribution +% - use_fzero_flag [logical] scalar, Use (matlab/octave's implementation of) fzero to solve for nu if true, use +% dynare's implementation of the secant method otherwise. +% - name [string] name of the parameter or random variable. +% +% OUTPUS +% - s [double] scalar, first hyperparameter. +% - nu [double] scalar, second hyperparameter. +% +% REMARK +% The call to the matlab's implementation of the secant method is here for testing purpose and should not be used. This routine fails +% more often in finding an interval for nu containing a signe change because it expands the interval on both sides and eventually +% violates the condition nu>2. -%@info: -%! @deftypefn {Function File} {[@var{s}, @var{nu} ]=} colon (@var{mu}, @var{sigma}, @var{type}, @var{use_fzero_flag}) -%! @anchor{distributions/inverse_gamma_specification} -%! @sp 1 -%! Computes the inverse Gamma (type 1 or 2) hyperparameters from the prior mean (@var{mu}) and standard deviation (@var{sigma}). -%! @sp 2 -%! @strong{Inputs} -%! @sp 1 -%! @table @ @var -%! @item mu -%! Double scalar, prior mean. -%! @item sigma -%! Positive double scalar, prior standard deviation. -%! @item type -%! Integer scalar equal to one or two, type of the Inverse-Gamma distribution. -%! @item use_fzero_flag -%! Integer scalar equal to 0 (default) or 1. Use (matlab/octave's implementation of) fzero to solve for @var{nu} if equal to 1, use -%! dynare's implementation of the secant method otherwise. -%! @end table -%! @sp 1 -%! @strong{Outputs} -%! @sp 1 -%! @table @ @var -%! @item s -%! Positive double scalar (greater than two), first hypermarameter of the Inverse-Gamma prior. -%! @item nu -%! Positive double scala, second hypermarameter of the Inverse-Gamma prior. -%! @end table -%! @sp 2 -%! @strong{This function is called by:} -%! @sp 1 -%! @ref{set_prior} -%! @sp 2 -%! @strong{This function calls:} -%! @sp 2 -%! @strong{Remark:} -%! The call to the matlab's implementation of the secant method is here for testing purpose and should not be used. This routine fails -%! more often in finding an interval for nu containing a signe change because it expands the interval on both sides and eventually -%! violates the condition nu>2. -%! -%! @end deftypefn -%@eod: - -% Copyright (C) 2003-2012 Dynare Team +% Copyright (C) 2003-2015 Dynare Team % % This file is part of Dynare. % @@ -61,16 +36,52 @@ function [s,nu] = inverse_gamma_specification(mu,sigma,type,use_fzero_flag) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -check_solution_flag = 1; +if nargin<3 + error('At least three input arguments are required!') +end + +if ~isnumeric(mu) || ~isscalar(mu) || ~isreal(mu) + error('First input argument must be a real positive scalar!') +end + +if ~isnumeric(sigma2) || ~isscalar(sigma2) || ~isreal(sigma2) || sigma2<=0 + error('Second input argument must be a real positive scalar!') +end + +if ~isnumeric(mu) || ~isscalar(mu) || ~ismember(type, [1, 2]) + error('Third input argument must be equal to 1 or 2!') +end + +if nargin==3 || isempty(use_fzero_flag) + use_fzero_flag = false; +else + if ~iscalar(use_fzero_flag) || ~islogical(use_fzero_flag) + error('Fourth input argument must be a scalar logical!') + end +end + +if nargin>4 && (~ischar(name) || size(name, 1)>1) + error('Fifth input argument must be a string!') +else + name = ''; +end + +if ~isempty(name) + name = sprintf(' for %s ', name); +else + name = ' '; +end + +if mu<=0 + error('The prior mean%smust be above the lower bound of the Inverse Gamma (type %d) prior distribution!', name, type); +end + +check_solution_flag = true; s = []; nu = []; -if nargin==3 - use_fzero_flag = 0; -end - -sigma2 = sigma^2; -mu2 = mu^2; +sigma = sqrt(sigma2); +mu2 = mu*mu; if type == 2; % Inverse Gamma 2 nu = 2*(2+mu2/sigma2); @@ -120,7 +131,7 @@ elseif type == 1; % Inverse Gamma 1 if abs(log(mu)-log(sqrt(s/2))-gammaln((nu-1)/2)+gammaln(nu/2))>1e-7 error('inverse_gamma_specification:: Failed in solving for the hyperparameters!'); end - if abs(sigma-sqrt(s/(nu-2)-mu^2))>1e-7 + if abs(sigma-sqrt(s/(nu-2)-mu*mu))>1e-7 error('inverse_gamma_specification:: Failed in solving for the hyperparameters!'); end end @@ -133,17 +144,61 @@ else end %@test:1 +%$ try +%$ [s, nu] = inverse_gamma_specification(.5, .05, 1); +%$ t(1) = true; +%$ catch +%$ t(1) = false; +%$ end %$ -%$ [s0,nu0] = inverse_gamma_specification(.75,.2,1,0); -%$ [s1,nu1] = inverse_gamma_specification(.75,.2,1,1); -%$ [s3,nu3] = inverse_gamma_specification(.75,.1,1,0); -%$ [s4,nu4] = inverse_gamma_specification(.75,.1,1,1); -%$ % Check the results. -%$ t(1) = dassert(s0,s1,1e-6); -%$ t(2) = dassert(nu0,nu1,1e-6); -%$ t(3) = isnan(s4); -%$ t(4) = isnan(nu4); -%$ t(5) = dassert(s3,16.240907971002265,1e-6);; -%$ t(6) = dassert(nu3,30.368398202624046,1e-6);; +%$ if t(1) +%$ t(2) = abs(0.5-sqrt(.5*s)*gamma(.5*(nu-1))/gamma(.5*nu))<1e-12; +%$ t(3) = abs(0.05-s/(nu-2)+.5^2)<1e-12; +%$ end %$ T = all(t); -%@eof:1 \ No newline at end of file +%@eof:1 + +%@test:2 +%$ try +%$ [s, nu] = inverse_gamma_specification(.5, .05, 2); +%$ t(1) = true; +%$ catch +%$ t(1) = false; +%$ end +%$ +%$ if t(1) +%$ t(2) = abs(0.5-s/(nu-2))<1e-12; +%$ t(3) = abs(0.05-2*.5^2/(nu-4))<1e-12; +%$ end +%$ T = all(t); +%@eof:2 + +%@test:3 +%$ try +%$ [s, nu] = inverse_gamma_specification(.5, Inf, 1); +%$ t(1) = true; +%$ catch +%$ t(1) = false; +%$ end +%$ +%$ if t(1) +%$ t(2) = abs(0.5-sqrt(.5*s)*gamma(.5*(nu-1))/gamma(.5*nu))<1e-12; +%$ t(3) = isequal(nu, 2); %abs(0.05-2*.5^2/(nu-4))<1e-12; +%$ end +%$ T = all(t); +%@eof:3 + +%@test:4 +%$ try +%$ [s, nu] = inverse_gamma_specification(.5, Inf, 2); +%$ t(1) = true; +%$ catch +%$ t(1) = false; +%$ end +%$ +%$ if t(1) +%$ t(2) = abs(0.5-s/(nu-2))<1e-12; +%$ t(3) = isequal(nu, 4); +%$ end +%$ T = all(t); +%@eof:4 \ No newline at end of file diff --git a/matlab/set_prior.m b/matlab/set_prior.m index 3d3aa458f..698e0a5dd 100644 --- a/matlab/set_prior.m +++ b/matlab/set_prior.m @@ -206,12 +206,8 @@ k2 = find(isnan(bayestopt_.p4(k))); bayestopt_.p3(k(k1)) = zeros(length(k1),1); bayestopt_.p4(k(k2)) = Inf(length(k2),1); for i=1:length(k) - if (bayestopt_.p1(k(i))bayestopt_.p4(k(i))) - error(['The prior mean of ' bayestopt_.name{k(i)} ' has to be above the lower (' num2str(bayestopt_.p3(k(i))) ') bound of the Inverse Gamma prior density!']); - end - [bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = ... - inverse_gamma_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)),bayestopt_.p2(k(i)),1,0) ; - bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 4) ; + [bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = inverse_gamma_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)), bayestopt_.p2(k(i)), 1, false, bayestopt_.name{k(i)}); + bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 4); end % uniform distribution @@ -231,11 +227,7 @@ k2 = find(isnan(bayestopt_.p4(k))); bayestopt_.p3(k(k1)) = zeros(length(k1),1); bayestopt_.p4(k(k2)) = Inf(length(k2),1); for i=1:length(k) - if (bayestopt_.p1(k(i))bayestopt_.p4(k(i))) - error(['The prior mean of ' bayestopt_.name{k(i)} ' has to be above the lower (' num2str(bayestopt_.p3(k(i))) ') bound of the Inverse Gamma II prior density!']); - end - [bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = ... - inverse_gamma_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)),bayestopt_.p2(k(i)),2,0); + [bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = inverse_gamma_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)), bayestopt_.p2(k(i)), 2, false, bayestopt_.name{k(i)}); bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 6) ; end