Computes the inverse Gamma hyperparameters from the prior mean and standard deviation.
0001 function [s,nu] = inverse_gamma_specification(mu,sigma,type,use_fzero_flag) 0002 % Computes the inverse Gamma hyperparameters from the prior mean and standard deviation. 0003 0004 %@info: 0005 %! @deftypefn {Function File} {[@var{s}, @var{nu} ]=} colon (@var{mu}, @var{sigma}, @var{type}, @var{use_fzero_flag}) 0006 %! @anchor{distributions/inverse_gamma_specification} 0007 %! @sp 1 0008 %! Computes the inverse Gamma (type 1 or 2) hyperparameters from the prior mean (@var{mu}) and standard deviation (@var{sigma}). 0009 %! @sp 2 0010 %! @strong{Inputs} 0011 %! @sp 1 0012 %! @table @ @var 0013 %! @item mu 0014 %! Double scalar, prior mean. 0015 %! @item sigma 0016 %! Positive double scalar, prior standard deviation. 0017 %! @item type 0018 %! Integer scalar equal to one or two, type of the Inverse-Gamma distribution. 0019 %! @item use_fzero_flag 0020 %! 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 0021 %! dynare's implementation of the secant method otherwise. 0022 %! @end table 0023 %! @sp 1 0024 %! @strong{Outputs} 0025 %! @sp 1 0026 %! @table @ @var 0027 %! @item s 0028 %! Positive double scalar (greater than two), first hypermarameter of the Inverse-Gamma prior. 0029 %! @item nu 0030 %! Positive double scala, second hypermarameter of the Inverse-Gamma prior. 0031 %! @end table 0032 %! @sp 2 0033 %! @strong{This function is called by:} 0034 %! @sp 1 0035 %! @ref{set_prior} 0036 %! @sp 2 0037 %! @strong{This function calls:} 0038 %! @sp 2 0039 %! @strong{Remark:} 0040 %! The call to the matlab's implementation of the secant method is here for testing purpose and should not be used. This routine fails 0041 %! more often in finding an interval for nu containing a signe change because it expands the interval on both sides and eventually 0042 %! violates the condition nu>2. 0043 %! 0044 %! @end deftypefn 0045 %@eod: 0046 0047 % Copyright (C) 2003-2012 Dynare Team 0048 % 0049 % This file is part of Dynare. 0050 % 0051 % Dynare is free software: you can redistribute it and/or modify 0052 % it under the terms of the GNU General Public License as published by 0053 % the Free Software Foundation, either version 3 of the License, or 0054 % (at your option) any later version. 0055 % 0056 % Dynare is distributed in the hope that it will be useful, 0057 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0058 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0059 % GNU General Public License for more details. 0060 % 0061 % You should have received a copy of the GNU General Public License 0062 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0063 0064 check_solution_flag = 1; 0065 s = []; 0066 nu = []; 0067 0068 if nargin==3 0069 use_fzero_flag = 0; 0070 end 0071 0072 sigma2 = sigma^2; 0073 mu2 = mu^2; 0074 0075 if type == 2; % Inverse Gamma 2 0076 nu = 2*(2+mu2/sigma2); 0077 s = 2*mu*(1+mu2/sigma2); 0078 elseif type == 1; % Inverse Gamma 1 0079 if sigma2 < Inf 0080 nu = sqrt(2*(2+mu2/sigma2)); 0081 if use_fzero_flag 0082 nu = fzero(@(nu)ig1fun(nu,mu2,sigma2),nu); 0083 else 0084 nu2 = 2*nu; 0085 nu1 = 2; 0086 err = ig1fun(nu,mu2,sigma2); 0087 err2 = ig1fun(nu2,mu2,sigma2); 0088 if err2 > 0 % Too short interval. 0089 while nu2 < 1e12 % Shift the interval containing the root. 0090 nu1 = nu2; 0091 nu2 = nu2*2; 0092 err2 = ig1fun(nu2,mu2,sigma2); 0093 if err2<0 0094 break 0095 end 0096 end 0097 if err2>0 0098 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...'); 0099 end 0100 end 0101 % Solve for nu using the secant method. 0102 while abs(nu2/nu1-1) > 1e-14 0103 if err > 0 0104 nu1 = nu; 0105 if nu < nu2 0106 nu = nu2; 0107 else 0108 nu = 2*nu; 0109 nu2 = nu; 0110 end 0111 else 0112 nu2 = nu; 0113 end 0114 nu = (nu1+nu2)/2; 0115 err = ig1fun(nu,mu2,sigma2); 0116 end 0117 end 0118 s = (sigma2+mu2)*(nu-2); 0119 if check_solution_flag 0120 if abs(log(mu)-log(sqrt(s/2))-gammaln((nu-1)/2)+gammaln(nu/2))>1e-7 0121 error('inverse_gamma_specification:: Failed in solving for the hyperparameters!'); 0122 end 0123 if abs(sigma-sqrt(s/(nu-2)-mu^2))>1e-7 0124 error('inverse_gamma_specification:: Failed in solving for the hyperparameters!'); 0125 end 0126 end 0127 else 0128 nu = 2; 0129 s = 2*mu2/pi; 0130 end 0131 else 0132 error('inverse_gamma_specification: unkown type') 0133 end 0134 0135 %@test:1 0136 %$ 0137 %$ [s0,nu0] = inverse_gamma_specification(.75,.2,1,0); 0138 %$ [s1,nu1] = inverse_gamma_specification(.75,.2,1,1); 0139 %$ [s3,nu3] = inverse_gamma_specification(.75,.1,1,0); 0140 %$ [s4,nu4] = inverse_gamma_specification(.75,.1,1,1); 0141 %$ % Check the results. 0142 %$ t(1) = dyn_assert(s0,s1,1e-6); 0143 %$ t(2) = dyn_assert(nu0,nu1,1e-6); 0144 %$ t(3) = isnan(s4); 0145 %$ t(4) = isnan(nu4); 0146 %$ t(5) = dyn_assert(s3,16.240907971002265,1e-6);; 0147 %$ t(6) = dyn_assert(nu3,30.368398202624046,1e-6);; 0148 %$ T = all(t); 0149 %@eof:1