# Probability of a run of events from past data

I would like to calculate the probability that there will be a car crash at a certain crossroad for three days in a row. I would like to compute that based on the data I have about the crashes in the same spot.

My data consists of 8 years of records. For each day we know whether there has been one or more crashes in a particular crossroad. The data looks like this:

NNNNNYN NNNYNYN YNNYYNN NNYYNYN NYNN....


It contains one N for each day in which there have been no crashes and one Y for each day in which there has been at least one crash.

I calculate the probability $P(E)$ of a car crash happening on a certain day with $$P(E) = \frac{num(Y)}{num(Y)+num(N)} = 0.23$$

What is the correct way to calculate the probability $P(E_n|E_{n-1})$, i.e the probability that there will be crash on day $n$ if there has been a crash in the day $n-1$?

I know that the events E are not independent, so I cannot just calculate it as $0.23\times0.23=0.05$, and that value looks off anyway. This is a summary of the data, the full dataset is available at http://pastie.org/10370493.

\begin{array} {|r|r|} \hline days & 2924 \\ \hline N & 2250 \\ \hline Y & 674 \\ \hline Y\ days\ followed\ by\ another\ Y\ day & 101 \\ \hline Y\ days\ followed\ by\ two\ another\ Y\ day & 15 \\ \hline \end{array}

Now, my question is, and how can I calculate the expected probability of tomorrow having a crash if I know that in each of the last N days a crash have happened?

• I would use a numeric approach. Lets see if a well behaved random number generator gives anywhere near the same number of sequences that your results do. – EngrStudent Aug 24 '15 at 10:57
• The data is definitely not random: slicing it in pieces of 300 data points shows a consistent, non seasonal, 0.2 frequency of days with crashes. But the data is available if you want to try. – gioele Aug 24 '15 at 11:12

So I did a simulation, and my perfect binomial lines up reasonably well with your real-world data.

Code:

#parameters
N <- 2924
p <- 0.23

m <- 50000

#initializations
ysum <- as.numeric(matrix(nrow =m,ncol = 1))
set.seed(1)

#wrapper loop
for (k in 1:m){

#draw samples
y <- rbinom(size=1,prob=p,n=N)

#count pairs (I'm sure there are faster ways)
yy <- 0
for (i in 2:N){
if ( (y[i]==1) & (y[i-1]==1)){
yy <- yy+1
}
}

#store pair-count into variable "ysum"
ysum[k] <- yy
}

#make quantile plot
qqnorm(ysum)
grid()

#show summary of values
summary(ysum)


The graph of the distribution (qqnorm) was:

The text output of the summary was:

> summary(ysum)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
94.0   145.0   154.0   154.6   164.0   224.0


In 50,000 attempts to get the lowest value, I barely get where your data does. It is a far tail, but it is within the tail. Your results are out near 4 standard deviations from the mean.

What you are looking at (imo) falls into long-run/short-run statistics. I have seen a few items about "longest" but it is interesting that your data is meeting the "shortest" so well. There is physics driving that. Why is it not meeting the mean? A betting person would expect it to line up with the mean, not with an extreme minimum.

The level of minimum is on the order of 3 in 50,000 or 1 in 16,600. I had to run ~16k re-runs to get a minimum as low as yours. That is very atypical. One in a few hundred, or even a few thousand might not be as strange. This only occurs on the order of 1 in 10,000 sequences. And your data is from years of real-world data.

Thoughts:

• If a wreck is bad enough to block the intersection then there will be less traffic.
• If a wreck is bad enough to hit the news, then folks might avoid that route.
• What if I am a few decimal points off in the probability and it is 0.225 not 0.23? That might move the mean over some. This is why I like using at least 4 decimal places to work with - so my answers are accurate to 2 or 3 decimal places.

Possible references:

And so I answer the "3 days i a row" by changing the code from looking at pairs to looking at triples. It is a good thing that I made sure that my sample size was sufficient to re-create the 2-days in a row. I don't know if I am "in the ballpark" but I am at least in the same county.

The distribution was:

the summary was:

> summary(ysum)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
11.00   30.00   35.00   35.49   40.00   72.00


If I were a betting man, and if the 2-day sequence had NOT been in the tail, then I would guess about 35 times there was a 3-day sequence. Note that 4 in a row is counted as 3-days twice.

So lets think .. if there are 2 in a row, how many are needed to make it 3 in a row? Just one more positive outcome. As long as they are independent draws, then the likelihood is right about 23%.

PS:
So I looked at arima models, and the best fit is a first order moving average. model = (0,0,1). This suggests that the runs aren't exactly independent.

Part of the solution:

You can use the Poisson Distribution to model the number of accidents.

The rate of accidents is given by: $\lambda = \frac{|Y|}{|Y|+|N|} \text{accidents}/\text{day}$

Let $X_t$ be a random variable representing the number of accidents occurring upto and including the $t^{th}$ day, then, $P(X_t = k) = \frac{e^{-\lambda t}(\lambda t)^k}{k!}$

Then, the probability of observing $k+1$ car crashes the next day of haing observed $k$ car crashes is:

\begin{align} P(X_{t+1}=(k+1)|X_t=k)&=\frac{P(X_{t+1}=k+1,X_t=k)}{P(X_t=k)}\\ &=\frac{P(X_t=k|X_{t+1}=k+1)P(X_{t+1}=k+1)}{P(X_t=k)}\\&= \frac{P(X_1=1)P(X_{t+1}=k+1)}{P(X_t=k)}\\ &= \frac{e^{-\lambda}\lambda}{1!} \times e^{-\lambda}\lambda \frac{t+1}{k+1}(1+\frac{1}{t})^k \end{align}