R Functions dgamma, pgamma, qgamma, and rgamma

Random variable \(X\) is distributed \(X \sim gamma(\alpha, \theta)\) with mean \(\mu=\alpha \theta\) and variance \(\sigma^2 = \alpha \theta^2\) if \(X\) is the interval until the \(\alpha^{th}\) successful event when the average interval is \(\theta\). The probability of an interval of \(X=x\) \(f(X=n) =\frac{{1}}{{\theta^2\Gamma(\alpha)}}x^{\alpha-1}e^{-x/\theta}\). The probability of an interval between \(a\) and \(b\) is \(F_X(x) =\frac{{1}}{{\theta^2\Gamma(\alpha)}}x^{\alpha-1}e^{-x/\theta}dx\).

R function dgamma(x, reate) is the probability of interval x until the \(\alpha^{th}\) successful event when successful events occur with rate =\(1/\theta\). R function pgamma(q, shape, scale, lower.tail) is the cumulative probability (lower.tail = TRUE for left tail, lower.tail = FALSE for right tail) of less than or equal to an interval of q prior to the \(\alpha^{th}\) success. R function rgamma(n, shape, scale) returns n random numbers from the gamma distribution X~gamma(alpha, theta). R function qgamma(p, shape, scale, lower.tail) is the interval at the qth percentile (lower.tail = TRUE).

Example

On average, someone sends a money order once per 15 minutes (\(\theta = .25\)). What is the probability someone sends \(\alpha = 10\)* money orders in less than \(x=3\) hours?*

alpha = 10
theta = 15 / 60
x = 3
# exact
pgamma(q = x, shape = alpha, scale = theta)
## [1] 0.7576078
# simulated
mean(rgamma(n = 10000, shape = alpha, scale = theta) <= x)
## [1] 0.7527
library(dplyr)
library(ggplot2)

data.frame(x = 0:1000 / 100, prob = pgamma(q = 0:1000 / 100, shape = alpha, scale = theta, lower.tail = TRUE)) %>%
  mutate(Interval = ifelse(x >= 0 & x <= 3, "0 to 3", "other")) %>%
ggplot(aes(x = x, y = prob, fill = Interval)) +
  geom_area(alpha = 0.3) +
  labs(title = "X ~ Gam(alpha = 10, theta = .25)",
       subtitle = "Probability of 10 events in X hours when the mean time to an event is .25 hours.",
       x = "Interval (x)",
       y = "Cum Probability")