Gamma Distribution

Section 8.4 Gamma Distribution

Extending the exponential distribution model developed above, consider a Poisson Process where you start with an interval of variable length X so that X measures the interval needed in order to obtain the rth success for some natural number r. Then

R = (0, \infty)

and the resulting distribution of X will be called a Gamma distribution.

🔗

Definition 8.4.1. Gamma Function.

🔗

Γ (t) = \int_{0}^{\infty} u^{t - 1} e^{- u} d u

🔗

Theorem 8.4.2. Gamma Function on the natural numbers.

🔗

For

n \in N,

Γ (n + 1) = n!

🔗

Proof.

Letting n be a natural number and applying integration by parts one time gives

\begin{aligned} Γ (n + 1) & = \int_{0}^{\infty} u^{n} e^{- u} d u \\ = - u^{n} \cdot e^{- u} |_{0}^{\infty} + n \int_{0}^{\infty} u^{n - 1} e^{- u} d u \\ = 0 - 0 + n Γ (n) \end{aligned}

Continuing using an inductive argument to obtain the final result.

🔗

To find the probability function for the gamma distribution, once again focus on the development of F(x). Assuming r is a natural number greater than 1 and noting that X measures the interval length needed in order to achieve the rth success

\begin{aligned} F (x) & = P (X \leq x) \\ = 1 - P (X > x) \\ = 1 - P (fewer than r successes in [0,x]) \\ = 1 - [\frac{(λ x)^{0} e^{- λ x}}{0!} + \frac{(λ x)^{1} e^{- λ x}}{1!} + . . . + \frac{(λ x)^{r - 1} e^{- λ x}}{(r - 1)!}] \\ = 1 - \sum_{k = 0}^{r - 1} \frac{(λ x)^{k} e^{- λ x}}{k!} \end{aligned}

🔗

where the discrete Poisson probability function is used on the interval [0,x]. The derivative of this function however is "telescoping" and terms cancel. Indeed,

