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:


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?

  • $\begingroup$ 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. $\endgroup$ Aug 24 '15 at 10:57
  • $\begingroup$ 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. $\endgroup$
    – 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.


N <- 2924
p <- 0.23

m <- 50000

ysum <- as.numeric(matrix(nrow =m,ncol = 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

#show summary of values

The graph of the distribution (qqnorm) was:

enter image description here

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.


  • 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:

enter image description here

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%.

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} $$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.