tmvtnorm's Introduction

tmvtnorm

tmvtnorm: A package for the Truncated Multivariate Normal Distribution

This package contains a number of useful methods for the truncated multivariate normal distribution. It considers random number generation with rejection and Gibbs sampling, computation of marginal densities as well as computation of the mean and covariance of the truncated variables.

For a more detailed introduction, see this RJournal (2010) paper tmvtnorm: A Package for the Truncated Multivariate Normal Distribution.

tmvtnorm's People

Contributors

Stargazers

Watchers

tmvtnorm's Issues

Error when using `rtmvnorm.gibbs` in one dimension

When working in one dimension, the Gibbs sampler seems to produce random samples with the wrong variance:

library(tmvtnorm)

set.seed(2023)

y_mean <- 0
y_var <- 0.1
y_trunc <- 0
y <- rnorm(1000, sd = sqrt(y_var))

# Moments of truncated distribution
mtmvnorm(mean = y_mean, sigma = y_var, lower = y_trunc)

$tmean
[1] 0.2523133

$tvar
           [,1]
[1,] 0.03633802

# Samples from truncated distribution
y_new <- rtmvnorm(sum(y > y_trunc), mean = y_mean, sigma = y_var,
                  lower = y_trunc, algorithm = "gibbs")
c(mean(y_new), var(y_new))

[1] 0.079204130 0.003461616

And we can see from the histogram that the new samples don't follow the desired distribution:

y[y > y_trunc] <- y_new
hist(y, breaks = 50)

Rejection sampling works as expected:

y_new_rej <- rtmvnorm(sum(y > y_trunc), mean = y_mean, sigma = y_var,
                      lower = y_trunc, algorithm = "rejection")
c(mean(y_new_rej), var(y_new_rej))

[1] 0.24295202 0.03373638

As does Gibbs sampling in 2 independent dimensions:

y_new_2d <- rtmvnorm(sum(y > y_trunc), mean = rep(y_mean, 2),
                     sigma = diag(y_var, 2), lower = rep(y_trunc, 2),
                     algorithm = "gibbs")
colMeans(y_new_2d)

[1] 0.2564023 0.2474838

var(y_new_2d)

            [,1]        [,2]
[1,]  0.03777465 -0.00102825
[2,] -0.00102825  0.04048425

I believe the issue is because rtnorm.gibbs treats the argument sigma as the standard deviation (as per its documentation), but whenever it is called (lines 318, 425, 521 & 523 of rtmvnorm.R), the variance sigma[1,1] is supplied. I think changing either the function or the calls to take the square root should resolve the issue.

(BTW - I am not sure if line 521 is supposed to have 1 / H[1,1] instead of 1 / sigma[1,1] as well)

Recommend Projects