# mcmc

### David Shera (view profile)

MCMC -- Markov Chain Monte Carlo Tools

wishrnd(Sc,nu)
```% WISHRND - Random Matrix from Wishart Distribution
%
%   [W] = wishrnd(Sc,nu)
%
% W = returned random symmetric positive definite matrix
%
% Sc = p x p symmetric, postitive definite "scale" matrix
% nu = "degrees of freedom" (when integer)
%    = "number of observation" (when integer)
%    = precision parameter (when non-integer)
%
% uses the Odell and Feiveson (1966) algorithm
% as printed in Kennedy and Gentle (1980)
%
% Note:
%   Different sources use different parameterizations.
%   See INVWISHRND for details.
%

function [W] = wishrnd(Sc,nu)

[p,p2] = size(Sc) ;
z = normrnd(0,1,p,p) ;

% y = chi2rnd(nu-(1:p))
% --- note, matlab's chi2 functions do not work for non-integer DF.
% --- thus, we use the gamma equivalent

y = gamrnd( (nu-(1:p))/2, 2 ) ;
b = NaN*Sc ;

% first element
b(1,1) = y(1) ;

if (p>1),
% rest of diagonal
for (j=2:p),
zz = z( 1:(j-1) ,j) ;
b(j,j) = y(j) + zz'*zz ;
end

%first row and column
for (j=2:p),
b(1,j) = z(1,j) * sqrt(y(1)) ;
b(j,1) = b(1,j) ;  % mirror
end
end

if p>2,
for (j=3:p),
for ( i=2:(j-1) ),
ix = 1:(i-1) ;
zki = z(ix,i) ;
zkj = z(ix,j) ;
b(i,j) = z(i,j)*sqrt(y(i)) + zki'*zkj ;
b(j,i) = b(i,j) ;   %mirror
end
end
end

[ A, nonposdef ] = chol(Sc) ;

if nonposdef,
disp('ERROR: wishrnd: parameter is not positive-definite')
Argument = A
end

W = A'*b*A ;
```