So I have an event that can occur on a daily basis with probabilty 1%. There is no dependence between days. Now I can calculate the probabality of the event occuring X times in Y days:

  COMBIN(Y,X) * (0.01 ^ X) * [0.99 ^ (Y - X)]

where COMBIN() is the number of combinations for X in Y. Now what I'd like to calculate is the probablity of X events in any Y window for the last Z days, e.g what is the probability of 5 events occuring in at least one 30 day window during the last 10 years. Any assistance much appreciated.


Note: Original question edited to add "at least" qualifier per MansT suggestion below.

  • $\begingroup$ If this is homework, it should have the homework tag. $\endgroup$ – Peter Flom Oct 1 '12 at 21:50
  • $\begingroup$ Not homework, real world problem cast in generalised terms... my stats is a bit rusty $\endgroup$ – Richard H Oct 3 '12 at 7:53
  • $\begingroup$ Do you mean "the probability of 5 events in at least one 30 day window for the last 10 years"? $\endgroup$ – MånsT Oct 10 '12 at 12:15
  • $\begingroup$ @MånsT - yes exactly $\endgroup$ – Richard H Oct 10 '12 at 12:23
  • 2
    $\begingroup$ Analytical answers are preferable, but given its not homework, I am curious why you don't use brute force Monte Carlo simulation. $\endgroup$ – phaneron Oct 12 '12 at 17:00

If the days are truly independent as you say then let us denote an event A as X events in any Y day period and denote P(A), the probability of X events in Y days:

P(A) = COMBIN(Y,X) * (0.01 ^ X) * [0.99 ^ (Y - X)]

This is true for any Y day window provided the probability of an event occurring is not time dependent.

Then you can apply the Binomial probability again to the event A. We have to be careful because now the event A is not independent for overlapping windows, i.e. for a window of length 10, the probability of event A on days 1-10 is not independent of days 2-11. Thus the following is assuming non-overlapping windows to ensure independence. Now, because you want the probability of at least one non-overlapping window having an event, it is easier to do:

1 - P(no overlapping windows satisfying event A).

Let n be the number of non-overlapping windows of size Y in time period Z. Then the probability of X events in a non-overlapping window of size Y over a time period Z is:

1 - COMBIN(n,0) * P(A)^0 * [1-P(A)]^(Y-0)

= 1 - n * [1-P(A)]^Y

  • $\begingroup$ This would answer the question if it were about months and not windows. I believe that a window is defined as 30 consecutive days (correct me if I'm wrong, @Richard). January 1-30 is a window, but so is January 2-31. The windows are not independent in that case. Indeed, window $i+1$ can only have $X$ events if window $i$ has at least $X-1$ events. $\endgroup$ – MånsT Oct 10 '12 at 16:24
  • $\begingroup$ MansT is correct, there is an implied assumption that the windows are independent. I've edited my response to make it clear that the solution is for non-overlapping windows. $\endgroup$ – adunaic Oct 12 '12 at 8:35
  • $\begingroup$ @aduniac - yes I am looking for a solution for over-lapping windows, not adjacent non-overlapping windows. $\endgroup$ – Richard H Oct 18 '12 at 8:20
  • $\begingroup$ @MansT - your understanding of the problem is correct: probability of X number of events in any 30 day window in the last 10 years. $\endgroup$ – Richard H Oct 18 '12 at 8:22

Here is some code, a brute-force approach, to get a ballpark for the probability.

In this I assume that the particular 30 day window does not matter. If there were 3 at the end of one month and 2 at the beginning of the next, I still count it as 5 in a row.


#number of runs
N <- 5e5

#max number of events per month (simulation tops out around 4)
n_max <- 10

#other parameters
n_years <- 10
days_per_year <- 365.24

#would be number of rows if year change was reset
n <- floor(n_years*days_per_year)

#sample is wrapped in this to handle workspace and memory
#find run lengths
my_rle <- rle(rbinom(n = n*N,prob = 0.01,size = 1))

#find run-lengths that refer to "TRUE" values sequences
idx <- which(my_rle$values ==1)

#subset out non-zero runs
cus_y <- my_rle$lengths[idx]

#pre-declare for loop
store<- numeric(length = n_max)

#put zeros into a single bin
store[1] <- sum(my_rle$lengths[-idx])/n/N

#find bin frequenceis
for (j in 1:(n_max-1)){
  store[j+1] <- length(which(cus_y==j,arr.ind=T))/n/N

#stage for plot and model
x    <- 0:(n_max-1)
y <- log10(store)

#subset to non-nan values
y1 <- y[1:5]
x1 <- seq(from=0,to=4,by=1)

#fit model
est <- lm(y1~x1)
summary( est)

x2 <- 5
y2 <- est$coefficients[1]+est$coefficients[2]*x2

#main plot
plot(x,y,ylab="Log10 frequency", xlab="run length",ylim = c(-10,0))


Here is the plot that I get

enter image description here

The fit gave me this summary:

> summary( est)

lm(formula = y1 ~ x1)

         1          2          3          4          5 
 0.0072500  0.0007404 -0.0025619 -0.0260976  0.0206690 

            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.01161    0.01528   -0.76    0.503    
x1          -1.99795    0.00624 -320.21 6.72e-08 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01973 on 3 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 1.025e+05 on 1 and 3 DF,  p-value: 6.717e-08

Using the coefficients gives -10.00137 as the expected log10 of frequency for runs of 5 in a row. This is ~1e-10. The estimated probability for a 5-element sequence at any time is 1e-10, or (0.01^5).


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.