Probability Distribution Function and Shape

The Gamma Distribution

A random variable X has a Gamma( , ) distribution if and only its probability density is given by

for x > 0

f ( x ) =

0 elsewhere

where > 0.

The following code will draw the density function for the Gamma( ) distribution for , but you are encouraged to try other values. Remember to make sure that they are both positive numbers.

> restart;

> with(plots):

> alpha:=2; beta:=3;

> f(x):=GammaPDF(alpha,beta,x);

> plot(f(x),x=0..20, title="Gamma(2,3) PDF");

To see the effect of on the shape of the Gamma( , ) distribution, the following animation will draw a series of probability density functions as varies from 0.1 to 2.1 by increments of 0.1 while holding constant at 3. Do you see why is a "shape" parameter?

> alpha[0]:=0.1; alpha_step:=0.1; beta:=3;

> for n from 0 to 20 do

> density[n]:=plot(GammaPDF(alpha[0]+n*alpha_step,beta,x),x=0..20):

> num:=convert(evalf(alpha[0]+n*alpha_step,3),string):

> tracker[n]:=textplot([14,0.3,`alpha is `.num],color=blue);

> P[n]:=display({density[n],tracker[n]}):

> od:

> display([seq(P[n], n=0..20)], insequence=true, title="Ranging Alpha");

This next bit of code will animate a Gamma( , ) with a dynamic , holding constant at 2. will vary from 0.2 to 3 at intervals of 0.2. Do you see why is a "scale" parameter?

> alpha:=2; beta[0]:=0.2; beta_step:=0.1;

> for n from 0 to 30 do

> density[n]:=plot(GammaPDF(alpha,beta[0]+beta_step*n,x),x=0..20):

> num:=convert(evalf(beta[0]+n*beta_step,3),string):

> tracker[n]:=textplot([14,1,`beta is `.num],color=blue);

> density[n]:=plot(GammaPDF(alpha,beta[0]+beta_step*n,x),x=0..20):

> P[n]:=display({density[n],tracker[n]}):

> od:

> display([seq(P[n], n=0..30)], insequence=true, title="Ranging Beta");