Introduction to ABC

Here we will perfom some basic Approximate Bayesian Computation (ABC) inference of the mean and the standard deviation of a Normal distribution. We will use a rejection sampling algorithm, and then we will use a simple MCMC algorithm.

First, let’s generate our data.

We are assuming that our data consists of 100 samples from a Normal distribution of mean 4.3 and sd 2.7.

number_of_data_points = 100
data =  rnorm(number_of_data_points, mean = 4.3, sd = 2.7)

summary(data)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -2.540   2.798   4.355   4.440   5.999  13.150 

Rejection sampling algorithm

In this algorithm, we will simulate a large number (n_iterations) of samples from Normal distributions of parameters mu and sigma themselves sampled from uniform distributions, and retain only those samples that “look like” the observed data. In practice, we will try three different ways of quantifying how much simulated samples look like the observed data. Two measures will be based on quantiles, and one measure will be based on the moments of the distribution (the mean and the variance).

The priors on mu and sigma, the parameters of the Normal distribution

We assume uniform priors from which we draw mu and sigma.

draw_mu <- function () {
  return (runif(1, min=0, max=10))
}
draw_sigma <- function () {
  return (runif(1, min=0, max=10))
}

The function to simulate a data set from the parameters mu and sigma

simulate_data <- function (number_of_data_points, mu, sigma) { 
  return(rnorm(number_of_data_points, mean = mu, sd = sigma))
}

A function to compare simulated and observed samples

# Function to compute the quantiles.
# We choose to use 3 quantiles.
compute_quantiles <- function(data) {
  return (quantile(data, probs=c(0.1, 0.5, 0.9)))
}
# First method to compare a simulated sample to the observed data
compare_quantiles_with_squared_distance <- function (true, simulated) {
  distance = sqrt(sum(mapply(function(x,y) (x-y)^2, true, simulated)))
  return(distance)
}

A function to accept or reject a sample based on its distance(s) to the observed data

# Accept or reject based on the first method to compare a simulated sample to the observed data
accept_or_reject_with_squared_distance <- function (true, simulated, acceptance_threshold) {
  distance = compare_quantiles_with_squared_distance(compute_quantiles(true), compute_quantiles(simulated))
  if((distance < acceptance_threshold) ) return(T) else return(F)
}

Combining all the above elements to get the full rejection sampler

We implement the function so that it can take any accept_or_reject function. This function will create a data frame with three columns: the first, named “accepted_or_rejected”, tells whether a given sample has been accepted or rejected by the algorithm. The second and third, “sampled_mus” and “sampled_sigmas”, contain the sampled values of the parameters.

sample_by_rejection <- function (true_data, n_iterations, acceptance_threshold, accept_or_reject_function) {
  number_of_data_points = length(true_data)
  accepted_or_rejected <- vector(length = n_iterations)
  sampled_mus <- vector(length = n_iterations, mode = "numeric")
  sampled_sigmas <- vector (length = n_iterations, mode = "numeric")
  for (i in 1:n_iterations){
    mu <- draw_mu()
    sigma <- draw_sigma()
    parameters = list("mu"=mu, "sigma"=sigma)
    simulated_data <- simulate_data(number_of_data_points, mu, sigma)
    accepted_or_rejected[i] = accept_or_reject_function(true_data, simulated_data, acceptance_threshold)
    sampled_mus[i] = mu
    sampled_sigmas[i] = sigma
  }
  return(data.frame(cbind("accepted_or_rejected" = accepted_or_rejected, "sampled_mus" = sampled_mus, "sampled_sigmas" = sampled_sigmas)))
}

Now we perform inference using our distance functions, and analyze the results

First, the squared distance

We run the analysis:

sampled_parameter_values_squared_distances = sample_by_rejection(data, 200000, 0.5, accept_or_reject_with_squared_distance)

Analysis of the output

How many samples have been accepted out of 200000?

sum(sampled_parameter_values_squared_distances$accepted_or_rejected)

[1] 332

Using Coda to summarize the samples and do some plots

# Useful library used here for plotting mostly.
library(coda)
rej_samples_squared_distances_as_mcmc = mcmc(sampled_parameter_values_squared_distances[which(sampled_parameter_values_squared_distances$accepted_or_rejected==1),c(2,3)])

summary(rej_samples_squared_distances_as_mcmc)

Iterations = 1:332
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 332 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                Mean     SD Naive SE Time-series SE
sampled_mus    4.249 0.3316  0.01820        0.01721
sampled_sigmas 2.834 0.2855  0.01567        0.01567

