trhoangdung / trun_mvnt Goto Github PK

Sampling from Truncated Multivariate Normal and t Distributions under linear inequality constraints

This module implements the Efficient sampling algorithm of truncated multivariate (scale) mixtures of normals under linear inequality constraints proposed by Li and Ghosh (2015) (doi:10.1080/15598608.2014.996690) under Python environment (analogous to the R package tmvmixnorm).

The two cores function rtmvn() and rtmvt() are useful to overcoming difficulties in simulating truncated normal and Student's t distribution respectively. Efficient rejection sampling for simulating truncated univariate normal distribution is also included in the modeule, which shows superiority in terms of acceptance rate and numerical stability compared to existing methods. An efficient function for sampling from truncated multivariate normal distribution subject to convex polytope restriction regions based on Gibbs sampler.

author: Ting Fung (Ralph) Ma

email: [email protected] (feel free to email me if you find any bug!)

Summary:

In short, this module can be used to generate random sample of truncated multivariate normal (and Student't). Note that the truncation can be in the form of many (linear) constraints.

Suppose we want to draw sample from p-dimensioanl normal ), and similarly for Student's t with df=,

subject to constranst(s): , where is , even if .

The core functions rtmvn() and rtmvt() can solve the problem well.

Examples:

import numpy as np
from trun_mvnt import rtmvn, rtmvt


########## Traditional problem  ##########
##### lower < X < upper #####
# So D = identity matrix

D = np.diag(np.ones(4))
lower = np.array([-1,-2,-3,-4])
upper = -lower
Mean = np.zeros(4)
Sigma = np.diag([1,2,3,4])

n = 10 # want 500 final sample
burn = 100 # burn-in first 100 iterates
thin = 1 # thinning for Gibbs


random_sample = rtmvn(n, Mean, Sigma, D, lower, upper, burn, thin) 
# Numpy array n-by-p as result!
random_sample

# Suppose you want a Student t

nu = 10
random_sample_t = rtmvt(n, Mean, Sigma, nu, D, lower, upper, burn, thin)
# Numpy array n-by-p as result!
random_sample_t

########## Non-full rank problem (more constraints than dimension) ##########
Mean = np.array([0,0])
Sigma = np.array([1, 0.5, 0.5, 1]).reshape((2,2)) # bivariate normal

D = np.array([1,0,0,1,1,-1]).reshape((3,2)) # non-full rank problem
lower = np.array([-2,-1,-2])
upper = np.array([2,3,5])

n = 500 # want 500 final sample
burn = 100 # burn-in first 100 iterates
thin = 1 # thinning for Gibbs

random_sample = rtmvn(n, Mean, Sigma, D, lower, upper, burn, thin) # Numpy array n-by-p as result!