Simulating Gamma distributed RV's using Matlab

I am a beginner in Matlab and I need an explanation for this. I am trying to generate random variables that are gamma distributed and compare them with the output of gamcdf. I edited a code that does the same purpose with the exponential distribution but it seems that i made some mistakes.

n = 10000;
% inline function: gamma distributed random variables %generated by the sum of exponentials
randGamma = inline('-sum(log(rand(1,n)))/a', 'n', 'a');
a = 2
x = randGamma(n,a);
t = linspace(0, 5/a, 500);
figure
plot(sort(x), (1:n)/n, 'rx'), hold on
plot(t, gamcdf(t,a), 'k')
axis([0, max(t), 0, 1.1])
title(['Gammadistribution a = ', num2str(a)])
end

in the figures I get x is not showing up. I actually have a doubt in

plot(sort(x), (1:n)/n, 'rx'), hold on.

Can anyone please help me with that. I will have to do such an implementation with other distributions. So understanding this one will help me a lot. Thank you in advance

Answer 1

[Paul]

For starters, probably shouldn't use inline. Use an anonymous function instead. And I'm going to change the variables involved to be consistent with Matlab's definitions.

randGamma = @(n,lamda) (-sum(log(rand(1,n)))./lamda);

But the bigger problem is that randGamma, as defined, only generates a single output for inputs n and a. Won't you have to call it many times to generate the samples you seek?

Did you really mean that the gamma-distributed random variable is supposed to be the sum of 10000 exponentially-distributed random variables?

Let's assume that random variable G is the sum of 10 i.i.d random variables with exponential distribution with mean 2

n = 10;
mu = 2;
lamda = 1/mu;

Now generate an array of G using randGamma and using exprnd

ntrials = 1000;
G1 = nan(ntrials,1);
G2 = G1;
for ii = 1:ntrials
G1(ii) = sum(exprnd(mu,1,n));
G2(ii) = randGamma(n,lamda);
end

Now compare the experimental and exact CDFs

cdfplot(G1);
hold on
cdfplot(G2)
a = n; b = 1/lamda;
plot(0:50,gamcdf(0:50,a,b))
legend('exprnd','randGamma','gamcdf')

[Tarek Mahjoub]

First thank you Paul for the great help. I actually used the

(-sum(log(rand(1,n)))./a)

because I used the inverse transform method to generate exponential random variables (where X = -1/a log(U) where U is uniform distributed random variable) . So it would be better to use this expression rather than exprnd.

[Paul]

Sure, I was just showing that your approach matches what would result from using exprnd. But I didn't want to use 'a' as the argument to randGamma becuse Matlab uses 'a' as the first parameter of the gamma distribution. Alos you could avoid the loop with

randGamma = @(n,lamda,ntrials) (-sum(log(rand(ntrials,n)),2)./lamda);

[Tarek Mahjoub]

Thanks a lot for your help