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
trans<-c(0.32,0.51)
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)

Part freelance data science guy, part Lattepapa. I'm interested in Bayesian modelling, data visualization and ever-better ways to teach people about stats. I used to be a university lecturer and before that a faceless bean-counter. "I don't have much, but what I have is gold."
View all posts by Robert

Yeah, it’s wordpress.com being a pain. If you edit a post it switches between the ascii symbols and the html in a rather unpredictable way. You’re right, those are mostly assignments and quotes. I’ll fix it, but you get the picture…

Assuming ” < ” is ” <- ". Fix rendering of code…

Yeah, it’s wordpress.com being a pain. If you edit a post it switches between the ascii symbols and the html in a rather unpredictable way. You’re right, those are mostly assignments and quotes. I’ll fix it, but you get the picture…

Hmmm, this looks like it could be useful. I’m too scared to edit this one now though. http://wordpress.stackexchange.com/questions/46482/how-to-stop-wordpress-mangling-r-syntax