# Simulating Gamma distributed RV's using Matlab

10 views (last 30 days)
Tarek Mahjoub on 26 Jul 2021
Commented: Tarek Mahjoub on 26 Jul 2021
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])
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

Paul on 26 Jul 2021
Edited: Paul on 26 Jul 2021
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 on 26 Jul 2021
Thanks a lot for your help