File Exchange

image thumbnail

Truncated multivariate normal

version 1.5 (7.2 KB) by

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

12 Downloads

Updated

View License

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.

Comments and Ratings (14)

Uri Cohen

Manny

Manny (view profile)

When I set a non negativity constraint, with a negative mean and a very small variance, I get the following:

Error using *
Inner matrix dimensions must agree.

Error in rmvnrnd (line 106)
elseif all(A'*mu' <= b')

Is it because it is reaching an infeasible area?

Javier Garcia

Pankaj

Pankaj (view profile)

Hi Tim
I run the file rmvnrnd_eg.m and received the following code, could you have a look at it. Thanks
%
Error using *
Inner matrix dimensions must agree.

Error in rmvnrnd (line 106)
elseif all(A'*mu' <= b')

Error in rmvnrnd_eg (line 11)
X3 = rmvnrnd(m, S, 100, [-1 0; 0 1], [-2; -2]);

Tim Benham

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

Tim Benham

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

Giorgos Minas

Hi,

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.

Best

Tim Benham

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

Kris Villez

Hi,

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.

Best,
Kris

Sergio

Sergio (view profile)

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.

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!

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

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.

Updates

1.5

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

1.4

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

1.3

Added example and improved interface.

1.2

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

MATLAB Release
MATLAB 7.9 (R2009b)
Acknowledgements

Inspired by: chebycenter(A,b,r0), Truncated Gaussian

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video

Win prizes and improve your MATLAB skills

Play today