Let us apply the Poisson distribution in form of an example. We will look at a one-hundred-year flood, a concept that is widely used in river engineering to design flood control measures.
Recall the mathematical notation of a Poisson random variable:
\[ P(X = x) = e^{-\lambda}\frac{\lambda x}{x!}, \qquad x = 0, 1, 2, . . . ,\]
where \(\lambda\) is a positive real number, which represents the average number of events occurring during a fixed time interval, and \(e \approx 2.7182818\).
A 100-year flood is a shorthand term for a flood with an annual exceedance probability of 1 % and an average recurrence interval of 100 years. The term can be confusing as is might suggest that a given flood occurs once every 100-years. However, this is not true. A flood with an annual exceedance probability of 1 % indicates that with a probability of 0.01 a flood of a magnitude corresponding to a 100-year flood occurs in any given year.
Putting this in the context of a Poisson distribution the expected value, \(E(x) = \lambda\), of such a flood event during the fixed interval of 100 years is set to \(\lambda= 100\times 0.01 = 1\). Consequently, the Poisson random variable \(X\) is the number of occurrences, which of course may take several values depending on the question. We may be interested in the probability that such a flood event will not occur during the 100-year interval \(P(x=0)\). Or, we want to know the probability of such a flood event occurring exactly once during the 100-year interval, thus \(P(x=1)\). We might also want to know the probability that two or more of such flood events occur during the 100-year interval, thus \(P(x \ge 2)\). Plugging these values in the equation from above yields:
\(\lambda = 1, x = 0,1,2,..,n\)
\[ P(X = 0) = e^{-1}\frac{1 \times 0}{0!}, \qquad \text{for } x = 0\] \[ P(X = 1) = e^{-1}\frac{1 \times 1}{1!}, \qquad \text{for } x = 1\] \[ P(X \ge 2) = \sum_{i=2}^n e^{-1}\frac{1 \times x_i}{x_i!}, \qquad \text{for } x_i = 2,3,...,n\]
We turn to R to do the calculations. We will make use of the
dpois()
and the ppois()
functions.
lambda <- 1 # set expected value
x <- c(0, 1) # set Poisson random variables of interest
dpois(x = x[1], lambda = lambda)
## [1] 0.3678794
dpois(x = x[2], lambda = lambda)
## [1] 0.3678794
# recall: lower.tail = TRUE, probabilities are P[X <= x], otherwise P[X > x]
ppois(q = x[2], lambda = lambda, lower.tail = FALSE)
## [1] 0.2642411
The results indicate, that the probability of no flood \(P(X = 0)\) with a magnitude corresponding to a 100-year flood occurring during a period of 100 years is 0.37. Interestingly, this is as likely as the occurrence of exactly one flood \(P(X = 1)\). The probability that two or more \((P(X \ge 2))\) of such flood events occur during the 100-year interval is 0.26 and thus lower. However, please note that the probability that two or more \((P(X \ge 2))\) of such flood events occur during the 100-year interval is approximately 26 %!
For a sanity check we sum up the probabilities \(P(X = 0)\), \(P(X = 1)\) and \(P(X \ge 2)\), which should yield 1:
dpois(x = x[1], lambda = lambda) +
dpois(x = x[2], lambda = lambda) +
ppois(q = x[2], lambda = lambda, lower.tail = FALSE)
## [1] 1
In order to gain a better intuition we plot the probabilities of the Poisson random variable \(x = 0,1,2,3,4,\ge 5\).
lambda <- 1 # set expected value
out <- NULL
out[1] <- dpois(x = 0, lambda = lambda)
out[2] <- dpois(x = 1, lambda = lambda)
out[3] <- dpois(x = 2, lambda = lambda)
out[4] <- dpois(x = 3, lambda = lambda)
out[5] <- dpois(x = 4, lambda = lambda)
out[6] <- ppois(q = 5, lambda = lambda, lower.tail = FALSE)
# plot
barplot(out,
ylab = "Probability P(X = x)",
xlab = "Number of floods (x)",
names = paste(c("0", "1", "2", "3", "4", ">= 5")),
main = 'Probability of the occurence of "100-year floods" \nduring a period of 100 years'
)
Exercise: At the bus stop Emmichstrasse near FU Campus Lankwitz the transport company schedules 6 buses per hour. Thus, we may regard 6 buses per hour as the expected value lambda. Let us look at this declared intention as a Poisson process and calculate the probability for X = 0, X = 1, …, X = 12 buses/h. Provide a graph for the results!
### your code here
lambda <- 6
prob_buses <- NA
for (x in 0:12) {
prob_buses[x + 1] <- dpois(x = x, lambda = lambda)
}
barplot(prob_buses,
ylab = "Probability P(X = x)", xlab = "Number of buses per hour (x)",
col = "#f7f90e", names = 0:12,
main = "Probability for the expected number of buses\n at the bus stop Emmichstrasse"
)
Citation
The E-Learning project SOGA-R was developed at the Department of Earth Sciences by Kai Hartmann, Joachim Krois and Annette Rudolph. You can reach us via mail by soga[at]zedat.fu-berlin.de.
Please cite as follow: Hartmann, K., Krois, J., Rudolph, A. (2023): Statistics and Geodata Analysis using R (SOGA-R). Department of Earth Sciences, Freie Universitaet Berlin.