Investigate synthetic likelihood as a principled objective for stochastic models

[wshlavacek]

Ask

Investigate whether PyBNF should adopt the synthetic likelihood approach for fitting stochastic models, and report back with a recommendation. This issue does not propose a design. It asks for the investigation and for enough of a written answer that we can decide whether to fund the work.

Background

For a stochastic model we cannot write down the probability of the data given the parameters. Computing it would require summing over every trajectory the simulator might have taken. So today a stochastic fit in PyBNF minimizes a distance between simulated output and measured data, which is a reasonable heuristic but is not a statistical model of anything.

Wood proposed a way around this. Choose some summary statistics of the data, such as the mean at each time point, the variance, a peak height, or a time to peak. At a given set of parameter values, run the simulator several times and compute those summaries for each run. Assume the summaries are approximately drawn from a multivariate normal distribution, and estimate its mean and covariance from those runs. The likelihood is then the density of the observed summaries under that fitted normal. An intractable likelihood is replaced by a tractable one estimated by simulation.

Price and colleagues later gave the Bayesian treatment, which is the relevant version if we want to sample rather than only optimize.

• Wood, S. N. (2010). Statistically efficient inference for nonlinear ecological dynamic systems. Nature 466, 1102 to 1104.
• Price, L. F., Drovandi, C. C., Lee, A., and Nott, D. J. (2018). Bayesian Synthetic Likelihood. Journal of Computational and Graphical Statistics 27(1), 1 to 11.

Why this is worth looking at

It would give us a real likelihood for stochastic models. Our samplers (am, dream, hmc) assume there is a likelihood to evaluate. A stochastic model does not give us one, so on those models the samplers are working with a noisy stand in. With a synthetic likelihood they would be sampling something well defined, and the information criteria would mean something too.

It treats the noise as part of the model rather than as a nuisance. Issues #659, #660, and #661 all treat the randomness in a stochastic simulation as a problem to be averaged away. A synthetic likelihood treats the spread across repeated simulations as part of what is being tested. A model that predicts the right average but the wrong amount of variability is a wrong model, and this approach would notice. Right now PyBNF cannot. For many rule-based models the variability is the whole point of the study, so fitting only the average throws away the most interesting part of the measurement.

Much of the machinery is already planned. Step 3 of #660 calls for tracking the average and the spread of an objective across repeated simulations of the same parameter set, in order to decide whether a ranking can be trusted. A synthetic likelihood needs the average and covariance of summary statistics across repeated simulations. That is close to the same bookkeeping, so the two should at least be designed with each other in mind.

Questions the investigation should answer

• How would summary statistics be chosen and specified? This is the main open question, since the choice decides how much information is kept. Would a user declare them, or would we derive them from the existing experiment and observable definitions?
• How many repeated simulations are needed per parameter set to estimate a covariance matrix reliably, and how does that scale as summaries are added? This is the main cost.
• Does the assumption that the summaries are normally distributed hold well enough on our models? It is often reasonable when the summaries are averages, but it should be checked rather than assumed.
• How would this sit next to the existing objective functions? Is it a new objective a user selects, or something deeper?
• What would it take for the samplers to use it?
• Would designing Epic: finish the scatter search template, then make the reference set aware of simulation noise #660 step 3 with this in mind cost anything now? If it is cheap to keep the option open and expensive to retrofit, that is worth knowing before that work starts.

Related

• Epic: finish the scatter search template, then make the reference set aware of simulation noise #660 step 3, which builds the repeated simulation bookkeeping this would draw on.
• Best fit from a stochastic model is chosen from a single noisy simulation, so the reported objective is too good #659, which addresses a different consequence of the same underlying fact, that a single simulation of a stochastic model is a noisy measurement.
• Build a parameter recovery benchmark for fitting stochastic rule-based models #663, the proposed benchmark, which would be how we judge whether this actually fits better.