**Robert Grant's stats blog » R**, and kindly contributed to R-bloggers)

I got a bit bored (sorry Mike), and wrote this. I didn’t take long (I tell you that not so much to cover my backside as to celebrate the majesty of R). First, I estimated probabilities of a day being rainy if the previous was dry, and of it being rainy if the previous day was rainy. I couldn’t be bothered with any thinking, so I used ABC, which basically is an incredibly simple and intuitive way of finding parameters than match data. I started with 16 rain days in July and 15 in August, in my home town of Winchester (which is also Swithun’s hood) from here, and that led to probabilities that I plugged into a simulator (it’s only AR1, but I’m no meteorologist). That ran 10,000 years of 40-day periods (there’s a sort of run-in of ten days to get to a stable distribution first; it’s basically a Markov chain), and not a single one had rain for 40 days.

It ain’t gon’ rain.

# Estmiate Winchester July/August rainy day transition probabilities # We'll use Approximate Bayesian Computation abciter<-1000 drytorain<-seq(from=0.15,to=0.35,by=0.01) raintorain<-seq(from=0.5,to=0.7,by=0.01) ldr<-length(drytorain) lrr<-length(raintorain) pb<-txtProgressBar(0,ldr*lrr*abciter,style=3) loopcount<-1 prox<-matrix(NA,nrow=ldr,ncol=lrr) for(i in 1:ldr) { for(j in 1:lrr) { trans<-c(drytorain[i],raintorain[j]) rainydays<-rep(NA,abciter) for(k in 1:abciter) { setTxtProgressBar(pb,loopcount) runin<-rep(NA,10) runin[1]<-rbinom(1,1,0.4) for (m in 2:10) { runin[m]<-rbinom(1,1,trans[runin[m-1]+1]) } days<-rep(NA,40) days[1]<-rbinom(1,1,trans[runin[10]+1]) for (m in 2:40) { days[m]<-rbinom(1,1,trans[days[m-1]+1]) } rainydays[k]<-sum(days) loopcount<-loopcount+1 } prox[i,j]<-sum(abs(rainydays-15.5)<1)/abciter rainydays<-rep(NA,abciter) } } close(pb) image(prox) # I'm going to go with P(rain | dry)=0.32, P(rain | rain)=0.51 # St Swithun's Day simulator prain<-c(0.32,0.51),nrow=2) iter<-10000 runs<-rep(NA,iter) pb<-txtProgressBar(0,iter,style=3) for (i in 1:iter) { setTxtProgressBar(pb,i) runin<-rep(NA,10) runin[1]<-rbinom(1,1,0.4) for (j in 2:10) { runin[j]<-rbinom(1,1,trans[runin[j-1]+1]) } days<-rep(NA,40) days[1]<-rbinom(1,1,trans[runin[10]+1]) for (j in 2:40) { days[j]<-rbinom(1,1,trans[days[j-1]+1]) } runs[i]<-max(rle(days)$lengths) } close(pb) print(paste("There were ",sum(runs==40), " instances of St Swithun's Day coming true, over ", iter," simulated years.",sep="")) hist(runs)</code>

**leave a comment**for the author, please follow the link and comment on their blog:

**Robert Grant's stats blog » R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...