Combining multiple posterior distributions I am new to Bayesian statistics, and thus have problems coming up with a solution for the following problem:
Using Approximate Bayesian Computation (ABC), I generate a posterior distribution from one set of observed summary statistics, and a set of about a million simulated summary statistics. The prior distribution is the variable value of the model parameter used in the simulations. The posterior distribution is used to obtain point estimates (e.g., mode, mean, etc.) of the underlying model parameter.
The biological problem I am working on is that the true parameter value cannot be estimated from one observation only. Hence, I have $N$ observed data sets, and obtain $N$ posterior distributions. In ABC, I use the same simulated data for each observation data set.
Now, my problem is to combine the $N$ posterior distributions in a way to estimate the true parameter value, thus capturing the information contained in all $N$ densities.
Any ideas?
 A: Given a measurable space $\Lambda$ of parameter values and a measurable space $X$ of observable values, a Bayesian model can be considered a function $\phi : \Lambda \times X \to \mathbb{R}$ where the restricted function at each parameter set is a probability density function on the observable variables. That is,
$$ \forall \lambda \in \Lambda, x \in X : \phi(\lambda, x) \geq 0 \\ \forall \lambda \in \Lambda : \int_X \phi(\lambda, x)\, d x = 1 $$
Most often you'll have $\Lambda \subseteq \mathbb{R}^m$ and $X \subseteq \mathbb{R}^n$; we're just representing all the parameters at once and all the variables found in a single observation at once.
A postulated distribution for the parameters is likewise a probability density function on just the parameters. That is, a function $\alpha : \Lambda \to \mathbb{R}$ with
$$ \forall \lambda \in \Lambda : \alpha(\lambda) \geq 0 \\ \int_\Lambda \alpha(\lambda) \, d\lambda = 1 $$
Then given a model $\phi$ and an observation $x$, a Bayesian update from prior distribution $\alpha_\mathrm{pre}$ to posterior distribution $\alpha_\mathrm{post}$ looks like
$$ \alpha_\mathrm{post}(\lambda) = \frac{\alpha_\mathrm{pre}(\lambda) \phi(\lambda,x)}
{\int_\Lambda \alpha_\mathrm{pre}(\mu) \phi(\mu,x)\, d\mu} $$
(If the denominator is zero, there was a bad assumption in the model and/or the prior: the set of parameters $\mu$ with $\alpha_\mathrm{pre}(\mu)>0$ where the prior considers parameter values possible and the set of $\mu$ with $\phi(\mu,x)>0$ where the model considers observation $x$ possible don't intersect.)
When we have multiple observations $x_1, \ldots, x_k$, this generalizes to
$$ \alpha_\mathrm{post}(\lambda) = \alpha_\mathrm{pre}(\lambda) \frac{\prod_{i=1}^k \phi(\lambda,x_i)}
{\int_\Lambda \alpha_\mathrm{pre}(\mu) \left[\prod_{i=1}^k \phi(\mu,x_i)\right] d\mu} $$
This is just what we get when updating parameter distributions in a sequence $\alpha_\mathrm{pre} = \alpha_0 \to \alpha_1 \to \cdots \to \alpha_k = \alpha_\mathrm{post}$, updating for one additional observation $x_i$ at a time. (Happily, the end result makes it obvious the order of those observations didn't matter.)
