---
title: "Poisson Processes"
author: "Math 365"
date: "February 29, 2016"
output: html_document
---

### 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).

```{r}
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
```

(2) The plot showing number of customers who have arrived by time $t$ should look roughly like a line with what slope?

## Real life example

(3) 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:

```{r}
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. Note that we will only use the Days data, not the Casualties information. We can only track time between accidents using the Poisson process, not the severity of each particular accident.

(4) Plot a histogram of the days between mining accidents to see if it looks roughly exponential:

```{r}
hist(f[,"Days"], freq=FALSE) # density plot
```

(5) 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.

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

(7) 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:

```{r}
f[1:100,"Days"]
```