Wishart distribution
Wishart distribution
Wishart(nu, S),
where nu are degrees of freedom of the Wishart distribution and S is its scale matrix. The same parametrization as in Gelman (2004) is assumed, that is, if W~Wishart(nu,S) then
E(W) = nu*S.
Prior to version 3.4-1 of this package, functions dWISHART and
rWISHART were called as dWishart and rWishart,
respectively. The names were changed in order to avoid conflicts with
rWishart from a standard package stats.
dWISHART(W, df, S, log=FALSE) rWISHART(n, df, S)
W |
Either a matrix with the same number of rows and columns as
|
n |
number of observations to be sampled. |
df |
degrees of freedom of the Wishart distribution. |
S |
scale matrix of the Wishart distribution. |
log |
logical; if |
The density of the Wishart distribution is the following
f(W) = (2^{nu*p/2} * pi^{p*(p-1)/4} * prod[i=1]^p Gamma((nu + 1 - i)/2))^{-1} * |S|^{-nu/2} * |W|^{(nu - p - 1)/2} * exp(-0.5*tr(S^{-1}*W)),
where p is number of rows and columns of the matrix W.
In the univariate case, Wishart(nu,S) is the same as Gamma(nu/2, 1/(2*S)).
Generation of random numbers is performed by the algorithm described in Ripley (1987, pp. 99).
Some objects.
A numeric vector with evaluated (log-)density.
If n equals 1 then a sampled symmetric matrix W is returned.
If n > 1 then a matrix with sampled points (lower triangles of
W) in rows is returned.
Arnošt Komárek arnost.komarek[AT]mff.cuni.cz
Gelman, A., Carlin, J. B., Stern, H. S., and Rubin, D. B. (2004). Bayesian Data Analysis, Second edition. Boca Raton: Chapman and Hall/CRC.
Ripley, B. D. (1987). Stochastic Simulation. New York: John Wiley and Sons.
set.seed(1977) ### The same as gamma(shape=df/2, rate=1/(2*S)) df <- 1 S <- 3 w <- rWISHART(n=1000, df=df, S=S) mean(w) ## should be close to df*S var(w) ## should be close to 2*df*S^2 dWISHART(w[1], df=df, S=S) dWISHART(w[1], df=df, S=S, log=TRUE) dens.w <- dWISHART(w, df=df, S=S) dens.wG <- dgamma(w, shape=df/2, rate=1/(2*S)) rbind(dens.w[1:10], dens.wG[1:10]) ldens.w <- dWISHART(w, df=df, S=S, log=TRUE) ldens.wG <- dgamma(w, shape=df/2, rate=1/(2*S), log=TRUE) rbind(ldens.w[1:10], ldens.wG[1:10]) ### Bivariate Wishart df <- 2 S <- matrix(c(1,3,3,13), nrow=2) print(w2a <- rWISHART(n=1, df=df, S=S)) dWISHART(w2a, df=df, S=S) w2 <- rWISHART(n=1000, df=df, S=S) print(w2[1:10,]) apply(w2, 2, mean) ## should be close to df*S (df*S)[lower.tri(S, diag=TRUE)] dens.w2 <- dWISHART(w2, df=df, S=S) ldens.w2 <- dWISHART(w2, df=df, S=S, log=TRUE) cbind(w2[1:10,], data.frame(Density=dens.w2[1:10], Log.Density=ldens.w2[1:10])) ### Trivariate Wishart df <- 3.5 S <- matrix(c(1,2,3,2,20,26,3,26,70), nrow=3) print(w3a <- rWISHART(n=1, df=df, S=S)) dWISHART(w3a, df=df, S=S) w3 <- rWISHART(n=1000, df=df, S=S) print(w3[1:10,]) apply(w3, 2, mean) ## should be close to df*S (df*S)[lower.tri(S, diag=TRUE)] dens.w3 <- dWISHART(w3, df=df, S=S) ldens.w3 <- dWISHART(w3, df=df, S=S, log=TRUE) cbind(w3[1:10,], data.frame(Density=dens.w3[1:10], Log.Density=ldens.w3[1:10]))
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.