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?