Gamma Distribution Probabilities using R

Gamma Distribution Probabilities using R

In this tutorial, you will learn about how to use dgamma(), pgamma(), qgamma() and rgamma() functions in R programming language to compute the individual probabilities, cumulative probabilities, quantiles and to generate random sample for Gamma distribution.

Before we discuss R functions for Gamma distribution, let us see what is Gamma distribution.

Gamma Distribution

Gamma distribution distribution is a continuous type probability distribution. Gamma distribution has found applications in many fields. For example, gamma distribution is suitable for describing waiting times between successive occurrences of a random events, survival times, system reliability, etc.

Let $X\sim G(\alpha,\beta)$. Then the probability distribution of $X$ is

$$ \begin{aligned} f(x)&= \begin{cases} \frac{1}{\beta^\alpha\Gamma(\alpha)}x^{\alpha -1}e^{-x/\beta}, & x > 0;\alpha, \beta > 0; \\ 0, & Otherwise. \end{cases} \end{aligned} $$

where $\alpha$ is the shape parameter and $\beta$ is the scale parameter of Gamma distribution.

Read more about the theory and results of Gamma distribution here.

Gamma probabilities using dgamma() function in R

For continuous probability distribution, density is the value of the probability density function at $x$ (i.e., $f(x)$).

The syntax to compute the probability density function for Gamma distribution using R is

dgamma(x,shape, rate=1, scale=1/rate)

where

• x : the value(s) of the variable and,
• shape : shape parameter of gamma distribution,
• rate : rate parameter of gamma distribution,
• scale : scale parameter of gamma distribution.

Note: Rate =1/scale is an alternative way to specify the scale parameter.

The dgamma() function gives the density for given value(s) x, shape and scale.

Numerical Problem for Gamma Distribution

To understand the four functions dgamma(), pgamma(), qgamma() and rgamma(), let us take the following numerical problem.

Gamma Distribution Example

The lifetime of certain equipment is described by a random variable $X$ whose distribution is Gamma with parameters $\alpha = 2$ and $\beta = 1/3$.

(a) Find the value of the density function at $x=2$.
(b) Plot the graph of Gamma probability distribution.
(c) Find the probability that the lifetime of equipment is at most 2 unit of time.
(d) Find the probability that the lifetime of equipment is at least 1 unit of time.
(e) Find the probability that the lifetime of equipment is less than 2.5 unit of time but greater than 1.5 unit of time.
(f) Plot the graph of cumulative Gamma probabilities.
(g) What is the value of $c$, if $P(X\leq c) \geq 0.70$?
(h) Simulate 1000 Gamma distributed random variables with $\alpha= 2$ and $\beta = 1/3$.

Let $X$ denote the lifetime of certain equipment. Given that $X\sim Gamma(\alpha=2, \beta=1/3)$.

Example 1: How to use dgamma() function in R?

To find the value of the density function at $x=2$ we need to use dgamma() function.

Let $X$ denote the lifetime of certain equipment. Here $X\sim Gamma(2,1/3)$.

First let us define the given parameters as

# shape parameter 
alpha <- 2
# scale parameter
beta <- 1/3

The probability density function of $X$ is

$$ \begin{aligned} f(x)&= 9xe^{-3x},\\ &\quad\text{for } x \geq 0. \end{aligned} $$

For part (a), we need to find the density function at $x=2$. That is $f(2)$.

(a) The value of the density function at $x=2$ is

$$ \begin{aligned} f(2)&= 9\times 2\times e^{-3\times 2}\\ &=18\times e^{-6}\\ &= 0.0446175 \end{aligned} $$

The above probability can be calculated using dgamma(0.35,shape=2,scale=1/3) function in R.

# Compute Gamma probability
result1 <- dgamma(2,shape=alpha,scale=beta)
result1

[1] 0.04461754

Example 2 Visualize Gamma probability distribution

Using dgamma() function we can compute Gamma distribution probabilities for given x, shape and scale. To plot the probability density function of Gamma distribution, we need to create a sequence of x values and compute the corresponding probabilities.

# create a sequence of x values
x <- seq(0,4, by=0.02)
## Compute the Gamma pdf for each x
px<- dgamma(x,shape=alpha,scale=beta)

(b) Visualizing Gamma Distribution with dgamma() function and plot() function in R:

The probability density function of Gamma distribution with given 2 and 0.3333333 can be visualized using plot() function as follows:

