From 2e16ac86c43559036ac124bd01623518a014362a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?=
<stephane.adjemian@univ-lemans.fr>
Date: Fri, 2 Dec 2011 12:12:46 +0100
Subject: [PATCH] Fixed bugs + cosmetic changes. (cherry picked from commit
a0e1d3b34fedf5d0c51c29d94ee8897028b9d78b)
Conflicts:
matlab/distributions/inverse_gamma_specification.m
---
matlab/distributions/ig1fun.m | 2 +-
.../inverse_gamma_specification.m | 29 ++++++++++++-------
2 files changed, 20 insertions(+), 11 deletions(-)
diff --git a/matlab/distributions/ig1fun.m b/matlab/distributions/ig1fun.m
index 25aa4624c..684ef3082 100644
--- a/matlab/distributions/ig1fun.m
+++ b/matlab/distributions/ig1fun.m
@@ -16,4 +16,4 @@ function err = ig1fun(nu,mu2,sigma2)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-err = 2*mu2*gamma(nu/2).^2-(sigma2+mu2)*(nu-2).*gamma((nu-1)/2).^2;
\ No newline at end of file
+err = log(2*mu2) - log((sigma2+mu2)*(nu-2)) + 2*( gammaln(nu/2)-gammaln((nu-1)/2) );
diff --git a/matlab/distributions/inverse_gamma_specification.m b/matlab/distributions/inverse_gamma_specification.m
index 867a42783..ea6a9981c 100644
--- a/matlab/distributions/inverse_gamma_specification.m
+++ b/matlab/distributions/inverse_gamma_specification.m
@@ -35,7 +35,12 @@ function [s,nu] = inverse_gamma_specification(mu,sigma,type,use_fzero_flag)
%! @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:
@@ -77,10 +82,9 @@ elseif type == 1; % Inverse Gamma 1
nu2 = 2*nu;
nu1 = 2;
err = ig1fun(nu,mu2,sigma2);
- err1 = ig1fun(nu1,mu2,sigma2);
err2 = ig1fun(nu2,mu2,sigma2);
- if sign(err2) % Too short interval.
- while nu2/nu<1000 % Shift the interval contaioning the root.
+ if err2>0 % Too short interval.
+ while nu2<500 % Shift the interval containing the root.
nu1 = nu2;
nu2 = nu2*1.01;
err2 = ig1fun(nu2,mu2,sigma2);
@@ -88,6 +92,9 @@ elseif type == 1; % Inverse Gamma 1
break
end
end
+ if err2>0
+ error('inverse_gamma_specification:: Failed in finding an interval containing a sign change! You should check that the prior variance is not too small compared to the prior mean...');
+ end
end
% Sove for nu using the secant method.
while abs(nu2-nu1) > 1e-14
@@ -125,14 +132,16 @@ else
end
%@test:1
-%$ %addpath ../matlab/distributions
-%$
-%$ % Define two dates
-%$ [s0,nu0] = inverse_gamma_specification(.75,.1,1,0);
-%$ [s1,nu1] = inverse_gamma_specification(.75,.1,1,1);
-%$
+%$ [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) = dyn_assert(s0,s1,1e-6);
%$ t(2) = dyn_assert(nu0,nu1,1e-6);
+%$ t(3) = isnan(s4);
+%$ t(4) = isnan(nu4);
+%$ t(5) = dyn_assert(s3,16.240907971002265,1e-6);;
+%$ t(6) = dyn_assert(nu3,30.368398202624046,1e-6);;
%$ T = all(t);
%@eof:1
\ No newline at end of file
--
GitLab