Model Selection – How to Calculate Weights in ABC-SMC for Bayesian Computation?

approximate-bayesian-computationmodel selectionmonte carlosequential-monte-carlo

I have been reading through the Tutorial on ABC rejection and ABC SMC for parameter estimation and model selection by Tina Toni and Michael P. H. Stumpf.
I can't work out how to calculate the weights for the SMC approach. Could anyone run through this step by step example to help me understand?

For t=0 Say I take a sample from a prior distribution based on uniform[1,3]

2.879864 2.684748 1.889464 2.945675 2.097058 1.003143 2.514226 2.242223 1.594360 2.764085 2.787965 1.052775 1.320575
2.108242 1.740970 2.214639 1.501381 2.234161 2.194186 1.331659

I then use these values to run some simulations and using a tolerance of 50% with a standard abc rejection method keep the following:

1.889464 1.003143 1.594360 1.052775 1.320575 1.740970 1.501381 1.331659 2.194186 2.097058

I add a bit of noise to each of these values using a Gaussian random walk:

1.9020030 0.9874041 1.6011953 1.0711497 1.2948880 1.7577606 1.5704593 1.3434156 2.2347718 2.1125749

How do I weight these in order to use them on a second set of simulations?

Best Answer

If you read the reference [4] in the tutorial, you can see a very detailed description of the ABC-SMC algorithm:

At the first iteration, the weights are all equal to one. At later iterations, because a new proposal which is a mixture of $K_t$ kernels is used instead of the prior, the weights are the ratios of the prior values to the mixtures values, as explained in S2.2.

There are several R packages implementing ABC-SMC, including easyABC, which, as an R package, should be easy (!) enough to install & experimen with.

Simply type install.packages("EasyABC") from within R.

There is for instance an ABC_sequential function that covers several sequential methods, including the one we proposed in Beaumont et al. (2009):

> ABC_Beaumont<-ABC_sequential(method="Beaumont", model=toy_model,
+ prior=toy_prior, nb_simul=n, summary_stat_target=sum_stat_obs,
+ tolerance_tab=tolerance)

"Method must be Beaumont, Drovandi, Delmoral, Lenormand or Emulation"

Related Solutions

Solved – Particle filter in Matlab – what is going wrong

The problem may stem from the fitness ratio you use in your code. The inverse of the distances may give very large differences in the resampling probabilities of the new generation of the samples. For example, let x stand for the rms values in your design and let them be

x = [0.1 1 1 3 10 50]

When you use the inverses, that leads to

y = (1./x) ./ norm(1./x)
y = 
    0.9896    0.0990    0.0990    0.0330    0.0099    0.0020

As you see, the first sample dominates the next generation. It may converge to wrong proposals very easily. Alternatively, you may use a sigmoid-like function. For example,

z = exp(-x) ./ norm(exp(-x))
z = 
    0.8659    0.3521    0.3521    0.0476    0.0000    0.0000

Now, the particles in the next generation will probably involve some samples similar to the second and third ones. So, the generation scheme will be more robust to erronous proposals.

Edit: Moreover, after dividing by the norms, the sum of the result is not 1. You might use sum instead.

y = (1./x) ./ sum(1./x)
y = 
    0.8030    0.0803    0.0803    0.0268    0.0080    0.0016

z = exp(-x) ./ sum(exp(-x))
z =
    0.5353    0.2176    0.2176    0.0295    0.0000    0.0000

Solved – ABC model selection

Nice question building on our work! Are you aware of the follow-up paper where we derive conditions on the summary statistic to achieve consistency in the Bayes factor? This may sound too theoretical but the consequence of the asymptotic results is quite straightforward:

Given a summary statistic $T$,

run an ABC algorithm based on $T$ for each model under evaluation ($i=1,..,I$) and estimate the parameters $\theta_i$ of those models by the ABC estimate $\hat\theta_i(T)$;
simulate the distribution of the statistic $T$ for each model and each estimated parameter, by a Monte Carlo experiment;
check whether the means $\mathbb{E}_{\hat\theta_i(T)}[T(X)]$ are all different by using step 2 with a sufficiently large number of iterations and, e.g., a t-test.

This procedure is not in the first version of the paper but should soon appear in the revised version

Best Answer

Related Solutions

Solved – Particle filter in Matlab – what is going wrong

Solved – ABC model selection

Related Question