abc.jl

using Distributions

#=

ABC rejection smapling

This code estimates the value of two parameters (optimum and tolerance)
for simulated data of fitness along a gradient of x values. The fitness is
given by a Gaussian function.

Parameters:

 n_samples - number of samples to generate
 n_tests   - number of tests to do
 epsilon   - value above which simulated parameters are rejected

Timothée Poisot <[email protected]> July 23, 2014

=#


const n_samples = 80
const x = rand(Uniform(-1, 1), n_samples)
const dev = rand(Uniform(-0.05, 0.05), length(x))
const n_tests = 10000
const epsilon = 2.0

fit = (x, x0, xi) -> exp(-(x-x0)*(x-x0)/(2*xi))

fits = zeros(n_samples)
@inbounds for i in [1:n_samples]
   fits[i] = fit(x[i], -0.2, 0.22) + dev[i]
end

x0 = rand(Uniform(-2, 2), n_tests)
xi = rand(Uniform(0.0, 1.0), n_tests)

x0_kept = zeros(0)
xi_kept = zeros(0)

@inbounds for test in [1:n_tests]
   t_fits = zeros(n_samples)
   for i in [1:n_samples]
      t_fits[i] = fit(x[i], x0[test], xi[test])
   end
   sum_of_squares = sum((fits-t_fits).*(fits-t_fits))
   if sum_of_squares <= epsilon
      append!(x0_kept, [x0[test]])
      append!(xi_kept, [xi[test]])
   end
end


@printf "x0 ~ {%4f %4f} - ξ ~ {%4f %4f}\n" mean(x0_kept) var(x0_kept) mean(xi_kept) var(xi_kept)