\begin{aligned} F^{'} (x) & = λ e^{- λ x} / 0! \\ - λ e^{- λ x} / 1! + λ x \cdot λ e^{- λ x} / 1! \\ - λ^{2} 2 x e^{- λ x} / 2! + λ^{2} x^{2} \cdot λ e^{- λ x} / 2! \\ - λ^{3} 3 x^{2} e^{- λ x} / 3! + λ^{3} x^{3} \cdot λ e^{- λ x} / 3! \\ . . . \\ - λ^{r - 1} (r - 1) x^{r - 2} e^{- λ x} / (r - 1)! + λ^{r - 1} x^{r - 1} \cdot λ e^{- λ x} / (r - 1)! \\ = λ^{r} x^{r - 1} e^{- λ x} / (r - 1)! \end{aligned}

🔗

where you can replace

(r - 1)! = Γ (r) .

🔗

Notice that for this random variable,

μ = λ T

can be obtained for the exponential distribution. For the Gamma distribution, the following takes

μ

to be the average interval till the first success and then modifies the corresponding Gamma parameters according to increasing values of r.

🔗

So, to summarize, you get the Gamma Probability Function below.

🔗

Definition 8.4.3. Gamma Probability Function.

🔗

X

is a random variable with real parameters

λ

and

α

with probability function

f (x) = \frac{x^{α - 1} \cdot e^{- \frac{x}{λ}}}{Γ (α) \cdot λ^{α}}

🔗

then

X

is said to have a Gamma Distribution.

🔗

The Gamma Distribution is a general one that has some important special cases. We will analyze the Gamma Distributions properties in general but a special case of the Gamma Distribution that extends the exponential distribution is especially appropriate at this point. Indeed, the distribution below considers the case when

α

is an integer.

🔗

Definition 8.4.4. Weibull Probability Function.

🔗

Suppose X measures the interval until the rth success and with parameter

λ

as the average interval until the 1st success. Then X with probability function

f (x) = \frac{x^{r - 1} \cdot e^{- \frac{x}{λ}}}{Γ (r) \cdot λ^{r}}

🔗

has a Weibul Distribution. (Note, this is a special case of the Gamma Distribution and is sometimes called the Gamma Distribution.)

🔗

Checkpoint 8.4.5.

It can be presumed that the life span in years of a certain brand of lightbulbs follows the Poisson process assumptions. Suppose that the mean life span till failure is known to be 0.7 years. To do long-term study, a series of light bulbs are arranged so that when the first one fails, the next one comes on, etc.

Determine the probability that it takes at most 3.5 years before the 4th bulb fails:

Answer.

0.73497

Solution.

This is a gamma distribution with r = 4 and

F (x) = 1 - \sum_{k = 0}^{3} \frac{(x / 0.7)^{k} e^{- x / 0.7}}{k!}

So,

\begin{aligned} P (X < 3.5) & = F (3.5) \\ = 1 - \frac{(3.5 / 0.7)^{0} e^{- 3.5 / 0.7}}{0!} \\ - \frac{(3.5 / 0.7)^{1} e^{- 3.5 / 0.7}}{1!} \\ - \frac{(3.5 / 0.7)^{2} e^{- 3.5 / 0.7}}{2!} - \frac{(3.5 / 0.7)^{3} e^{- 3.5 / 0.7}}{3!} \\ = 1 - e^{- 5} - 5 \cdot e^{- 5} - \frac{25}{2} \cdot e^{- 5} - \frac{125}{6} \cdot e^{- 5} \\ = 1 - e^{- 5} (1 + 5 + 25 / 2 + 125 / 6) \approx 0.73497 . \end{aligned}

🔗

Example 8.4.6. Router Requests Revisited Again.

For the third time, let’s consider a router which, over time, has been shown to receive on average 1000 requests in any given 10 minute period during regular working hours and you want to know the likelihood that it takes more than 4 seconds in order to receive the 5th request. As you have already seen, it takes on average

\frac{10}{1000} = \frac{1}{100} = 0.01

minutes to receive the first request so we use that again here. If X were to measure the time interval until the fifth actual request comes in, then the Gamma distribution would be a good model using

f (x) = \frac{x^{5 - 1} \cdot e^{- \frac{x}{0.01}}}{Γ (5) \cdot {0.01}^{5}}

The question above asks for

P (X > 4 seconds) = P (X > \frac{4}{60}) = 1 - F (\frac{4}{60}) .

Therefore

P (X > 4 seconds) = 1 - F (\frac{4}{60}) \approx 0.205627 .

Again, since X is a continuous variable you must integrate to compute probabilities. This will require integration by parts or you can use the F(x) from the derivation above. Here, let’s just let Sage do the integration for us noting that

Γ (5) = 4! = 24 .

You can compute the needed integral using the interactive cell immediately below.


    
        
xxxxxxxxxx
 
1
var('x')
2
f = x^4*e^(-x/0.01)/(24*0.01^5)
3
prob_complement = integrate(f,x,0,4/60)
4
print('For the interval from 0 to 4/60, the probability above is')
5
print(n(1-prob_complement))

    
    
    
    
        
            
                Language:
                
            
        
    
    




    
    
        
        Messages

🔗

Theorem 8.4.7. Verify Gamma Probability function.

🔗

\int_{0}^{\infty} \frac{x^{r - 1} e^{- x / λ}}{Γ (r) λ^{r}} d x = 1

🔗

Proof.

This is kind of a tough integral to do. Let’s just evaluate the sage code below.


    
        
xxxxxxxxxx
 
1
# Gamma Distribution
2
var('x,lam,r')
3
assume(lam>0)
4
assume(r,'integer')
5
assume(r>1)
6
f(x) = x^(r-1)*e^(-x/lam)/(gamma(r)*lam^r)
7
S = integral(f,x,0,oo).full_simplify()
8
F = ("$ \\int_0^{\\infty} \\frac{x^{r-1}"+
9
    "e^{-x/ \\lambda}}{\\Gamma(r) \\lambda^r} dx = %s $"%str(S))
10
show(html(F))

    
    
    
    
        
            
                Language:
                
            
        
    
    




    
    
        
        Messages

🔗

You can graph the gamma distribution’s probability function for various parameters below. Notice as r increases the curve becomes increasingly bell-shaped whicle changing the mean only shifts the curve around.


    
        
xxxxxxxxxx
 
1
# Gamma Distribution Graphing
2
var('x,lam,r')
3
assume(lam>0)
4
assume(r,'integer')
5
@interact
6
def _(r=[2,3,6,12,24],mu=slider(1,12,1,5,label='$$\\lambda$$')):
7
    f(x) =x^(r-1)*e^(-x/lam)/(gamma(r)*lam^r)
8
    plot(f,x,0,r*15).show(figsize=[5,3])

    
    
    
    
        
            
                Language:
                
            
        
    
    




    
    
        
        Messages

🔗

Theorem 8.4.8. Properties of the Gamma Distribution.

🔗

For the gamma distribution from an underlying Poisson process with

λ

found by solving the Poisson mean =

λ T,

μ = r λ

σ^{2} = r λ^{2}

γ_{1} = \frac{2}{\sqrt{r}}

γ_{2} = \frac{6}{r} + 3

🔗

Proof.

Evaluate the interactive cell below to let Sage do the heavy lifting.

🔗

Derivation of mean, variance, skewness, and kurtosis. Pick "alpha" for the general formulas.


    
        
xxxxxxxxxx
 
1
# Gamma Distribution
2
var('x,lam,r,alpha')
3
assume(lam>0)
4
assume(alpha,'integer')
5
assume(alpha>1)
6
@interact
7
def _(r=[2,3,6,9,alpha]):
8
    f(x) =x^(r-1)*e^(-x/lam)/(gamma(r)*lam^r)
9
    mean = integral(x*f,x,0,oo).full_simplify()
10
    M2 = integral(x^2*f,x,0,oo).full_simplify()
11
    M3 = integral(x^3*f,x,0,oo).full_simplify()
12
    M4 = integral(x^4*f,x,0,oo).full_simplify()
13
    
14
    print('Mean = ',mean)
15
    
16
    v = (M2-mean^2).full_simplify().factor()
17
    print('Variance = ',v)
18
    stand = sqrt(v)
19
    
20
    sk = (((M3 - 3*M2*mean + 2*mean^3))/stand^3).full_simplify()
21
    print('Skewness = ',sk)
22
    
23
    kurt = (M4 - 4*M3*mean + 6*M2*mean^2 -3*mean^4).factor()/stand^4
24
    print('Kurtosis = ',(kurt-3).factor(),'+3')

    
    
    
    
        
            
                Language:
                
            
        
    
    




    
    
        
        Messages

🔗

The interactive cell below can be used to compute the distribution function for the gamma distribution for various input values. If you desire to let r get bigger than the slider allows, feel free to edit the cell above and evaluate again.


    
        
xxxxxxxxxx
 
1
# Gamma Distribution Calculator
2
var('x,lam,r')
3
pretty_print('Enter the number of successes r desired,')
4
pretty_print('the given lambda, and the value of X to get F(X)')
5
@interact(layout=dict(top=[['lambda','b']],bottom=[['r']]))
6
def _(r=slider(1,10,1,2, label="$$ r = $$"),
7
      lam = input_box(2,label="$$\\lambda = $$",width=10),
8
      b=input_box(2,label="$$ X = $$",width=10)):
9
    f(x) =x^(r-1)*e^(-x/lam)/(gamma(r)*lam^r)
10
    p = integral(f,x,0,b)
11
    
12
    pretty_print(html('$$\\text{Probability} = %s'%str(latex(p))
13
                      +' \\approx %s$$'%str(p.n(digits=5))))

    
    
    
    
        
            
                Language:
                
            
        
    
    




    
    
        
        Messages


    
        
xxxxxxxxxx
 
1
r=3                # the number of successes desired
2
mu1 = 3            # the mean till first must be given
3
mu = mu1*r
4
sdev = sqrt(r)*mu1  # the formula for the standard deviation
5
6
M = mu*3   # the space is infinite but we just go out 3 standard deviations
7
X = 0:M    # quantiles for the space R of the random variable 
8
9
Ppois <- function(x){dgamma(x, shape=r, scale=mu1 )}  # create the probability function over X
10
11
curve(Ppois, from=0, to=M, xlab="X", col="blue", lwd=3,
12
 main="Gamma Sampling vs Gamma Curve vs Approximating 'Bell Curve'") 
13
Pnormal <- function(X){dnorm(X, mean=mu, sd=sdev)}   # to overlap a bell curve
14
curve(Pnormal, col="red", lwd=2, add=TRUE) 
15
16
Psample = rgamma(10^6, shape=r, scale=mu1)  # to create a histogram, sample a lot
17
# Xtop=max(Psample)          # for scaling the x-axis. Shift by 1/2 below.
18
hist(Psample, prob=TRUE, add=TRUE)

    
    
    
    
        
            
                Language:
                
            
        
    
    




    
    
        
        Messages

Essentials of Mathematical Probability and Statistics

Search Results:

Section 8.4 Gamma Distribution

Definition 8.4.1. Gamma Function.

Theorem 8.4.2. Gamma Function on the natural numbers.

Proof.

Definition 8.4.3. Gamma Probability Function.

Definition 8.4.4. Weibull Probability Function.

Checkpoint 8.4.5.

Example 8.4.6. Router Requests Revisited Again.

Theorem 8.4.7. Verify Gamma Probability function.

Proof.

Theorem 8.4.8. Properties of the Gamma Distribution.

Proof.