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