St Swithun’s Day simulator

July 15, 2015
By

(This article was first published on 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>

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



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers


Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)