## Plot the Gamma probability dist
plot(x,px,type="l",xlim=c(0,4),ylim=c(0,max(px)),
     lwd=3, col="darkred",ylab="f(x)")
title("PDF of Gamma (alpha = 2, beta= 1/3)")

Gamma cumulative probability using pgamma() function in R

The syntax to compute the cumulative probability distribution function (CDF) for Gamma distribution using R is

pgamma(q,shape, rate=1, scale=1/rate)

where

• q : the value(s) of the variable,
• shape : shape parameter of gamma distribution,
• rate : rate parameter of gamma distribution,
• scale : scale parameter of gamma distribution.

Using this function one can calculate the cumulative distribution function of Gamma distribution for given value(s) of q (value of the variable x), shape and scale.

Example 3: How to use pgamma() function in R?

In the above example, for part (c), we need to find the probability $P(X\leq 2)$.

(c) The probability that the lifetime is at most 2 unit of time is

$$ \begin{aligned} P(X\leq 2) &=\int_0^{2} f(x)\; dx. \end{aligned} $$

## Compute cumulative Gamma probability
result2 <- pgamma(2,shape=alpha,scale=beta)
result2

[1] 0.9826487

Example 4: How to use pgamma() function in R?

In the above example, for part (d), we need to find the probability $P(X \geq 1)$.

To calculate the probability that a random variable $X$ is greater than a given number, one can use the option lower.tail=FALSE in pgamma() function.

Above probability can be calculated easily using pgamma() function with argument lower.tail=FALSE as

$P(X \geq 1) =\int_{1}^\infty f(x)\; dx$= pgamma(1,shape=alpha,scale=beta,lower.tail=FALSE)

or by using complementary event as

$P(X \geq 1) = 1- P(X\leq 1)$= 1- pgamma(1,shape=alpha,scale=beta)

# compute cumulative Gamma probabilities
# with lower.tail False
pgamma(1,shape=alpha,scale=beta,lower.tail=FALSE)

[1] 0.1991483

(d) The probability that the lifetime is at least 1 unit of time is

$$ \begin{aligned} P(X\geq 1) &=\int_{1}^\infty f(x)\; dx\\ &=0.1991483. \end{aligned} $$

# Using complementary event
1-pgamma(1,shape=alpha,scale=beta)

[1] 0.1991483

Example 5: How to use pgamma() function in R?

One can also use pgamma() function to calculate the probability that the random variable $X$ is between two values.

(e) The probability that the lifetime of equipment is less than 2.5 unit of time but greater than 1.5 unit of time can be written as $P(1.5 < X < 2.5)$.

$$ \begin{aligned} P(1.5 < X < 2.5) &= P(X< 2.5) -P(X < 1.5)\\ &= 0.9952988 - 0.9389005\\ &= 0.0563983 \end{aligned} $$

The above probability can be calculated using pgamma() function as follows:

a <- pgamma(2.5,shape=alpha,scale=beta)
b <- pgamma(1.5,shape=alpha,scale=beta)
result3 <- a - b
result3

[1] 0.05639826

Example 6: Visualize the cumulative Gamma probability distribution

Using pgamma() function we can compute Gamma cumulative probabilities (CDF) for given x, shape1 and shape2. To plot the CDF of Gamma distribution, we need to create a sequence of x values and compute the corresponding cumulative probabilities.

# create a sequence of x values
x <- seq(0,4, by=0.02)
## Compute the Gamma pdf for each x
Fx <- pgamma(x,shape=alpha,scale=beta)

(f) Visualizing Gamma Distribution with pgamma() function and plot() function in R:

The cumulative probability distribution of Gamma distribution with given x, shape1 and shape2 can be visualized using plot() function as follows:

## Plot the Gamma  probability dist
plot(x,Fx,type="l",xlim=c(0,4),ylim=c(0,1),
     lwd=3, col="darkred",ylab="F(x)")
title("CDF of Gamma (alpha = 2, beta= 1/3)")

Gamma Distribution Quantiles using qgamma() in R

The syntax to compute the quantiles of Gamma distribution using R is

qgamma(p,shape,rate=1,scale=1/rate)

where

• p : the value(s) of the probabilities,
• shape : shape parameter of gamma distribution,
• rate : rate parameter of gamma distribution,
• scale : scale parameter of gamma distribution.

The function qgamma(p,shape,scale) gives $100*p^{th}$ quantile of Gamma distribution for given value of p, shape and scale.