2. Quantiles for each variable:

                2.5%   25%   50%   75% 97.5%
sampled_mus    3.623 4.000 4.270 4.474 4.848
sampled_sigmas 2.339 2.629 2.802 2.998 3.442

plot(rej_samples_squared_distances_as_mcmc)

Is there some correlation between subsequent iterations?

autocorr.plot(rej_samples_mean_and_variance_as_mcmc)

Of course not! The results would be the same with the other rejection samplers.

Trying two other distance functions

# Second method to compare a simulated sample to the observed data
compare_quantiles_with_median_and_spread <- function (true_quantiles, simulated_quantiles) {
  distances=vector(length=2, mode="numeric")
  distances[1] <- abs(true_quantiles[2]-simulated_quantiles[2])
  true_spread <- true_quantiles[3]-true_quantiles[1]
  simulated_spread <- simulated_quantiles[3]-simulated_quantiles[1]
  distances[2] <- abs(true_spread - simulated_spread)
  return(distances)
}
# Third method to compare a simulated sample to the observed data
compare_distributions_with_mean_and_variance <-function (true, simulated) {
  # comparison with the observed summary statistics
  diffmean <- abs(mean(true) - mean(simulated))
  diffsd <- abs(sd(true) - sd(simulated))
  return(c(diffmean, diffsd))
}

Two other accept or reject functions based on these distances

# Accept or reject based on the second method to compare a simulated sample to the observed data
accept_or_reject_with_median_and_spread <- function (true, simulated, acceptance_threshold) {
  distances = compare_quantiles_with_median_and_spread(compute_quantiles(true), compute_quantiles(simulated))
  if((distances[1] < acceptance_threshold) & (distances[2] < 4*acceptance_threshold) ) return(T) else return(F)
}
# Accept or reject based on the third method to compare a simulated sample to the observed data
accept_or_reject_with_mean_variance <- function (true, simulated, acceptance_threshold) {
  differences = compare_distributions_with_mean_and_variance(true, simulated)
  if((differences[1] < acceptance_threshold) & (differences[2] < 2*acceptance_threshold)) return(T) else return(F)
}

We run the analyses with these distance functions

Second, the median and spread

We run the analysis:

sampled_parameter_values_median_and_spread = sample_by_rejection(data, 200000, 0.1, accept_or_reject_with_median_and_spread)

Third, the mean and variance

We run the analysis:

sampled_parameter_values_mean_and_variance = sample_by_rejection(data, 200000, 0.1, accept_or_reject_with_mean_variance)

How many samples have been accepted?

sum(sampled_parameter_values_median_and_spread$accepted_or_rejected)

[1] 122

sum(sampled_parameter_values_mean_and_variance$accepted_or_rejected)

[1] 139

Summarizing the posterior distributions

rej_samples_median_and_spread_as_mcmc = mcmc(sampled_parameter_values_median_and_spread[which(sampled_parameter_values_median_and_spread$accepted_or_rejected==1),c(2,3)])
rej_samples_mean_and_variance_as_mcmc = mcmc(sampled_parameter_values_mean_and_variance[which(sampled_parameter_values_mean_and_variance$accepted_or_rejected==1),c(2,3)])

summary(rej_samples_median_and_spread_as_mcmc)

Iterations = 1:122
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 122 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                Mean     SD Naive SE Time-series SE
sampled_mus    4.359 0.3732  0.03378        0.03751
sampled_sigmas 2.824 0.3063  0.02773        0.03222

2. Quantiles for each variable:

                2.5%   25%   50%   75% 97.5%
sampled_mus    3.629 4.089 4.374 4.621 5.120
sampled_sigmas 2.279 2.622 2.791 3.063 3.381

summary(rej_samples_mean_and_variance_as_mcmc)

Iterations = 1:139
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 139 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                Mean     SD Naive SE Time-series SE
sampled_mus    4.442 0.2953  0.02505        0.02505
sampled_sigmas 2.840 0.2128  0.01805        0.01805

2. Quantiles for each variable:

                2.5%   25%   50%   75% 97.5%
sampled_mus    3.802 4.250 4.438 4.633 4.983
sampled_sigmas 2.500 2.679 2.829 2.953 3.313

Plotting the posterior distributions

plot(rej_samples_median_and_spread_as_mcmc)

plot(rej_samples_mean_and_variance_as_mcmc)

