Sample from multivariate exponential distribution

I'd like to generate random vectors according to a multivariate exponential distribution, that is with a pdf f: R^n->R given by

for appropriate normalisation constant c_e. I wondered if there was built in functionality to do this, if not how would I go about doing this manually?

[the cyclist]

Are the x variables independent, or correlated?

Answer 1

[Roger Stafford]

You can make use of the ‘gammaincinv’ function for this problem:

   X = randn(m,n);
   X = X.*repmat(gammaincinv(rand(m,1),n)./sqrt(sum(X.^2,2)),1,n);

The resulting array X has m rows, each consisting of n coordinates in R^n space with the requested distribution.

(I’m assuming here that ‘gammaincinv’ will accept the scalar n for its second argument. If not, the n will have to be repeated m times using ‘repmat’. I have also assumed that the n coordinate variables are to be statistically independent.)

[Roger Stafford]

I should probably have given an explanation for the above answer.

For each element of X = randn(m,n) the pdf is 1/sqrt(2*pi)*exp(-x^2/2). Since they are regarded as independent, the n-dimensional pdf of each i-th row in X is the product of its elements’ separate pdf’s:

1/sqrt(2*pi)^n*prod(exp(-x(i,:)^2/2)) =
1/sqrt(2*pi)^n*exp(-sum(x(i,:)^2)/2)

and is therefore only a function of r, the distance from the origin in R^n. That means that each row of

X./sqrt(sum(X.^2,2))

is located entirely on the unit n-sphere surface and is uniformly distributed throughout that surface.

The cumulative distribution function cdf of the solid sphere of radius R, as per your request, is to be:

cdf(r<=R) = K*integral from 0 to R of (exp(-r)*S*r^(n-1) dr
          = gammainc(R,n)

where S is the surface measure of a unit n-sphere and where K is the appropriate normalizing constant. To generate this distribution using the rand function in the usual way, we need to use the inverse of gammainc:

gammaincinv(rand,n)

All of this finally gives us the second line in the given answer:

X = X.*repmat(gammaincinv(rand(m,1),n)./sqrt(sum(X.^2,2)),1,n);

[ayoub bouayaben]

@Roger Stafford hello, can you explain me please how can i modify the code that you mentioned before if i want to sample a vector of N elements from a gaussian distribution which is given by :

p(y) = exp(-1/2 *(y - a).' * C* ( y - a) ) with C is the cov matrix with y in R^n

Thank you for your help !

[the cyclist]

Can you use just use the mvnrnd function for that?