#ST407 Monte Carlo Methods 2012/2013
#Practical 7: From Data Augmentation to ABC
#Question 1
#Part b
rlazy<-function(n, y){
#A very inefficient way of sampling from the posterior of theta
#Via data augmentation/ABC
#n = number of trials in the binomia, set up
#y= the actual observed number of successes in our Binomial experiment
#theta in this case is the probability of success
yp<--1 #Initial setup to ensure we enter the while loop at least once
while(yp!=y){
theta<-rbeta(n=1, shape1=2, shape2=2) #Sample one proposed theta from beta prior (the n rbeta argument is NOT the same n used at the beginning)
yp<-rbinom(n=1, size=n, prob=theta) #Sample a simulate Binomial outcome (synthetic data), based on this sample value of theta
#If yp!=y, we discard all proposed samples, and we repeat everything in this loop again
}
theta #Return the value of theta where yp=y
}
system.time(rlazy(n=10, y=10))
system.time(rlazy(n=1000, y=10))
system.time(rlazy(n=10000, y=10))
system.time(rlazy(n=100000, y=10))
#System times are dependent on the specific computer you are running it on.
#What you should observe is that as the number of trials increases, the required time to perform the function increases.
#This is because the probability of 10 successes occurring in large n trials slowly becomes smaller.
#And also that the associated small theta required, is not as probable under the prior chosen
#There is thus less chance that we'll select and simulate a theta that will obtain such an event
#A slightly more sophisticated (and generally advisable) approach is to repeat the function such that we get more than one theta, say 10, to capture the uncertainty regarding theta.
#This also reduces Monte Carlo variability a little.
#We could similarly average the system times to produce a single theta for each scenario.
rlazyrep<-function(n, y, r){
#If we want to return r samples under the rabc function
#Adam's way of doing things, using a for loop
#reps<-vector("numeric", r) #Storage
#for(i in 1:r){
# reps[i]<-rlazy(n=10, y=10)
#}
#Christopher's way, avoiding the use of a for loop
reps<-replicate(r, rlazy(n=n, y=y)) #replicate just performs the same function however many times you need it to
reps
}
#If we want to obtain 100 posterior values of theta, r=100
theta.r<-100
system.time(rlazyrep(n=10, y=10, r=theta.r))/theta.r
system.time(rlazyrep(n=1000, y=10, r=theta.r))/theta.r
system.time(rlazyrep(n=10000, y=10, r=theta.r))/theta.r
system.time(rlazyrep(n=100000, y=10, r=theta.r))/theta.r
#Question 2, a more efficient approach
rabc<-function(n,y, k=0.5){
#A slightly more efficient approach in obtaining a posterior for theta
#n and y are the same as before (number of trails, and number of successes)
#k = the required closeness between the synthetic and real data to accept the proposed theta
yp<- -2*k #Again, a initilisation procedure to ensure the while loop runs
while(abs(yp-y)>k){ #Run the while loop until the condition in step 3 is not satisfied
theta<-rbeta(n=1, shape1=2, shape2=2) #Sample one proposed theta from beta prior (the n rbeta argument is NOT the same n used at the beginning)
yp<-rbinom(n=1, size=n, prob=theta) #Sample a simulate Binomial outcome (synthetic data), based on this sample value of theta
}
theta #Accept and return the theta if the synthtetic and real data are close enough
}
rabcrep<-function(n, y, k, r){
#If we want to return r samples under the rabc function
#Adam's way of doing things, using a for loop
#reps<-vector("numeric", r) #Storage
#for(i in 1:r){
# reps[i]<-rabc(n,y,k)
#}
#Christopher's way, avoiding the use of a for loop
reps<-replicate(r, rabc(n=n,y=y,k=k)) #replicate just performs the same function however many times you need it to
reps
}
#If we want to obtain 100 posterior values of theta, r=100
theta.r<-100
theta.k<-1
system.time(rabcrep(n=10, y=10, k=theta.k, r=theta.r))/theta.r
system.time(rabcrep(n=1000, y=10, k=theta.k, r=theta.r))/theta.r
system.time(rabcrep(n=10000, y=10, k=theta.k, r=theta.r))/theta.r
system.time(rabcrep(n=100000, y=10, k=theta.k, r=theta.r))/theta.r
#You already notice that for the same value of k=1, the timiing are much quicker
#This is because we are being less restrictive in that we need to exactly obtain the real data and are throwing away less proposals
#Now with different values of k
theta.k<-2
system.time(rabcrep(n=10, y=10, k=theta.k, r=theta.r))/theta.r
system.time(rabcrep(n=1000, y=10, k=theta.k, r=theta.r))/theta.r
system.time(rabcrep(n=10000, y=10, k=theta.k, r=theta.r))/theta.r
theta.k<-3
system.time(rabcrep(n=10, y=10, k=theta.k, r=theta.r))/theta.r
system.time(rabcrep(n=1000, y=10, k=theta.k, r=theta.r))/theta.r
system.time(rabcrep(n=10000, y=10, k=theta.k, r=theta.r))/theta.r
theta.k<-4
system.time(rabcrep(n=10, y=10, k=theta.k, r=theta.r))/theta.r
system.time(rabcrep(n=1000, y=10, k=theta.k, r=theta.r))/theta.r
system.time(rabcrep(n=10000, y=10, k=theta.k, r=theta.r))/theta.r
theta.k<-5
system.time(rabcrep(n=10, y=10, k=theta.k, r=theta.r))/theta.r
system.time(rabcrep(n=1000, y=10, k=theta.k, r=theta.r))/theta.r
system.time(rabcrep(n=10000, y=10, k=theta.k, r=theta.r))/theta.r
#What you should observe for different k is that the required computational time decreases quite dramatically.
#This is not surprising as we are being less restrictive with respect to our acceptance proecure.
#However, the values of theta returned are not as accurate from the relaxed condition.
#part c
theta.n<-100; theta.r<-5000
abc.k1<-rabcrep(theta.n, y=10, k=1, r=theta.r)
abc.k2<-rabcrep(theta.n, y=10, k=2, r=theta.r)
abc.k3<-rabcrep(theta.n, y=10, k=3, r=theta.r)
abc.k4<-rabcrep(theta.n, y=10, k=4, r=theta.r)
abc.k5<-rabcrep(theta.n, y=10, k=5, r=theta.r)
x.vec<-seq(0,1, len=1000)
density.vec<-dbeta(x.vec, 10+2, theta.n+2-10)
hist(abc.k1, probability=TRUE); points(x.vec, dbeta(x.vec, 10+2, theta.n+2-10), col="red")
hist(abc.k2, probability=TRUE); points(x.vec, dbeta(x.vec, 10+2, theta.n+2-10), col="red")
hist(abc.k3, probability=TRUE); points(x.vec, dbeta(x.vec, 10+2, theta.n+2-10), col="red")
hist(abc.k4, probability=TRUE); points(x.vec, dbeta(x.vec, 10+2, theta.n+2-10), col="red")
hist(abc.k5, probability=TRUE); points(x.vec, dbeta(x.vec, 10+2, theta.n+2-10), col="red")
#The prior and the likelihood are conjugate so that the posterior has the form of a Beta distribution.
#You should notice that as k increases, the histogram of the posterior samples of theta become worse and worse.
#You should see a reduction in the mass in concentration int eh peak and will lead to a slight reduction in strength of any conclusions (more diffuse and less efficient)
#Thus a trade-off between speed and approximation as k increases.