Generating random Gamma distribution

i would like to generate random gamma distributuon with variate size = 100, alpha=10, beta = 2.

but i dont want to use toolbox like this

gammarandom = gamrnd(10, 2, [100, 1]);

i want to generate randim number by using a gamma formula (not toolbox) .

can you please help me? thanks a lot.

[mohamed ali]

I want to Generating random Generalized gamma Distribution Data in Matlab How?

Answer 1

[John D'Errico]

The trick is to understand how a transformation can solve the problem. This is a standard way to produce random numbers with some given distribution. You need to know the CDF of the distribution, and be able to invert that CDF. For example, suppose I wish to generate normally distributed random variates? Yes, I know you can use normrnd, and randn is even easier. But suppose we wanted to the work ourselves?

The idea is to start with a uniform sample.

N = 1e6;
U = rand(1,N);

So U is uniformly distributed on the interval (0,1).

Now, we map those points backward through the inverse normal CDF. I can do it using norminv, but I can also do it using erfinv, since they are with a linear transformtion of being the same mapping. Thus

R = norminv(U);

mean(R)
ans =
    -0.001839
    
var(R)
ans =
       1.0007   
       
hist(R,100)

In fact, as N increases, we will see that R has mean zero and variance 1. The vector R is normally distributed, as if I had used randn to do so.

I did not really need to use norminv here. Perhaps I can convince you of that, by a simple expedient.

R = sqrt(2)*erfinv(U*2 - 1);

mean(R)
ans =
    -0.001839
    
var(R)
ans =
       1.0007

As I said, norminv is just like erfinv, but wearing a mask.

The idea is, we can sample from some other distribution by using an inverse mapping through the desired CDF. THis works nicely as long as we have that CDF, and it is efficient to perform the inverse mapping.

Next, we can consider how to do the same for a gamma variate. Still easy peasy.

We start the same way, with a uniform sample.

N = 1e6;
U = rand(1,N);

As it turns out, the gamma function does the work for us, in one of its various incarnations - the inverse incomplete gamma.

https://en.wikipedia.org/wiki/Gamma_distribution

On that page, we can learn the gamma PDF and the gamma CDF are very nicely related to the incomplete gamma function itself. So the pdf is

f(x;alpha,beta) = beta^alpha * x^(alpha-1) * exp(-beta*x)/gamma(alpha)

The CDF is as simple, being effectively the incomplete gamma integral, when properly normalized to have unit area. So now, for alpha and beta parameters in the gamma distribution, we can form gamma variates easily by the inverse mapping.

a = 2;
b = 3;
R = b * gammaincinv(U,a);

Did that work? It certainly looks gamma-ish. I know, that is not a word. So sue me.

hist(R,100)

>> mean(R)
ans =
       5.9933
       
>> var(R)
ans =
       17.992
       
>> a*b
ans =
     6
     
>> a*b^2
ans =
    18

The mean and variance certainly seem to be approaching the theoretical parameters.

I hope you will accept this to be the solution to your question.

Just for kicks, suppose we wanted to generate some distribution where we don't have a nice custom made CDF to invert? Perhaps I'll generate a triangular distribution using the above scheme? Yes, we can do that quite easily by taking the sum of two uniform random variates. But this time I want to do it the hard way.

The desired CDF is a piecewise thing now. On the interval [0,1] we will have the pdf as simply x. On the interval [1,2] it drops down again linearly to 0. If my pencil and paper is correct, that makes it look like 2-x on the second half.

syms x
P = piecewise(x <= 0,0,x > 0 & x <= 1,x,x > 1 & x <= 2,2-x,x>2,0);

int(P,[0,2])
ans =
1

fplot(P,[0,2])

It integrates to 1, so it is well defined as a pdf. The CDF is easy enough to create.

C = int(P,[0,x])
C =
piecewise(x == 0, 0, x == 1, 1/2, 2 <= x, 1, x <= 1 & x ~= 0, x^2/2, x in Dom::Interval(1, [2]), 1/2 - ((x - 1)*(x - 3))/2)
fplot(C,[0,2])

Surprise! It works. A nice, well defined CDF. In this case a piecewise quadratic function, thus on the intervals

[0,1] : x^2/2

[1,2] : 1/2 - ((x - 1)*(x - 3))/2)

It reaches 1 at the upper end point, when x==2.

How do we work with this to now sample from a triangular random variable? Start from the very same place!

N = 1e6;
U = rand(1,N);

For u less than 1/2, we invert the lower segment, between 1/2 and 1, we invert the upper segment.

syms u
Cinv = piecewise(u <= 1/2,sqrt(2*u),u>1/2,2 - sqrt(2)*sqrt(1 - u))
fplot(Cinv,[0,1])

Which still looks good as the inverse CDF.

R = zeros(size(U));
ind = U <= 1/2;
R(ind) = sqrt(2*U(ind));
R(~ind) = 2 - sqrt(2)*sqrt(1 - U(~ind));
hist(R,100)

It worked as it should - a triangular random sampling, done using the inverse CDF method.

[Zeynep Toprak]

hmmm thank you for explanation. but how can I link what you have stated with the gamma distribution?

[John D'Errico]

I'm not sure what you are asking. Click on the link. That takes you directly to a Wikipedia page for the gamma distribution. Or paste the link into your web browser. Either should work.