abc.R

# @file: abc.R
# @author: Carl Boettiger <cboettig@gmail.com>
# @date: 11 April, 2011
# Description: A trivial example of Approximate Bayesian Computing
# We will simulate under a Gaussian process to determine 

# The target data:
TrueMean <- 7     ## irrelevant note: "True" is rather a frequentist name for this
TrueSD <- 1
Npts <- 100
X <- rnorm(Npts, TrueMean, TrueSD)

# Set summary stat function
S <- function(X_prime) c(mean(X_prime), sd(X_prime)) 

# set the target 
s <- c(ObservedMean <- mean(X), ObservedSD <- sd(X))

# This should be done m times, and we then do a regression of these pairs of means and sds: 
m <- 1000

sims <- sapply(1:m, function(i){
# Prior distributions on both parameters 
PriorMean <- abs(rnorm(1, mean=2, sd=10))
PriorVar <- abs(rnorm(1, mean=1, sd=100))

# simulate the process
X_prime <- rgamma(Npts, PriorMean, sqrt(PriorVar))
# Evaluate summary statistics



## Output is summary statistic and the parameter that achieved that value
c(summary=S(X_prime)[1], parameter=PriorMean)
})

#Tolerance
delta <- 50

# Select only those points near the estimated value of the statistic, 
# where we can expect a linear relationship between the statistic and the 
# parameter
good <- sims[ , abs(sims[1,]-S(X)[1]) <  delta]

par(mfrow=c(1,2)) #2x1 plotting window
plot(t(sims))

# Compute our regression
fit <- lm(good[1,]~good[2,])
b <- fit$coeff[2]

posterior = good[2,] - b*(good[1,] - S(X)[1])

# The posterior distribution for the parameter is just the residuals 
plot(density( posterior ))