The objective of this lab is to explore basic properties and behavior of Poisson processes.
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
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.
hist(f[,"Days"], freq=FALSE) # density plot
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.
Plot the cumulative sum of the ``Days" data. Does this look like a time-homogeneous process?
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