rMVNorm {bayesSurv}

Sample from the multivariate normal distribution

Description

According to the parametrization used, sample from the multivariate normal distribution.

The following parametrization can be specified

standard

In this case we sample from either \mathcal{N}(\mu, \Sigma) or from \mathcal{N}(\mu, Q^{-1}).

canonical

In this case we sample from \mathcal{N}(Q^{-1}b,\;Q^{-1}).

Generation of random numbers is performed by Algorithms 2.3-2.5 in Rue and Held (2005, pp. 34-35).

Usage

  rMVNorm(n, mean=0, Sigma=1, Q, param=c("standard", "canonical"))

Arguments

n

number of observations to be sampled.

mean

For param="standard", it is a vector \mu of means. If length(mean) is equal to 1, it is recycled and all components have the same mean.

For param="canonical", it is a vector b of canonical means. If length(mean) is equal to 1, it is recycled and all components have the same mean.

Sigma

covariance matrix of the multivariate normal distribution. It is ignored if Q is given at the same time.

Q

precision matrix of the multivariate normal distribution.

It does not have to be supplied provided Sigma is given and param="standard".

It must be supplied if param="canonical".

param

a character which specifies the parametrization.

Value

Matrix with sampled points in rows.

Author(s)

Arnošt Komárek [e-mail]

References

Rue, H. and Held, L. (2005). Gaussian Markov Random Fields: Theory and Applications. Boca Raton: Chapman and Hall/CRC.

See Also

rnorm, Mvnorm.

Examples

### Mean, covariance matrix, its inverse
### and the canonical mean
mu <- c(0, 2, 0.5)
L <- matrix(c(1, 1, 1,  0, 0.5, 0.5,  0, 0, 0.3), ncol=3)
Sigma <- L %*% t(L)
Q <- chol2inv(t(L))
b <- Q %*% mu

print(Sigma)
print(Q)
print(Sigma %*% Q)

### Sample using different parametrizations
set.seed(775988621)
n <- 10000

### Sample from N(mu, Sigma)
xx1 <- rMVNorm(n=n, mean=mu, Sigma=Sigma)
apply(xx1, 2, mean)
var(xx1)

### Sample from N(mu, Q^{-1})
xx2 <- rMVNorm(n=n, mean=mu, Q=Q)
apply(xx2, 2, mean)
var(xx2)

### Sample from N(Q^{-1}*b, Q^{-1})
xx3 <- rMVNorm(n=n, mean=b, Q=Q, param="canonical")
apply(xx3, 2, mean)
var(xx3)