The Poisson Distribution

Scott Greenhalgh

2024-10-03
R has several built-in functions for the Poisson distribution. They’re listed in a table below
along with brief descriptions of what each one does.

We will begin our demo with 𝑟𝑝𝑜𝑖𝑠( ). First we’ll use it to make one (𝑛 = 1) randomly
generated observation of a random variable that follows the Poisson distribution and has
an average number of 10 successes (𝜆 = 10) per period.
rpois(1,10)

## [1] 3

This single observation isn’t very interesting on its own. So, the first somewhat interesting
thing we’ll do with 𝑟𝑝𝑜𝑖𝑠( ) is generate some data that we can use to plot this distribution
to see what it looks like. We’re going to generate 1,000 random observations with the same
value for 𝜆
[Link](2)

poisson_data <- [Link]('data' = rpois(1000, 10))

poisson_data %>% ggplot() +

geom_histogram(aes(x = data,
y = stat(count / sum(count))),
color = 'black',
binwidth = 1) +
geom_vline(xintercept = 10,
size = 1,
linetype = 'dashed',
color = 'red') +
theme_bw() +
labs(x = 'Number of successes per period',
y = 'Proportion',
title = '1,000 samples of Pois(lambda = 10)')
𝑑𝑝𝑜𝑖𝑠( ) and 𝑝𝑝𝑜𝑖𝑠( ) work the same way as their counterparts from the binomial
distribution. We’ll see them in action in the following practical examples.

Predicting the number of babies born in a hospital

Data from the maternity ward in a certain hospital shows that there is a historical average
of 4.5 babies born in this hospital every day Roberts et al. What is the probability that 6
babies will be born in this hospital tomorrow?
First, let’s calculate the theoretical probability of this event using 𝑑𝑝𝑜𝑖𝑠( ). The number of
successes we’re considering is 6, so we will set 𝑥 = 6. Additionally, this historical average
of 4.5 babies per day is our value for 𝜆, so we will set 𝜆 = 6,
dpois(6, 4.5)

## [1] 0.1281201

The theoretical probability of 6 babies being born tomorrow if the historical average is 4.5
is about 13%.
Now let’s simulate births in this hospital for a year (𝑛 = 365) using 𝑟𝑝𝑜𝑖𝑠( ) and compare
the proportion of days in which there were 6 births to the theoretical probability we
calculated above. We will also visualize this result.
[Link](2)

babies <- [Link]('data' = rpois(365, 4.5))

babies %>% ggplot() +
geom_histogram(aes(x = data,
y = stat(count / sum(count)),
fill = data == 6),
binwidth = 1,
color = 'black',) +
scale_x_continuous(breaks = 0:10) +
labs(x = 'Number of babies born per period',
y = 'Proportion',
title = '365 simulated births in a hospital with Pois(lambda = 4.5)')
+
theme_bw()

babies %>% dplyr::summarize(six_babies = sum(babies$data == 6) / n())

## six_babies
## 1 0.1150685

The simulated result of about 11.5% is pretty close to our theoretical probability of about
13%.
What about the probability of more than 6 babies being born?
ppois(6,4.5, [Link]=FALSE)

## [1] 0.1689494
This theoretical probability is about 16.9%. We’ll learn about what ‘[Link] = FALSE’
means a bit later in the course, specifically when we start to talk about cumulative
probability functions, i.e. 𝑃(𝑋 < 𝑥), instead of just 𝑃(𝑋 = 𝑥), which are called probability
density functions (or probability mass functions).
What about the corresponding proportion in our simulation?
babies %>% ggplot() +
geom_histogram(aes(x = data,
y = stat(count / sum(count)),
fill = data > 6),
binwidth = 1,
color = 'black',) +
scale_x_continuous(breaks = 0:10) +
labs(x = 'Number of babies born per period',
y = 'Proportion',
title = '365 simulated births in a hospital with Pois(lambda = 4.5)')
+
theme_bw()

babies %>% dplyr::summarize(six_babies = sum(babies$data > 6) / n())

## six_babies
## 1 0.1863014

