Code covered by the BSD License  

Highlights from
Truncated multivariate normal

5.0

5.0 | 1 rating Rate this file 34 Downloads (last 30 days) File Size: 6.17 KB File ID: #34402
image thumbnail

Truncated multivariate normal

by

 

01 Jan 2012 (Updated )

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

| Watch this File

File Information
Description

   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
   2.9e-4.

Acknowledgements

Truncated Gaussian inspired this file.

MATLAB release MATLAB 7.9 (R2009b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (5)
02 Aug 2014 Sergio  
06 Feb 2013 jenka

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.

05 Feb 2013 jenka

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!

03 Feb 2013 Tim Benham

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),'.')

31 Jan 2013 jenka

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.

Updates
03 Feb 2013

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

04 Feb 2013

Added example and improved interface.

Contact us