Infer parameters with ABC with non-uniform prior

Question

Edit (thanks to Xi'an) :

My data consists of $n$ realization of a specific experiment with $t$ time points and $2$ types of data measured in each time point. I summarized this data by computing the mean of each data type in each time points. This $2t$ values will be my summary statistics.

I have a model with 4 parameters that I wish to test to see if it can fit this data and predict another (further time points). However, my model is stochastic. This means that a given set of parameters will result in different simulated summary statistics. Because I'm not interested in the variance observed in the data (remember that my summary statistics are just the means of the observed data), for each parameter set I run my model 30 times and consider the average at each timepoint to be the simulated summary statistic.

Another characteristic of this model is that specific combinations of parameters will give rise to similar simulated summary statistics. Al thought in the limit of infinite number of simulation for each parameter set I expect the simulated summary statistics to be different (meaning that the parameters are identifiable), with a smaller number of runs and comparing with data that included errors the distance between the observed and simulated summary statistics could be very similar between different parameters. As an example, after running a simple ABC and selecting the best 5% I observed that between in the posterior distribution the correlation between 2 parameters is 0.95.

Because running the model is very time consuming, I can not explore the parameter space freely. My solution is to use an approach similar to sequential ABC, where at each step I use the results from the previous steps to compute the new area to search. I know that this approach does not guaranty convergence, but I think I had made it work (we can discuss this in another topic if anyone is interested). By the way, do you have a solution for this?

At this level of description, there is nothing wrong with using another parameterisation, provided it is one-to-one and that the prior is transformed by the Jacobian formula.

To be honest I did not understood this. What I did was the follwing: Assuming $p_i$ are the parameters in the model and $P_i$ are the transformed versions that I use to explore the parameter space:

$p_2 = e^{P_2}$

$p_3 = e^{P_3}$

$p_4 = \frac{1}{200(1+e^{P_4})}$

$p_1 = \frac{1}{200(1+e^{P_4})(1+e^{P_1})}$

The reasons for the previous equations are the limits I wanted to impose on the model parameters $p$ while being able to search in an unlimited scenario $P$.

Now, given a non-uniform sampling of the parameter space, how do I get:

• Point estimate of the parameters
• Confidence intervals for parameters
• A sampler from which I could draw parameters values that would obey the observed constrains

I selected the top 5% simulations and fit a multivariate normal distribution. By the way, for some reason the post-processing solution (Beaumont and neural network) implemented in abc (R) are not working fine (meaning they predict parameters that are worse than my solution.

How do I take into account the prior? Can you provided an R based solution?

Answer 1

Possible issues with the approach:

[Warning: the OP completely rewrote the question after I made those points, so most quotes relate to an earlier version!]

my model is stochastic. This means that a given set of parameters will result in different simulated summary statistics.

This is the standard setting for statistical inference, so I do not understand why this is stressed in the question. In an ABC algorithm, the simulated summary statistics are random variables and as such make the method valid (asymptotically).

My strategy was to simulate my model several times to each parameter set and use the average as the representation of the result for that parameter set.

This construction of a summary statistics must mimic the one for the original data, meaning you must have the same "several times" for the original data. When you mention 30 times, is $n=30$?

I didn't use a uniform prior to the parameters, but instead keep updating it to search is the areas with best fit.

Using an arbitrary distribution for the parameters is fine, provided one corrects for this change by an importance sampling weight. However, adapting the simulation of the parameters based on earlier ABC steps offers no convergence guarantee.

because I wanted to impose some constrains in the parameters I use transformed versions of then (like log or others).

At this level of description, there is nothing wrong with using another parameterisation, provided it is one-to-one and that the prior is transformed by the Jacobian formula.

This, plus something intrinsic to the model, lead to a correlation between some parameters (one could get similar results with different combinations of some parameters).

This is hard to understand. Correlation between the parameters is either a consequence of correlation in the prior or of dependences in the likelihood. Your parenthesised explanation is however off-key: finding two different values of the parameter with the same posterior value or with similar simulations has nothing to do with correlation.

the distance between the observed and simulated summary statistics could > be very similar between different parameters.

This is natural given the stochastic nature of the data. If one dataset could uniquely identify the parameter behind it, there would be no need for statistics.

Now, how do I compute the posterior distribution for my parameters?

ABC produces a sample from the ABC posterior. This posterior has a density that cannot be computed in closed form.

My simplistic strategy was the following (a) Accept the best 5% simulations (b) Fit a multivariate normal distribution to the parameters

The first bit is regular ABC, the second bit is a further approximation that does not seem necessary, as one could instead use a kernel approximation. The prior need be incorporated via further importance weighting.

As an option, I redid the previous step using as weighting the score

What is the meaning of score in this context? The distance? In which case the post-processing solution of Beaumont et al. (2003) should be considered.

Can you provide[d] an R based solution?

No, both on principle (I can help you understand better the approach, not solve the problem) and because your description is too vague to write an R code.

given a non-uniform sampling of the parameter space

You should provide a description of this non-uniform sampling on $p$ or $P$. From the available description, I can only guess you are using a certain proposal distribution, say $h(p)$, instead of the correct prior, say $g(p)$. In which case, each simulation should be [importance] weighted by $g(p)/h(p)$.