The simulated proportion of about 18.6% is pretty close to the theoretical proportion
above.
Simulating deaths by horse kick of Prussian cavalry
One of the most famous studies based on the Poisson distribution was by Ladislaus
Bortkiewicz, a Polish economist and statistician, in his book The Law of Small Numbers.
This book actually contained two studies: one about deaths by horse kicks of Prussian
cavalry soldiers and one about child suicides in Prussia Dobrow et al. The former is far
better known than the latter, probably because its topic is far less grim.
In his study, Bortkiewicz considered 20 years of data for 10 corps (groups) of Prussian
cavalry soldiers. Over this period there were 122 total deaths by horse kick among these
soldiers. Bortkiewicz divided the data into 20 individual periods for each group of soldiers,
for a total of 20 x 10 = 200 corps years. The average number of deaths by horse kick was
121 / 200 = 0.61, which means that lambda = 0.61. The data from his study is shown below.
horsekicks<-data("VonBort")
head(horsekicks)

## [1] "VonBort"

fisherYESVonBort <- VonBort[VonBort$fisher == 'yes',] #fisher column just

indicates if it was used in a study by Fisher in 1925. For simplicity, we'll
select 'yes' entries
head(fisherYESVonBort)

## deaths year corps fisher

## 3 0 1875 II yes
## 4 0 1875 III yes
## 5 0 1875 IV yes
## 6 0 1875 V yes
## 8 1 1875 VII yes
## 9 1 1875 VIII yes

Next, we will compare the experimentally derived data collected by von Bortkiewicz to the
theoretical data predicted by the model of a Poisson process.
#defining a couple colors for plotting purposes
greenFill_VB <- rgb(0.2,0.9,0.2,0.7)
grayFill_VB <- rgb(0.2,0.2,0.2,0.7)

#modelling a poisson process only takes 1 parameter, the mean

lambda_hat_fYes <- mean(fisherYESVonBort$deaths)

xvals <- 0:4

nY <- length(fisherYESVonBort$deaths)
expected_fYes <- nY*dpois(xvals,lambda_hat_fYes )

#plot of just the 'yes' data

fYES_histo <- ggplot( [Link]( deaths=fisherYESVonBort$deaths ),
aes(x=deaths) ) +
geom_histogram(binwidth=.5, color="green", fill=greenFill_VB)
#the following combines the horse kick data w/ what is expected from the
Poisson dist.
fYes_histoData <- layer_data( fYES_histo )
data2fit_YES <- fYes_histoData$count[seq(1,length(fYes_histoData)/2+1,2)]
yesFisher_ObsGF <- [Link](data_subset = "observed",count = data2fit_YES,
xvals = 0:4)
yesFisher_ExpGF <- [Link](data_subset = "expected",count = expected_fYes,
xvals = 0:4)
[Link] <- rbind( yesFisher_ObsGF, yesFisher_ExpGF )
head([Link])

## data_subset count xvals

## 1 observed 109.0000 0
## 2 observed 65.0000 1
## 3 observed 22.0000 2
## 4 observed 3.0000 3
## 5 observed 1.0000 4
## 6 expected 108.6702 0

Using the merged data set, we can visually see if the horse kick deaths data resembles a
Poisson distribution by plotting a histogram of the death data beside the what was is
predicted by Poisson distribution with the same mean:
#plot data & format
Poisson_fY <- ggplot([Link], aes(x = xvals, y = count, fill =
data_subset)) +
geom_col(position = "dodge") +
ggtitle("Fisher = 'yes'") +
xlab("Annual death rate") +
ylab("Count") +
scale_color_manual(values = c("green")) +
scale_fill_manual(values = c(greenFill_VB, grayFill_VB)) +
annotate("text", x = 3, y = 90, label = " ", color = "green", size = 5)
print(Poisson_fY)
Comparing the expected and observed columns tells us the total number of deaths closely
follows a Poisson distribution with 𝜆 = 0.61.

Simulating costs of car accidents

The following question was taken from Probability in with Applications in R by Robert
Dobrow.
Suppose that the number of accidents per month at a busy intersection in the center of a
certain city is 7.5. This event follows a Poisson distribution and
lambda <- 7.5

Every time an accident occurs at this intersection, the city government has to pay about
$25,000 to clean up the area. What is the average cost of these accidents per year?
This question is a lot easier than it probably sounds. We know that the average number of
accidents per month is 7.5. We also know that there are 12 months in a year, so the average
number of accidents per year is just the product of these two numbers. Finally, since we
also know the average cost per accident, the average cost of accident clean-up per year for
this city is just the product of these three numbers.
7.5*12*25000

