rank constraints in multivariate normals

sample_constraints.R

## Define MVTN parameters
mu = c(-1,0,1)
Sigma = matrix(
  c(1,  0,  0,
    0,  1, .5,
    0, .5,  1 
  ), 3,3)
p = length(mu)

library(mvtnorm)

## Define some test data
N = 10000
X = mvtnorm::rmvnorm(N, mu, Sigma)


# Build rank constraint matrix:
# pairwise adjacent differences are positive,
# with an additional redundant row to make it 
# square
# x is a vector with the ordering we want
Dmat = function(x){
  p = length(x)
  D = matrix(0, p, p)
  for(i in 1:(p-1)){
    D[i, order(x)[i+1]] = 1
    D[i, order(x)[i]] = -1
  }
  ## Last row is redundant but we don't need a full rank matrix
  D[p, which.max(x)] = 1
  D[p, which.min(x)] = -1
  D
}

r = c(2,1,3)
D = Dmat(r)

## Sample using the rank constraint
zz = tmvmixnorm::rtmvn(10000, mu, Sigma,
      D = D,
      lower = matrix(0, p), 
      ## It appears the upper bound cannot be infinity, so choose "wisely"
      upper = matrix(10, p),
      int = NULL, burn = 10, thin = 1)



#### Testing below

## Check the conditional distributions

pairs(zz)

# Work out which observations are consistent with the contraint
ind = apply(apply(X, 1, rank) == r, 2, all)
  
## Plot against one another
qqplot(zz[,1], X[ind, 1])
abline(0,1)

qqplot(zz[,2], X[ind, 2])
abline(0,1)

qqplot(zz[,3], X[ind, 3])
abline(0,1)

## Can I recover the marginals using the orderings from each sample?
X2 = X * NA

for(i in 1:N){
  D = Dmat(X[i,])
  X2[i, ] = tmvmixnorm::rtmvn(1, mu, Sigma,
                              D = D,
                              lower = matrix(0, p), 
                              ## It appears the upper bound cannot be infinity, so choose "wisely"
                              upper = matrix(10, p),
                              int = NULL, burn = 10, thin = 1)
}

colMeans(X2)
cov(X2)

qqplot(X2[,1], X[, 1])
abline(0,1)

qqplot(X2[,2], X[, 2])
abline(0,1)

qqplot(X2[,3], X[, 3])
abline(0,1)