Code covered by the BSD License  

Highlights from
Truncated multivariate normal

4.5 | 2 ratings Rate this file 53 Downloads (last 30 days) File Size: 7.2 KB File ID: #34402 Version: 1.5
image thumbnail

Truncated multivariate normal


Tim Benham (view profile)


01 Jan 2012 (Updated )

Generates pseudo-random vectors drawn from the truncated multivariate normal distribution.

| Watch this File

File Information

   X = rmvnrnd(MU,SIG,N,A,B) returns in N-by-P matrix X a
   random sample drawn from the P-dimensional multivariate normal
   distribution with mean MU and covariance SIG truncated to a
   region bounded by the hyperplanes defined by the inequalities Ax<=B.
   [X,RHO,NAR,NGIBBS] = rmvnrnd(MU,SIG,N,A,B) returns the
   acceptance rate RHO of the accept-reject portion of the algorithm
   (see below), the number NAR of returned samples generated by
   the accept-reject algorithm, and the number NGIBBS returned by
   the Gibbs sampler portion of the algorithm.
   rmvnrnd(MU,SIG,N,A,B,RHOTHR) sets the minimum acceptable
   acceptance rate for the accept-reject portion of the algorithm
   to RHOTHR. The default is the empirically identified value


Truncated Gaussian and Chebycenter(A,B,R0) inspired this file.

MATLAB release MATLAB 7.9 (R2009b)
MATLAB Search Path
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (10)
16 Aug 2015 Tim Benham

Tim Benham (view profile)

Giorgos, I have updated the function to correct a bug in the initialization of the Gibbs sampler. Please try it.

Comment only
16 Aug 2015 Tim Benham

Tim Benham (view profile)

Thanks for trying my function Giorgos. Would it be possible for you to send me an example of the failure case?

Comment only
15 Aug 2015 Giorgos Minas


I used the function and generally worked well but in somewhat 'hard conditions' I experienced the issue explained below. First the 'hard conditions' are that I used the function to derive trajectories of a multivariate stochastic process in which all variables should be positive, but the dynamics driving the process can get the trajectory to a position where the next step has negative mean in most variables and all variables are highly correlated.

In those conditions, rmvnrnd returned inf in one variable and nan to the others, the next steps then obviously also all nans. Having no time to fix the problem, I just dropped those trajectories and start again, but as in one run of 3000 trajectories the issue occurred in 600 situations, this is a great cost, as you can imagine.

It would be great to see an attempt of a good fix to this problem as otherwise the function works fine.


Comment only
14 Dec 2014 Tim Benham

Tim Benham (view profile)

True, I fixed that in the R version. I will incorporate your fix into the MATLAB version.

Comment only
14 Dec 2014 Kris Villez


I just tested and used this package. I noticed that the Gibbs sampler does not necessarily start in a feasible point (within the polygon). For large-dimensional problems, this can make the chance of arriving in the feasible space by random sampling extremely low.

I made a small modification to the code by using the constrained maximum likelihood solution as the first sample. This helps to avoid this situation.

Otherwise, I think this is a great piece of code. Thanks for sharing.


02 Aug 2014 Sergio

Sergio (view profile)

06 Feb 2013 jenka

jenka (view profile)

Also, I think it would significantly improve this code if we had an option to change the size N for each variable p. Let's say for x_1 I would want N=100, for x_2 I would want N=55.

Comment only
05 Feb 2013 jenka

jenka (view profile)

Hi Tim,
yes, thank you for your reply. However, in your GIBBS sampling, it appears that the number of iterations only depends on the sample size N. Let's say the user wants to generate 100 sample from multivariate truncated normal. Then the number of iterations based on your code for Gibbs sampling is always going to be 100.
Also, I read the paper you reference in your code. However, I am not sure that accept-rejection method is from that paper. Perhaps, you could add another reference? Also, could you please add more comments to that part of the code? Thanks!

Comment only
03 Feb 2013 Tim Benham

Tim Benham (view profile)

m = [2 0];
S = [1 0.9; 0.9 1];
X = rmvnrnd(m,S,100,[-eye(2);eye(2)], [0; 0; 1; 1]);
clf; plot(X(:,1),X(:,2),'.')

Comment only
31 Jan 2013 jenka

jenka (view profile)

Could you please provide an example of using your code? I have a vector MU and matrix SIG. I want my realizations to be bounded between [0, 1]. Thanks.

Comment only
03 Feb 2013 1.2

Changed interface to support omission of A and b. Added example script rmvnrnd_eg.

04 Feb 2013 1.3

Added example and improved interface.

16 Aug 2015 1.4

1. Correct problem with initialization of the Gibbs sampler.
2. Behave correctly when no constraints supplied.

26 Aug 2015 1.5

Uses chebycenter to find a feasible point to seed the Gibbs sampler if the the random generation does not find one.

Contact us