## [1] 2250000
In a typical year, this city can expect to pay about $2.25 million in accident clean-up costs
for this intersection.
One thing that we should remember, however, is that we are talking about a random
variable which follows a certain distribution. This means that there will always be some
random variation in annual accident costs for this city. To get an idea of how much accident
costs can vary, we’re going to run a simulation.
First, we’ll simulate the annual accident cost for one year. We’ll use 𝑛 = 12 because 𝜆 = 7.5
represents the average number of accidents per month, and we want to simulate 12
months. We will save these results as a variable called accidents.
Next, we’ll multiply each element inside of accidents by 25,000 in order to calculate the
average cost per month of accident clean-up. Finally, we will add these monthly costs
together to get the total annual cost.
[Link](2)

accidents <- rpois(12, 7.5)

cost <- sum(25000 * accidents)

cost

## [1] 2325000

In this simulation, the total cost of cleaning up after accidents was about $2.32 million.
Now we’re going to use replicate() to simulate accident costs at this intersection for 1,000
years. We’re using a high number so that we can get a good look at what this distribution
looks like.
[Link](2)

cost_sim <- [Link]('data' = replicate(1000, sum(25000 * rpois(12, 7.5))))

mean(cost_sim$data)

## [1] 2259225

In this simulation, the mean cost of accident clean-up is about $2.26 million, which is quite
close to the theoretical total. This cost is marked on the histogram below with a dashed red
line.
cost_sim %>% ggplot() +
geom_histogram(aes(x = data,
y = stat(count / sum(count))),
color = 'black',
binwidth = 100000) +
geom_vline(xintercept = mean(cost_sim$data),
color = 'red',
size = 1,
 linetype = 'dashed') +
scale_x_continuous(breaks = c(1500000,
2000000,
2259225,
2500000,
2500000,
3000000)) +
labs(x = 'Annual accident cost',
y = 'Proportion',
title = 'Distribution of 1,000 simulated years of car accident costs',
subtitle = 'Pois(lambda = 7.5)') +
theme_bw()

How big is this

town? Is this the most dangerous intersection in terms of accident frequency? How does it
compare to others in the town? Unfortunately we can’t answer any of these questions. But
spending over $2 million in a typical year to deal with accidents at a single intersection is a
sign that something needs to be done about the design of that intersection to decrease the
frequency of those accidents because at the very least, they are a drain on the city’s
finances.

Bomb hits over London during WWII

The following setting is very general. Suppose n balls are thrown into 𝑛/𝜆 bowls so that
each ball has an equal chance of landing in any bowl. If a ball lands in a bowl, call it a “hit.”
The chance that a ball hits a particular bowl is 1/(𝑛/𝜆) = 𝜆/𝑛. Keeping track of whether or
not each ball hits that bowl, the successive hits form a Bernoulli sequence, and the number
of hits has a binomial distribution with parameters 𝑛 and 𝜆/𝑛. If n is large, the number of
balls in each bowl is approximated by a Poisson distribution with parameter 𝑛 ∗ (𝜆/𝑛) = 𝜆.
Many diverse applications can be fit into this ball and bowl setting. In his classic analysis of
Nazi bombing raids on London during World War II, William Feller (1968) modeled bomb
hits (balls) using a Poisson distribution. The city was divided into 576 small areas (bowls)
of 1/4 km squared. The number of areas hit exactly 𝑘 times was counted. There were a
total of 537 hits, so the average number of hits per area was 537/576 = 0.9323.
The data from this study is shown in the table below.
library(readxl)
LondonBombingWWIIdata <- read_excel("[Link]")
head(LondonBombingWWIIdata)

## # A tibble: 6 × 3
## hits observed_hits expected_hits
## <chr> <dbl> <dbl>
## 1 0 229 227.
## 2 1 211 211.
## 3 2 93 98.6
## 4 3 35 30.7
## 5 4 7 7.1
## 6 5+ 1 1.6

