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.

$\begin{equation*} \Gamma(t) = \int_0^{\infty} u^{t-1} e^{-u} du \end{equation*}$

🔗

Theorem 8.4.2. Gamma Function on the natural numbers.

For $n \in \mathbb{N}\text{,}$

$\begin{equation*} \Gamma(n+1) = n! \end{equation*}$

🔗

Proof.

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

\begin{align*} \Gamma(n+1) & = \int_0^{\infty} u^n e^{-u} du\\ & = -u^n \cdot e^{-u} \big |_0^{\infty} + n \int_0^{\infty} u^{n-1} e^{-u} du \\ & = 0 - 0 + n \Gamma(n) \end{align*}

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{align*} F(x) & = P(X \le x)\\ & = 1 - P(X \gt x)\\ & = 1 - P(\text{fewer than r successes in [0,x]})\\ & = 1 - \big [ \frac{(\lambda x)^0 e^{-\lambda x}}{0!} + \frac{(\lambda x)^1 e^{-\lambda x}}{1!} + ... + \frac{(\lambda x)^{r-1} e^{-\lambda x}}{(r-1)!} \big ]\\ & = 1 - \sum_{k=0}^{r-1} \frac{(\lambda x)^k e^{-\lambda x}}{k!} \end{align*}$

🔗

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{align*} F'(x) & = \lambda e^{-\lambda x}/0!\\ & - \lambda e^{-\lambda x}/1! + \lambda x \cdot \lambda e^{-\lambda x}/1!\\ & - \lambda^2 2x e^{-\lambda x}/2! + \lambda^2 x^2 \cdot \lambda e^{-\lambda x}/2!\\ & - \lambda^3 3x^2 e^{-\lambda x}/3! + \lambda^3 x^3 \cdot \lambda e^{-\lambda x}/3!\\ & . . .\\ & - \lambda^{r-1} (r-1)x^{r-2} e^{-\lambda x}/(r-1)! + \lambda^{r-1} x^{r-1} \cdot \lambda e^{-\lambda x}/(r-1)!\\ & = \lambda^r x^{r-1} e^{-\lambda x}/(r-1)! \end{align*}$

🔗

where you can replace

$(r-1)! = \Gamma(r)\text{.}$

🔗

Notice that for this random variable,

$\mu = \lambda T$ can be obtained for the exponential distribution. For the Gamma distribution, the following takes

$\mu$ 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.

If X measures the interval until the rth success and $\mu$ as the average interval until the 1st success, then X with probability function

$\begin{equation*} f(x) = \frac{x^{r-1} \cdot e^{-\frac{x}{\mu}}}{\Gamma(r) \cdot \mu^r} \end{equation*}$

has a Gamma Distribution.

🔗

Checkpoint 8.4.4.

🔗

Example 8.4.5. 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

$\begin{equation*} f(x) = \frac{x^{5-1} \cdot e^{- \frac{x}{0.01}}}{\Gamma(5) \cdot 0.01^5} \end{equation*}$

The question above asks for

$\begin{equation*} P(X \gt 4 \text{seconds}) = P(X \gt \frac{4}{60} ) = 1 - F(\frac{4}{60}). \end{equation*}$

Therefore

$\begin{equation*} P(X \gt 4 \text{seconds}) = 1 - F(\frac{4}{60}) \approx 0.205627. \end{equation*}$

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 $\Gamma(5) = 4! = 24\text{.}$ You can compute the needed integral using the interactive cell immediately below.

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

🔗

Theorem 8.4.6. Verify Gamma Probability function.

$\begin{equation*} \int_0^{\infty} \frac{x^{r-1} e^{-x/ \mu}}{\Gamma(r) \mu^r} dx = 1 \end{equation*}$

🔗

Proof.

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

xxxxxxxxxx
 
# Gamma Distribution
var('x,mu,r')
assume(mu>0)
assume(r,'integer')
assume(r>1)
f(x) = x^(r-1)*e^(-x/mu)/(gamma(r)*mu^r)
S = integral(f,x,0,oo).full_simplify()
F = ("$ \\int_0^{\\infty} \\frac{x^{r-1}"+
    "e^{-x/ \\mu}}{\\Gamma(r) \\mu^r} dx = %s $"%str(S))
show(html(F))

🔗

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

🔗

Theorem 8.4.7. Properties of the Gamma Distribution.

For the gamma distribution from an underlying Poisson process with mean $\lambda\text{,}$

$\begin{equation*} \mu = r \lambda \end{equation*}$

$\begin{equation*} \sigma^2 = r \lambda^2 \end{equation*}$

$\begin{equation*} \gamma_1 = \frac{2}{\sqrt{r}} \end{equation*}$

$\begin{equation*} \gamma_2 = \frac{6}{r} + 3 \end{equation*}$

🔗

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
 
# Gamma Distribution
var('x,mu,r,alpha')
assume(mu>0)
assume(alpha,'integer')
assume(alpha>1)
@interact
def _(r=[2,3,6,9,alpha]):
    f(x) =x^(r-1)*e^(-x/mu)/(gamma(r)*mu^r)
    mean = integral(x*f,x,0,oo).full_simplify()
    M2 = integral(x^2*f,x,0,oo).full_simplify()
    M3 = integral(x^3*f,x,0,oo).full_simplify()
    M4 = integral(x^4*f,x,0,oo).full_simplify()
    
    pretty_print('Mean = ',mean)
    
    v = (M2-mean^2).full_simplify().factor()
    pretty_print('Variance = ',v)
    stand = sqrt(v)
    
    sk = (((M3 - 3*M2*mean + 2*mean^3))/stand^3).full_simplify()
    pretty_print('Skewness = ',sk)
    
    kurt = (M4 - 4*M3*mean + 6*M2*mean^2 -3*mean^4).factor()/stand^4
    pretty_print('Kurtosis = ',(kurt-3).factor(),'+3')

🔗

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
 
# Gamma Distribution Calculator
var('x,mu,r')
pretty_print('Enter the number of successes r desired,')
pretty_print('the given mean, and the value of X to get F(X)')
@interact(layout=dict(top=[['mu','b']],bottom=[['r']]))
def _(r=slider(1,10,1,2, label="$$ r = $$"),
      mu = input_box(2,label="$$\\mu = $$",width=10),
      b=input_box(2,label="$$ X = $$",width=10)):
    f(x) =x^(r-1)*e^(-x/mu)/(gamma(r)*mu^r)
    p = integral(f,x,0,b)
    
    pretty_print(html('$$\\text{Probability} = %s'%str(latex(p))
                      +' \\approx %s$$'%str(p.n(digits=5))))

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