rm(list=ls(all=TRUE)) # clear all variables ###################################################################### ########Problem 1: Binomial example ###################################################################### set.seed(123) #Here are some parameters you'll paste them all in to initialize the ABC but then you'll change them one at at time in parts a,b,c,d: n=20 #number of observations N=100 #generated sample size true.p=.75 #true parameter value data=rbinom(n,1,true.p) epsilon=0 #Tolerance #Here is the distance function you'll paste in to initialize the ABC but then you'll change it in part e: rho=function(y,x) abs(sum(y)-sum(x))/n #Keep this section in the clip board and paste it in the R console to re-run as needed while you change parameters: alpha.hyper = 1 #Hyper parameter beta.hyper = 1 #Hyper parameter p=numeric(N) #"Good" proposals stored in p for(i in 1:N){ d=epsilon+1 while(d>epsilon) { proposed.p=rbeta(1,alpha.hyper,beta.hyper) x=rbinom(n,1,proposed.p) d=rho(data,x) } p[i]= proposed.p } ######## Plots of the ABC - posterior and true posterior grid0 = seq(0, 1, length.out = 1000) dev.new(width = 7.5, height = 6) par(mar = c(4,4,2,2), oma = c(2,2,2,2)) hist(p, main = paste("Tolerance: ", epsilon, ", N:", N, ", n:",n, sep = ""), nclass = 25, freq = FALSE, xlim = c(0,1)) den1 <- density(p) lines(den1\$x,den1\$y,col = "black",lwd = 4) abline(v = true.p, col = "red", lty = 2, lwd = 4) lines(grid0,dbeta(grid0,sum(data)+1,n-sum(data)+1),col = "blue",lwd = 4) legend("topleft", c("True Posterior","Estimated Posterior",paste("True p: ", true.p, sep = "")), col = c("blue","black","red"), lty = c(1,1,2),lwd = 3) #Stop copying here for Problem 1 ###################################################################### ########Problem 2: Gaussian example ###################################################################### set.seed(123) #Here are some parameters you'll paste them all in to initialize the ABC but then you'll change them one at at time in parts a,b: n=25 #number of observations N=100 #particle sample size true.mu = 0 #true parameter value (considered unknown) sigma = 1 #true parameter value (considered known) #Here is the distance function you'll paste in to initialize the ABC but then you'll change it in part c: #******WARNING: when you change the distance function, it is a good idea to start by running it with a larger epsilon rho=function(y,x) abs(sum(y)-sum(x))/n epsilon=0.05 #Tolerance #Keep this section in the clip board and paste it in the R console to re-run as needed while you change parameters: mu.hyper = 0 #Hyper parameter sigma.hyper = 10 #Hyper parameter data=rnorm(n,true.mu,sigma) mu=numeric(N) #"Good" proposals stored in mu for(i in 1:N){ d= epsilon +1 while(d>epsilon) { proposed.mu=rnorm(1,0,sigma.hyper) #<--prior draw x=rnorm(n, proposed.mu, sigma) d=rho(data,x) } mu[i]= proposed.mu } ######## Plots of the ABC - posterior and true posterior grid0 = seq(-1, 1, length.out = 1000) mu1 = (mu.hyper/sigma.hyper^2 + sum(data)/sigma^2)/(1/sigma.hyper^2 + n/sigma^2) sigma1 = sqrt(1/(1/sigma.hyper^2 + n/sigma^2)) dev.new(width = 7.5, height = 6) par(mar = c(4,4,2,2), oma = c(2,2,2,2)) hist(mu, main = paste("Tolerance: ", epsilon, ", N:", N, ", n:",n, sep = ""), nclass = 25, freq = FALSE, xlim = c(-1,1)) den1 <- density(mu) lines(den1\$x,den1\$y,col = "black",lwd = 4) abline(v = true.mu, col = "red", lty = 2, lwd = 4) lines(grid0,dnorm(grid0,mu1,sigma1),col = "blue",lwd = 4) legend("topright", c("True Posterior","Estimated Posterior","True mu"), col = c("blue","black","red"), lty = c(1,1,2),lwd = 3) #Stop copying here for Problem 2 ###################################################################### ########Problem 3: Gaussian example - sequential ###################################################################### set.seed(123) #Here you can change the number of time steps in the sequential algorithm. Paste it as is to initialize the ABC but then you'll change it part b: time.steps = 10 #Here is the distance function you'll paste in to initialize the ABC but then you'll change it in part e: rho=function(y,x) abs(sum(y)-sum(x))/n #Keep this section in the clip board and paste it in the R console to re-run as needed while you change parameters: n=25 #number of observations N=100 #particle sample size true.mu = 0 #true parameter value (considered unknown) sigma = 1 #true parameter value (considered known) mu.hyper = 0 #Hyper parameter sigma.hyper = 10 #Hyper parameter data=rnorm(n,true.mu,sigma) epsilon = 1 #Starting tolerance weights = matrix(1/N,time.steps,N) mu=matrix(NA,time.steps,N) d=matrix(NA,time.steps,N) for(t in 1:time.steps){ print(t) if(t==1){ for(i in 1:N){ d[t,i]= epsilon +1 while(d[t,i]>epsilon) { proposed.mu=rnorm(1,0,sigma.hyper) #<--prior draw x=rnorm(n, proposed.mu, sigma) d[t,i]=rho(data,x) } mu[t,i]= proposed.mu }} else{ epsilon = c(epsilon,quantile(d[t-1,],.75)) mean.prev <- sum(mu[t-1,]*weights[t-1,]) var.prev <- sum((mu[t-1,] - mean.prev)^2*weights[t-1,]) for(i in 1:N){ d[t,i]= epsilon[t]+1 while(d[t,i]>epsilon[t]) { sample.particle <- sample(N, 1, prob = weights[t-1,]) proposed.mu0 <- mu[t-1, sample.particle] proposed.mu <- rnorm(1, proposed.mu0, sqrt(2*var.prev)) x <- matrix(rnorm(n,proposed.mu, sigma),n,1) d[t,i]=rho(data,x) } mu[t,i]= proposed.mu mu.weights.denominator <- sum(weights[t-1,]*dnorm(proposed.mu, mu[t-1,],sqrt(2*var.prev))) mu.weights.numerator <- dnorm(proposed.mu,0,sigma.hyper) weights[t,i] <- mu.weights.numerator/mu.weights.denominator }} weights[t,] <- weights[t,]/sum(weights[t,]) } ######## Plots of the ABC - posterior and true posterior library(fields) ###<----------------------------Install the library `fields,' and then load it here (for the tim.colors) col1 = tim.colors(time.steps) grid0 = seq(-1, 1, length.out = 1000) mu1 = (mu.hyper/sigma.hyper^2 + sum(data)/sigma^2)/(1/sigma.hyper^2 + n/sigma^2) sigma1 = sqrt(1/(1/sigma.hyper^2 + n/sigma^2)) time.step <-1 #<----------------------- Change the time step here den1 <- density(mu[time.step,],weights = weights[time.step,]) plot(den1\$x,den1\$y,"l", col = col1[time.step],lwd = 3, xlab = "mu", ylab = "Density", ylim = c(0,2.5)) abline(v = true.mu, col = "red", lty = 2, lwd = 3) den1 <- density(mu[time.step,],weights = weights[time.step,]) lines(den1\$x,den1\$y,"l", col = "magenta",lwd = 8) lines(grid0,dnorm(grid0,mu1,sigma1),col = "black",lwd = 5) legend("topright", c("True Posterior","ABC Posterior", "True mu"), col = c("black","magenta","red"), lty = c(1,1,2),lwd = 3) ######## To plot several (or all) the time steps ABC - posteriors together, run the code below: library(fields) ###<----------------------------Install the library `fields,' and then load it here (for the tim.colors) col1 = tim.colors(time.steps) grid0 = seq(-1, 1, length.out = 1000) mu1 = (mu.hyper/sigma.hyper^2 + sum(data)/sigma^2)/(1/sigma.hyper^2 + n/sigma^2) sigma1 = sqrt(1/(1/sigma.hyper^2 + n/sigma^2)) dev.new(width = 7.5, height = 6) par(mfrow = c(1,1), mar = c(4,4,2,2), oma = c(2,2,2,2)) den1 <- density(mu[1,],weights = weights[1,]) plot(den1\$x,den1\$y,"l", col = col1[1],lwd = 3, xlab = "mu", ylab = "Density", ylim = c(0,2.25)) for(which.t in 2:time.steps){ den1 <- density(mu[which.t,],weights = weights[which.t,]) lines(den1\$x,den1\$y,col = col1[which.t],lwd = 3) } abline(v = true.mu, col = "red", lty = 2, lwd = 3) den1 <- density(mu[time.steps,],weights = weights[time.steps,]) lines(den1\$x,den1\$y,"l", col = "magenta",lwd = 8) lines(grid0,dnorm(grid0,mu1,sigma1),col = "black",lwd = 5) legend("topright", c("True Posterior","ABC Posterior", "True mu"), col = c("black","magenta","red"), lty = c(1,1,2),lwd = 3) ######## Plots the decreasing tolerances plot(epsilon, pch = 19, lwd = 3, xlab = "Time step", ylab = "Tolerance") #Stop copying here for Problem 3