TMVN {mixAK} | R Documentation |
Random generation for the truncated multivariate normal distribution.
The mean and covariance matrix of the original multivariate normal distribution
are mean
and Sigma
. Truncation limits are given by
a
, b
, type of truncation is given by trunc
.
This function uses a Gibbs algorithm to produce a Markov chain whose stationary distribution is the targeted truncated multivariate normal distribution, see Geweke (1991) for more details. Be aware that the sampled values are not i.i.d.!
rTMVN(n, mean=c(0, 0), Sigma=diag(2), a, b, trunc, xinit)
mean |
a numeric vector of the mean of the original multivariate normal distribution. |
Sigma |
covariance matrix of the original multivariate normal distribution. |
a |
a numeric vector of the same length as mean of
truncation limits 1. |
b |
a numeric vector of the same length as mean of
truncation limits 2. |
trunc |
a numeric vector of the same length as mean
describing the type of truncation in each margin.
If trunc is not given, it is assumed that it is equal to
4. Note that a , b and trunc
must have the same length, with exception that
b does not have to be supplied if all trunc values 0,
1, 2 or 4.
|
xinit |
a numeric vector of the same length as mean with
the initial value for the Gibbs sampler. If it is not supplied, the
function determines itself the initial value. |
n |
number of observations to be sampled. |
A matrix with the sampled values (Markov chain) in rows.
Arnošt Komárek arnost.komarek[AT]mff.cuni.cz
Geweke, J. (1991). Efficient simulation from the multivariate normal and Student-t distributions subject to linear constraints and the evaluation of constraint probabilities. Computer Sciences and Statistics, 23, 571–578.
set.seed(1977) exam2 <- function(n, mu, sigma, rho, a, b, trunc) { Sigma <- matrix(c(sigma[1]^2, rho*sigma[1]*sigma[2], rho*sigma[1]*sigma[2], sigma[2]^2), nrow=2) x <- rTMVN(n, mean=mu, Sigma=Sigma, a=a, b=b, trunc=trunc) x1.gr <- seq(mu[1]-3.5*sigma[1], mu[1]+3.5*sigma[1], length=100) x2.gr <- seq(mu[2]-3.5*sigma[2], mu[2]+3.5*sigma[2], length=100) z <- cbind(rep(x1.gr, 100), rep(x2.gr, each=100)) dens.z <- matrix(dMVN(z, mean=mu, Sigma=Sigma), ncol=100) MEAN <- round(apply(x, 2, mean), 3) SIGMA <- var(x) SD <- sqrt(diag(SIGMA)) RHO <- round(SIGMA[1,2]/(SD[1]*SD[2]), 3) SD <- round(SD, 3) layout(matrix(c(0,1,1,0, 2,2,3,3), nrow=2, byrow=TRUE)) contour(x1.gr, x2.gr, dens.z, col="darkblue", xlab="x[1]", ylab="x[2]") points(x[,1], x[,2], col="red") title(sub=paste("Sample mean = ", MEAN[1], ", ", MEAN[2], ", Sample SD = ", SD[1], ", ", SD[2], ", Sample rho = ", RHO, sep="")) plot(1:n, x[,1], type="l", xlab="Iteration", ylab="x[1]", col="darkgreen") plot(1:n, x[,2], type="l", xlab="Iteration", ylab="x[2]", col="darkgreen") return(x) } x1 <- exam2(1000, mu=c(-1, 1), sigma=c(1, sqrt(2)), rho=0, a=c(-6, -9), b=c(4, 11), trunc=c(3, 3)) x2 <- exam2(1000, mu=c(-1, 1), sigma=c(1, sqrt(2)), rho=0.7, a=c(-6, -9), b=c(4, 11), trunc=c(3, 3)) x3 <- exam2(1000, mu=c(-1, 1), sigma=c(1, sqrt(2)), rho=0.7, a=c(-100, -100), b=c(100, 100), trunc=c(3, 3)) x4 <- exam2(1000, mu=c(-1, 1), sigma=c(1, sqrt(2)), rho=-0.7, a=c(-6, -9), b=c(4, 11), trunc=c(3, 3)) x5 <- exam2(1000, mu=c(-1, 1), sigma=c(1, sqrt(2)), rho=-0.9, a=c(-6, -9), b=c(4, 11), trunc=c(3, 3)) x6 <- exam2(1000, mu=c(-1, 1), sigma=c(1, sqrt(2)), rho=0.7, a=c(0, 0), trunc=c(0, 2)) x7 <- exam2(1000, mu=c(-1, 1), sigma=c(1, sqrt(2)), rho=0.7, a=c(-1, 1), trunc=c(0, 2)) x8 <- exam2(1000, mu=c(-1, 1), sigma=c(1, sqrt(2)), rho=0.7, a=c(-1, 1), trunc=c(1, 2)) x9 <- exam2(1000, mu=c(-1, 1), sigma=c(1, sqrt(2)), rho=0.7, a=c(-1.5, 0.5), b=c(-0.5, 1.5), trunc=c(3, 3)) x10 <- exam2(1000, mu=c(-1, 1), sigma=c(1, sqrt(2)), rho=0.7, a=c(-1.5, 0.5), b=c(-0.5, 1.5), trunc=c(4, 3))