restrictedMVN-package: Sampler from multivariate normal with affine constraints

In restrictedMVN: Multivariate Normal Restricted by Affine Constraints

Description

The package implements a fast gibbs sampler for the multivariate nomral with affine constraints. For the d-dimensional Z~Normal(mu,sigma), the linear_part matrix A in d x r, and offset vector b in 1 x r define a multivariate normal with affine constraints in {Z| A*Z<= b}.

Details

Sampling is implemented in the main function, sample_from_constraints. It is parameterized by the parameters of the normal (mean_param and covariance), parameters of the restriction (linear_part and offset), and the number of samples ndraw. The user also needs to specify an initial point that satisfies the constraints. thresh2constraints is a helper function that translates coordinate-wise truncations into the affine form.

Author(s)

Jonathan Taylor and Yuval Benjamini.

Maintainer: Yuval Benjamini <[e-mail]>

See Also

The package was originally part of the github selective-inference code base.

Examples

constr = thresh2constraints(3, lower = c(0.2,0.2,0.2))
covariance = matrix(c(1,0.5,0,0.5,1,0.5,0,0.5,1),nc=3)
 
samp = sample_from_constraints(linear_part = constr$linear_part,
                                    offset= constr$offset,
                                    mean_param = c(0,0,0),
                                    covariance = covariance,
                                    initial_point = c(1,1,1), 
                                    ndraw=20000,
                                    burnin=2000)


# all points should be >= 0.2 
stopifnot(all(samp>=0.2))

mean_restricted = colMeans(samp)

# compare to rejection of multivariate normals
library("MASS")
full_samp = mvrnorm(n=100000,mu = c(0,0,0),Sigma = covariance)
# Add restrictions:
pass_restrictions = apply(full_samp, 1, 
          function(x){all(constr$linear_part%*% x - constr$offset <=0 )})
          
cond_samp = full_samp[pass_restrictions,]
mean_restricted_rej = colMeans(cond_samp)

stopifnot(all(abs(mean_restricted - mean_restricted_rej)<=0.05))

Example output

restrictedMVN documentation built on May 1, 2019, 7:39 p.m.