Code covered by the BSD License

### Highlights from Truncated multivariate normal

4.66667
4.7 | 3 ratings Rate this file 13 Downloads (last 30 days) File Size: 7.2 KB File ID: #34402 Version: 1.5

# Truncated multivariate normal

### Tim Benham (view profile)

01 Jan 2012 (Updated )

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

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 and Chebycenter(A,B,R0) inspired this file.

MATLAB release MATLAB 7.9 (R2009b)
MATLAB Search Path
`/`
29 Jun 2016 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?

Comment only
28 Apr 2016 Javier Garcia

### Javier Garcia (view profile)

25 Feb 2016 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]);

Comment only
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

### Giorgos Minas (view profile)

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

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

### Kris Villez (view profile)

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

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

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.