Poisson Processes

The objective of this lab is to explore basic properties and behavior of Poisson processes.

Poisson process simulations

  1. Run several simulations using different values of the rate parameter \(\lambda\). Visually compare the density function \(f(t)=\lambda e^{-\lambda t}\) (plotted as a solid curve) to the histogram for the simulated interarrival times. You should find that the average number of arrivals per unit time during the simulation is close to the expected value \(\lambda\) (using the cumulative sum plot).
lambda <- 2 # rate parameter
maxcustomers <- 1000
T <-rexp(maxcustomers,rate=lambda) # samples interarrival times from exponential distributon
hist(T,xlab="Interarrival time T",
     freq=FALSE,breaks=50,main="Distribution of interarrival times") # plots density
t <- seq(0,10/lambda,length.out=50)
points(t,lambda*exp(-lambda*t),type="l") # plots original density function on top of histogram

  1. The step function plot showing number of customers who have arrived by time \(t\) should look roughly like a line with what slope?

Real life example

  1. Download the file CoalMineAccidents.csv from the course webpage and then load into RStudio. The setwd command sets your working directory, which needs to be the folder where you saved CoalMineAccidents.csv. The following code loads the data into RStudio as a data frame:
setwd("~/Documents/WebPage/Math365Spring2016/RmarkdownFiles") # update this to location of file
f <- read.table(file="CoalMineAccidents.csv", sep=",", header=TRUE)

This file has data on days between mining accidents in the UK. The first listed accident occurred on March 15, 1851, and the last one occurred on January 12, 1918. We can model this data as a Poisson process.

  1. Plot a histogram of the days between mining accidents to see if it looks roughly exponential:
hist(f[,"Days"], freq=FALSE) # density plot

  1. What would you estimate for \(\lambda\)? A good estimate can be made by simply calculating the average number of accidents per day from the data. Add the density function for your \(\lambda\) estimate to the histogram to see how well it seems to fit the data.

  2. Plot the cumulative sum of the ``Days" data. Does this look like a time-homogeneous process?

  3. To improve the rate parameter estimate, break the data into two time periods, splitting at the point where the rate parameter appears to have significantly changed. Estimate \(\lambda\) for each time period. Below is how to specify a range like the first 100 accidents:

f[1:100,"Days"]
##   [1]   0 157 123   2 124  12   4  10 216  80  12  33  66 232 826  40  12
##  [18]  29 190  97  65 186  23  92 197 431  16 154  95  25  19  78 202  36
##  [35] 110 276  16  88 225  53  17 538 187  34 101  41 139  42   1 250  80
##  [52]   3 324  56  31  96  70  41  93  24  91 143  16  27 144  45   6 208
##  [69]  29 112  43 193 134 420  95 125  34 127 218   2   0 378  36  15  31
##  [86] 215  11 137   4  15  72  96 124  50 120 203 176  55  93  59