MCMC ABC from Florian Hartig (https://www.r-bloggers.com/a-simple-approximate-bayesian-computation-mcmc-abc-mcmc-in-r/)

 
 
# we want to use ABC to infer the parameters that were used. 
# we sample from the same model and use mean and variance
# as summary statistics. We return true for ABC acceptance when
# the difference to the data is smaller than a certain threshold
 
meandata <- mean(data)
standarddeviationdata <- sd(data)
 
ABC_acceptance <- function(par){
   
  # prior to avoid negative standard deviation
  if (par[2] <= 0) return(F) 
   
  # stochastic model generates a sample for given par
  samples <- rnorm(10, mean =par[1], sd = par[2])
 
  # comparison with the observed summary statistics
  diffmean <- abs(mean(samples) - meandata)
  diffsd <- abs(sd(samples) - standarddeviationdata)
  if((diffmean < 0.1) & (diffsd < 0.2)) return(T) else return(F)
}
 
 
# we plug this in a standard MCMC, 
# with the metropolis acceptance exchanged for the ABC acceptance
 
run_MCMC_ABC <- function(startvalue, iterations, thinning_interval=1, size_of_move=0.7){
 
    chain = array(dim = c(iterations+1,2))
    chain[1,] = startvalue
 
    for (i in 1:iterations){
         
        # proposalfunction
        proposal = rnorm(2,mean = chain[i,], sd= c(size_of_move,size_of_move))
         
        if(ABC_acceptance(proposal)){
            chain[i+1,] = proposal
        }else{
            chain[i+1,] = chain[i,]
        }
    }
    return(mcmc(chain, thin=thinning_interval))
}
 
posterior <- run_MCMC_ABC(c(4,2.3),300000, 1, 0.7)
plot(posterior)

summary(posterior)

Iterations = 1:300001
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 300001 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

      Mean     SD Naive SE Time-series SE
[1,] 4.553 1.5178 0.002771        0.17851
[2,] 3.361 0.7498 0.001369        0.04763

2. Quantiles for each variable:

      2.5%   25%   50%   75% 97.5%
var1 1.395 3.652 4.428 5.283 7.972
var2 2.056 2.797 3.341 3.842 4.851

Is there some correlation between subsequent iterations?

autocorr.plot(posterior)

Yes, there is a huge amount of correlation!

What is the effective sample size of the sample obtained with this MCMC?

effectiveSize(posterior)

     var1      var2 
 72.29052 247.79597 

which is obviously much less than the 300 000 iterations that have been run…

Let’s see how thinning affects the results:

posterior <- run_MCMC_ABC(c(4,2.3),300000, 500, 0.7)
plot(posterior)

summary(posterior)

Iterations = 1:1.5e+08
Thinning interval = 500 
Number of chains = 1 
Sample size per chain = 300001 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

      Mean    SD Naive SE Time-series SE
[1,] 4.258 1.015 0.001852         0.0780
[2,] 3.552 1.258 0.002297         0.1329

2. Quantiles for each variable:

      2.5%   25%   50%   75% 97.5%
var1 2.358 3.497 4.258 4.929 6.340
var2 2.008 2.635 3.205 4.001 6.916

autocorr.plot(posterior)

effectiveSize(posterior)

     var1      var2 
169.21251  89.59986 

Let’s see how changing the move width affects the results, first making it smaller:

posterior <- run_MCMC_ABC(c(4,2.3),300000, 1, 0.5)
plot(posterior)

summary(posterior)

Iterations = 1:300001
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 300001 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

      Mean    SD Naive SE Time-series SE
[1,] 4.308 1.054 0.001924        0.10503
[2,] 3.390 0.905 0.001652        0.07798

2. Quantiles for each variable:

      2.5%   25%   50%   75% 97.5%
var1 1.903 3.630 4.384 5.072 6.059
var2 1.976 2.702 3.255 3.965 5.304

autocorr.plot(posterior)

effectiveSize(posterior)

    var1     var2 
100.6622 134.6898 

Let’s see how changing the move width affects the results, now making it larger:

posterior <- run_MCMC_ABC(c(4,2.3),300000, 1, 0.9)
plot(posterior)

summary(posterior)

Iterations = 1:300001
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 300001 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

      Mean    SD Naive SE Time-series SE
[1,] 4.413 1.062 0.001939        0.07247
[2,] 3.394 1.016 0.001855        0.07581

2. Quantiles for each variable:

      2.5%   25%   50%   75% 97.5%
var1 2.299 3.727 4.380 5.059 6.791
var2 2.018 2.647 3.161 3.881 6.001

autocorr.plot(posterior)

effectiveSize(posterior)

    var1     var2 
214.7822 179.6936