rand_inverse_wishart

 function G = rand_inverse_wishart(m, v, H_inv_upper_chol)
 rand_inverse_wishart  Pseudo random matrices drawn from an
 inverse Wishart distribution
 G = rand_inverse_wishart(m, v, H_inv_upper_chol)
 Returns an m-by-m matrix drawn from an inverse-Wishart distribution.

 INPUTS:
     m:          dimension of G and H_inv_upper_chol.
     v:          degrees of freedom, greater or equal than m.
     H_inv_chol: upper cholesky decomposition of the inverse of the
                 matrix parameter.
                 The upper cholesky of the inverse is requested here
                 in order to avoid to recompute it at every random draw.
                 H_inv_upper_chol = chol(inv(H))
 OUTPUTS:
     G:          G ~ IW(m, v, H) where H = inv(H_inv_upper_chol'*H_inv_upper_chol)
                 or, equivalently, using the correspondence between Wishart and
                 inverse-Wishart: inv(G) ~ W(m, v, S) where 
                 S = H_inv_upper_chol'*H_inv_upper_chol = inv(H)
  
 SPECIAL REQUIREMENT
     none

0001 function G = rand_inverse_wishart(m, v, H_inv_upper_chol)
0002 
0003 % function G = rand_inverse_wishart(m, v, H_inv_upper_chol)
0004 % rand_inverse_wishart  Pseudo random matrices drawn from an
0005 % inverse Wishart distribution
0006 % G = rand_inverse_wishart(m, v, H_inv_upper_chol)
0007 % Returns an m-by-m matrix drawn from an inverse-Wishart distribution.
0008 %
0009 % INPUTS:
0010 %     m:          dimension of G and H_inv_upper_chol.
0011 %     v:          degrees of freedom, greater or equal than m.
0012 %     H_inv_chol: upper cholesky decomposition of the inverse of the
0013 %                 matrix parameter.
0014 %                 The upper cholesky of the inverse is requested here
0015 %                 in order to avoid to recompute it at every random draw.
0016 %                 H_inv_upper_chol = chol(inv(H))
0017 % OUTPUTS:
0018 %     G:          G ~ IW(m, v, H) where H = inv(H_inv_upper_chol'*H_inv_upper_chol)
0019 %                 or, equivalently, using the correspondence between Wishart and
0020 %                 inverse-Wishart: inv(G) ~ W(m, v, S) where
0021 %                 S = H_inv_upper_chol'*H_inv_upper_chol = inv(H)
0022 %
0023 % SPECIAL REQUIREMENT
0024 %     none
0025 %
0026 
0027 % Copyright (C) 2003-2009 Dynare Team
0028 %
0029 % This file is part of Dynare.
0030 %
0031 % Dynare is free software: you can redistribute it and/or modify
0032 % it under the terms of the GNU General Public License as published by
0033 % the Free Software Foundation, either version 3 of the License, or
0034 % (at your option) any later version.
0035 %
0036 % Dynare is distributed in the hope that it will be useful,
0037 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0038 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0039 % GNU General Public License for more details.
0040 %
0041 % You should have received a copy of the GNU General Public License
0042 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0043 
0044 X = randn(v, m) * H_inv_upper_chol; 
0045 
0046 
0047 % At this point, X'*X is Wishart distributed
0048 % G = inv(X'*X);
0049 
0050 % Rather compute inv(X'*X) using the SVD
0051 [U,S,V] = svd(X, 0);
0052 SSi = 1 ./ (diag(S) .^ 2);
0053 G = (V .* repmat(SSi', m, 1)) * V';