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}.$$

Usage

dWishart(W, df, S, log=FALSE)
rWishart(n, df, S)

Arguments

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

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`

Some objects.

`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.

`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.

`Examples`

Run this codeset.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]))
Run the code above in your browser using DataLab