Your first step should always be to plot. Here is a plot of when your occurrences happened, a histogram of the waiting times, and a time series plot of the waiting times:
The first thing that jumps out at us is that the last waiting time was enormously long, corresponding to the single spike at the far right of the histogram - this last waiting time was 660 seconds, all the others no more than 137 seconds.
The second observation (see the right-hand panel) is that waiting times apparently were already increasing before the very last occurrence. The variance of the waiting times was also increasing.
It seems to me like the first thing you should do is to investigate just what caused these two effects. Understanding this should inform your subsequent analysis and prediction. Did something change over time, and if so, is the change persistent, or will matters revert to the previous state, or something else?
You could then try modeling the waiting times with an appropriate distribution. The exponential distribution is commonly used for that, but since your data are discrete, you might want to model them as negative binomial (which describes the waiting times in a Bernoulli process). You can include predictors using negative binomial regression, see the textbook Negative Binomial Regression by Hilbe.
Alternatively, you could use standard time series forecasting techniques to forecast the likely next waiting time. For instance, you could fit an Exponential Smoothing model to the sequence of waiting times and forecast that out, which will give you a forecast of 188 seconds for the expected next waiting time, along with prediction intervals. You can even simulate from this series to get a probability density for when the next occurrence will happen:
However, do note that understanding your data is definitely more important than finding the most sophisticated model!
R code:
occurrences <- c(25,26,34,56,72,76,76,86,88,88,106,148,151,151,161,195,200,
214,215,215,220,231,243,245,247,257,263,265,288,295,314,339,342,342,353,
368,407,413,436,447,469,470,472,505,513,557,566,598,609,623,663,
676,683,687,776,789,850,875,921,1058,1078,1167,1255,1292,1952)
waiting_times <- diff(occurrences)
par(mfrow=c(1,3),las=1)
plot(occurrences,rep(1,length(occurrences)),type="h",lwd=2,
yaxt="n",ylab="",xlab="Second",main="Occurrences",ylim=c(0,1.3))
hist(waiting_times,breaks=seq(-0.5,max(waiting_times)+0.5))
plot(waiting_times,type="l")
library(forecast)
model <- ets(waiting_times)
forecast(model,h=1)
set.seed(1)
sims <- replicate(10000,simulate(model,nsim=1,bootstrap=TRUE))
table_sims <- hist(sims,breaks=seq(floor(min(sims)),ceiling(max(sims))),plot=FALSE)
plot(occurrences,rep(1,length(occurrences)),type="h",lwd=2,
yaxt="n",ylab="",xlab="Second",main="Occurrences with forecast",
ylim=c(0,1.3),xlim=c(0,max(occurrences)+max(sims)))
points(max(occurrences)+table_sims$mids,table_sims$density,
type="h",col="grey",lwd=3)
