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?
 A: 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:


*

*http://gato-docs.its.txstate.edu/mathworks/DistributionOfLongestRun.pdf

*https://math.stackexchange.com/questions/59738/probability-for-the-length-of-the-longest-run-in-n-bernoulli-trials
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.
A: 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}
$$  
