St Swithun’s Day simulator

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)
Advertisements

3 Comments

Filed under R

3 responses to “St Swithun’s Day simulator

  1. 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…

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s