mixAK (version 5.4)

Wishart: Wishart distribution

Description

Wishart distribution $$\mbox{Wishart}(\nu, \boldsymbol{S}),$$ where \(\nu\) are degrees of freedom of the Wishart distribution and \(\boldsymbol{S}\) is its scale matrix. The same parametrization as in Gelman (2004) is assumed, that is, if \(\boldsymbol{W}\sim\mbox{Wishart}(\nu,\,\boldsymbol{S})\) then $$\mbox{E}(\boldsymbol{W}) = \nu \boldsymbol{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.

Usage

dWISHART(W, df, S, log=FALSE)
rWISHART(n, df, S)

Value

Some objects.

Arguments

W

Either a matrix with the same number of rows and columns as S (1 point sampled from the Wishart distribution) or a matrix with ncol equal to ncol*(ncol+1)/2 and n rows (n points sampled from the Wishart distribution for which only lower triangles are given in rows of the matrix W).

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 TRUE, log-density is computed

Value for dWISHART

A numeric vector with evaluated (log-)density.

Value for rWISHART

If n equals 1 then a sampled symmetric matrix W is returned.

If n > 1 then a matrix with sampled points (lower triangles of \(\boldsymbol{W}\)) in rows is returned.

Author

Arnošt Komárek

Details

The density of the Wishart distribution is the following $$ f(\boldsymbol{W}) = \Bigl(2^{\nu\,p/2}\,\pi^{p(p-1)/4}\,\prod_{i=1}^p \Gamma\bigl(\frac{\nu + 1 - i}{2}\bigr)\Bigr)^{-1}\, |\boldsymbol{S}|^{-\nu/2}\,|\boldsymbol{W}|^{(\nu - p - 1)/2}\, \exp\Bigl(-\frac{1}{2}\mbox{tr}(\boldsymbol{S}^{-1}\boldsymbol{W})\Bigr), $$ where \(p\) is number of rows and columns of the matrix \(\boldsymbol{W}\).

In the univariate case, \(\mbox{Wishart}(\nu,\,S)\) is the same as \(\mbox{Gamma}(\nu/2, 1/(2S)).\)

Generation of random numbers is performed by the algorithm described in Ripley (1987, pp. 99).

References

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.

See Also

rWishart.

Examples

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]))