how to generate random numbers from a distribution which is the convolution of exponential distributions?

Hello,

I would like to simulate two independent laser speckle intensities. Based on J. Goodman's book (Speckle Phenomena in Optics, p.39), the formula for the intensity distribution is the convolution of two exponential distributions:

p(i) = i/(I^2) exp(-i/I)

with I being the average intensity.

Is there already a function doing this? If not, how should I proceed to write it myself?

Thanks for the answer!

Cheers Guillaume

Answer 1

[Shoaibur Rahman]

i = 0:0.1:10; I1 = 1; I2 = 1; % define as per your requirment
p1 = i/(I1^2).*exp(-i/I1); % first exponential function
p2 = i/(I2^2).*exp(-i/I2); % second exponential function

y = conv(p1,p2); % convolution of two functions

N = length(y);
r = randi(N,[1 N]); % generate random indices, uniform in this case

yourRandomNumber = y(r) % access y data at random indices

Answer 2

[Roger Stafford]

By integrating from 0 to infinity it can readily be determined that the cdf for your density p(x) = x/I^2*exp(-x/I) would be

   u = F(x) = 1-(1+x/I)*exp(-x/I)

To generate a random variable with this cdf we need the inverse of this function. There may be one already lurking somewhere in the Statistics Toolbox, but I am not aware of it. However, the inverse can be computed using the 'lambertw' function. Its inverse is given by the matlab function 'lambertw' as:

   x = (-lambertw(-1,(u-1)/exp(1))-1)*I;

where the initial -1 argument in the lambertw function indicates the appropriate branch for this many-branched function. Consequently you can generate your desired random variable with the single line:

   x = (-real(lambertw(-1,(rand(n,1)-1)/exp(1)))-1)*I;

where n is the desired number of samples with this distribution. (Note: I have taken the real part of the lambertw output here because, at least for my copy of this function, it can produce very tiny imaginary parts due presumably to round-off error.)