The $p^{th}$ quantile is the smallest value of Gamma random variable $X$ such that $P(X\leq x) \geq p$.

It is the inverse of pgamma() function. That is, inverse cumulative probability distribution function for Gamma distribution.

Example 7: How to use qgamma() function in R?

In part (g), we need to find the value of $c$ such a that $P(X\leq c) \geq 0.70$. That is we need to find the $70^{th}$ quantile of given Gamma distribution.

alpha <- 2
beta <- 1/3
prob <- 0.70

# compute the quantile for Gamma  dist
qgamma(0.70,shape=alpha, scale=beta)

[1] 0.8130722

The $70^{th}$ percentile of given Gamma distribution is 0.8130722.

Visualize the quantiles of gamma Distribution

The quantiles of gamma distribution with given p, shape=alpha and scale=beta can be visualized using plot() function as follows:

p <- seq(0,1,by=0.02)
qx <- qgamma(p,shape=alpha,scale=beta)
# Plot the Quantiles of Gamma  dist
plot(p,qx,type="l",lwd=2,col="darkred",
     ylab="quantiles",
main="Quantiles of Gamma(alpha= 2,beta = 1/3)")

Simulating Gamma random variable using rgamma() function in R

The general R function to generate random numbers from Gamma distribution is

rgamma(n,shape,rate=1,scale=1/rate)

where,

• n : the sample observations,
• shape : shape parameter of gamma distribution,
• rate : rate parameter of gamma distribution,
• scale : scale parameter of gamma distribution.

The function rgamma(n,shape,scale) generates n random numbers from Gamma distribution with given shape and scale.

Example 8: How to use rgamma() function in R?

In part (h), we need to generate 1000 random numbers from Gamma distribution with given $shape = 2$ and $scale=1/3$.

(h) We can use rgamma(1000,shape,scale) function to generate random numbers from Gamma distribution.

## initialize sample size to generate
n <- 1000
# Simulate 1000 values From Gamma  dist
x_sim <- rgamma(n,shape=alpha,scale=beta)

The below graphs shows the density of the simulated random variables from Gamma Distribution.

## Plot the simulated data
plot(density(x_sim),xlab="Simulated x",ylab="density",
     lwd=5,col="darkred",
     main="Simulated data from Gamma(2,1/3) dist")

If you use same function again, R will generate another set of random numbers from $Gamma(2,1/3)$.

# Simulate 1000 values From Gamma  dist
x_sim_2 <- rgamma(n,shape=alpha,scale=beta)

## Plot the simulated data
plot(density(x_sim_2),xlab="Simulated x",ylab="density",
     lwd=5,col="blue",
     main="Simulated data from Gamma(2,1/3) dist")

For the simulation purpose to reproduce same set of random numbers, one can use set.seed() function.

# set seed for reproducibility
set.seed(1457)
# Simulate 1000 values From Gamma  dist
x_sim_3 <- rgamma(n,shape=alpha,scale=beta)

## Plot the simulated data
plot(density(x_sim_3),xlab="Simulated x",ylab="density",
     lwd=5,col="darkred",
     main="Simulated data from Gamma(2,1/3) dist")

set.seed(1457)
# Simulate 1000 values From Gamma  dist
x_sim_4 <- rgamma(n,shape=alpha,scale=beta)

## Plot the simulated data
plot(density(x_sim_4),xlab="Simulated x",ylab="density",
     lwd=5,col="darkred",
     main="Simulated data from Gamma(2,1/3) dist")

Since we have used set.seed(1457) function, R will generate the same set of Gamma distributed random numbers.

To learn more about other discrete and continuous probability distributions using R, go through the following tutorials:

Discrete Distributions Using R

Binomial distribution in R
Poisson distribution in R
Geometric distribution in R
Negative Binomial distribution in R
Hypergeometric distribution in R

Continuous Distributions Using R

Uniform distribution in R
Exponential distribution in R
Log-Normal distribution in R
Normal distribution in R
Beta distribution in R
Cauchy distribution in R
Laplace distribution in R
Logistic distribution in R
Weibull distribution in R

Endnote

In this tutorial, you learned about how to compute the probabilities, cumulative probabilities and quantiles of Gamma distribution in R programming. You also learned about how to simulate a Gamma distribution using R programming.

To learn more about R code for discrete and continuous probability distributions, please refer to the following tutorials:

Probability Distributions using R

Let me know in the comments below, if you have any questions on Gamma Distribution using R and your thought on this article.