Simulating Gamma distributed RV's using Matlab

3 views (last 30 days)
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

Accepted Answer

Paul
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')
  3 Comments
Paul
Paul on 26 Jul 2021
Edited: Paul on 26 Jul 2021
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);

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!