This table is interpreted in a way similar to the one about horse kicking deaths. The first
column represents the number of balls (bombs) that landed in one of the bowls (1/4 km
square areas). The second and third columns represent the number of areas in which a
certain number of bombs landed and the number of areas in which it was expected that
number of bombs would land according to the Poisson distribution, respectively. Like the
study about horse kick deaths, the observed and expected values are quite close because
the random variable being examined appears to follow a Poisson distribution.
Now we’re going to run a simulation with this data that’s based on one by Robert Dobrow.
The first thing we’re going to do is create some variables which match the ones described
in the quotes from Dobrow at the beginning of this section. n will represent the number of
bombs (balls) dropped on this section of London and u will represent the number of 1/4
square kilometer sections of the city that were subject to bombing (bowls). We will assume
n is sufficiently large to define 𝜆 = 𝑛/𝑢. The bowls variable is a vector of length 𝑢 in which
each element (bowl) represents one quarter square kilometer section of the city that was
subject to bombing. Each element (bowl) is initialized with a value of 0 because before the
bombing starts, 0 bombs have landed in each section (bowl).
n <- 537
u <- 576
lambda = n / u
bowls <- rep(0, u)
We will use a for loop for this simulation. The loop iterates through a sequence of numbers
from 1 to 𝑛 = 537, once for each bomb that was dropped in the section of London that was
targeted. The first line of the loop selects one number at random from a sequence of
numbers from 1 to 𝑢 = 576. (Remember that all of these numbers have starting values of 0.)
The number that is selected represents the section of the city where the bomb will land.
The second line of the loop increases the number of bombs that have landed on this part of
the city by 1 to indicate that a bomb has landed there during this round of iteration. The
process runs a total of 537 times, once for each bomb that was dropped.
for (i in 1:n){
i <- sample(1:u, 1)
bowls[i] <- bowls[i] + 1
}

Now let’s have a glance at our results. Imagine that the printout below is an accurate spatial
representation of this section of London. Each one of the numbers below represents one
quarter square kilometer section of the area that was targeted and how many bombs
landed there. Some places were luckier than others.
head(bowls)

## [1] 0 0 3 0 1 0

We can briefly summarize this data using

table(bowls)

## bowls
## 0 1 2 3 4 5
## 221 219 105 21 5 5

How do these simulated totals compared to what we would expect according to the Poisson
distribution? Getting all of this data into a summary dataframe will be somewhat
complicated, but the process is mostly familiar.
bombing_sim <- [Link]('sim_data' = bowls)

bombing_sim_summary <- [Link]('bombs_landed' = c('0', '1', '2', '3', '4',

'5+'))

for (col in colnames(bombing_sim)){

zero <- sum(bombing_sim[[col]] == 0)
one <- sum(bombing_sim[[col]] == 1)
two <- sum(bombing_sim[[col]] == 2)
three <- sum(bombing_sim[[col]] == 3)
four <- sum(bombing_sim[[col]] == 4)
five_plus <- sum(bombing_sim[[col]] >= 5)

column <- [Link]('simulated_total' = c(zero,

one,
 two,
three,
four,
five_plus))

bombing_sim_summary <- cbind(bombing_sim_summary, column)

}

# calculation of theoretical totals

zero_to_four <- dpois(0:4, lambda) * u
five_plus <- ppois(4, lambda, [Link] = FALSE) * u

theoretical_total <- [Link]('theoretical_total' = c(zero_to_four,

five_plus))

bombing_sim_summary <- cbind(bombing_sim_summary, round(theoretical_total,

2))

bombing_sim_summary

## bombs_landed simulated_total theoretical_total

## 1 0 221 226.74
## 2 1 219 211.39
## 3 2 105 98.54
## 4 3 21 30.62
## 5 4 5 7.14
## 6 5+ 5 1.57

Notice how the theoretical probabilites were calculated. dpois() was used for the first five,
but the last one required ppois().
We’ll conclude with a visual summary of our bombing simulation results. The red dashed
line is drawn at our value for 𝜆.
bombing_sim %>% ggplot() +
geom_bar(aes(bowls),
color = 'black') +
geom_vline(xintercept = lambda,
color = 'red',
linetype = 'dashed',
size = 1) +
scale_x_continuous(breaks = 0:max(bowls)) +
labs(x = 'Bombs landed (x)',
y = 'Number of sections in which x bombs landed',
title = 'Bombing simulation results with lambda = 0.932') +
theme_bw()