Simulations based noisy likelihood function

Question

I have a problem where I have a measured data vector $D$ with Gaussian uncertainties (covariance matrix $\Sigma$). I am now trying to model this data with a generative model with parameters $\phi$. The problem is that my model is simulation based.

For a given $\phi$, and given random seed s, and given number of particles in the simulation n, I can produce the model prediction $M(\phi, s, n)$ and this model prediction can be compared with the data vector D. The issue here is that the model prediction vector will have some uncertainty (due to randomness of the simulation and finite number of particles). We can assume that this uncertainty can be characterized by a Gaussian distribution with known Covariance $S(\phi, s, n)$

The question is what is the best way to write the likelihood and sample the posterior on $\phi$? Can anyone point at some papers that maybe discuss this kind of problem and how to think about it ? ( I am broadly aware of ABC techniques, but it is not obvious how any of those would be an improvement over what we do now)

For the moment the approach we use is basically ensure that n is large enough (to have smaller simulation errors) and then basically we just ignore the randomness in $M(\phi,s,n)$ and our log-likelihood is just $$ log L(\phi) = -\frac{1}{2}(D-M(\phi))^T \Sigma^{-1} (D-M(\phi))$$ which we sample using MCMC. (the resulting likelihood function is somewhat noisy, but the hope is that it all averages out in the resulting posterior on $\phi$)

Are there any alternatives to this ?

Answer 1

First, note that for a fixed seed $s$, you would not have a random output from the simulator, but always the same output, so the covariance $S$ would be zero. So I will remove the dependence of $S$ on $s$. Also, I will subsume $n$ into $\phi$, so that we are left with just the parameter vector $\phi$ and the covariance $S$ of the simulator only depends on $\phi$: $S = S(\phi)$.

The likelihood of given data D for a model $M(\phi)$, which is completely described by the parameters $\phi$, is by definition the probability of $D$ under this model: $$ logL(D | \phi) := \log p(D| M(\phi)). $$ So the covariance that you have to use in the formula for $logL$ is the covariance $S(\phi)$ of the model $M(\phi)$: fitting a model with MLE means you presume that the data $D$ is generated by your model $M(\phi)$ and then you find the $\phi$ that maximizes the belonging probability; you don't fit a model on data that you presume has been created by another model.

Now your model, which is the simulator $M(\phi)$, is Gaussian with both the mean $\mu(\phi)$ and the variance $S(\phi)$ being functions of $\phi$. The log-likelihood that you want to maximize is: $$ logL(D|\phi) = -\frac{1}{2} \big\langle(D - \mu(\phi)), \: S(\phi) (D - \mu(\phi))\big\rangle. $$ Note that $\Sigma$ is not entering into your model.

Since the model is a complex simulator, I presume that the functions $\mu(\phi)$ and $S(\phi)$ are not known. Thus, if you wanted to find the MLE estimation of $\phi$, you would have to resort to numerical optimization.

However, if you want the posterior of $\phi$, since you mentioned in the OP that you want to sample $\phi$), you need to devise an appropriate prior $p_{prior}(\phi)$ for $\phi$, compute the unnormalized posterior (define $L := exp(logL)$): $$ p_{post}(\phi) \propto L(D|\phi) p_{prior}(\phi), $$ and then sample from that, e.g. with some MCMC method.

Without any more information, it is impossible to say whether this is more or less appropriate than any of those many ABC methods.

Answer 2

If I understand your problem correctly, you have some data $D$ that is assumed to have a Gaussian distribution with an unknown mean $\mu(\phi)$ and known covariance $\Sigma$, namely

$$ D = \mu(\phi) + \mathcal N (0,\Sigma) $$

As well as a simulation output $M(\phi,n)$ which is also assumed to normally distributed with known covariance $S(\phi,n)$

$$ M(\phi,n) = \mu(\phi) + \mathcal N (0, S(\phi,n)) $$

(Where we can think of $\mu(\phi)$ as being equal to $\mu(\phi) = M(\phi,\infty)$ , since $S(\phi,\infty)=0$ )

But notice that this just implies that

$$ D = M(\phi,n) + \mathcal N (0,\Sigma + S(\phi,n)) $$

So your likelihood should be

$$ \log \mathcal L(\phi) = -\frac{1}{2}(D - M(\phi,n))^T(\Sigma + S(\phi,n))^{-1}(D - M(\phi,n) )$$

Now you can sample it using MCMC as you did, or try to optimize it using something like simulated annealing or Bayesian optimization, or any combination of the above (get close to the global minimum and sample around it).