Parameter inference

Parameterizing complex models is still surprisingly difficult. it is both fortunate that many different tools exist today to approach model parameterization. This task has a variety of names with slightly different associated meaning, but they all address the same task. Namely, estimating the set of parameters that best describe some data \(Y\) under the (mathematical) assumptions of the model. This process can be called parameter inference, parameter estimation, optimization, calibration, fitting. While there are some subtle differences between the terms, they can be used more or less interchangibly. Relevant differences arise from the tools being used. The terms could be classified based on the probabilistic nature of the approach. While calibration and optimization are rather associated with deterministic models, which should have a unique optimal solution; estimation and inference are rather associated with approaches where the data are assumed to originate from a probability distribution.

For simplicity we will use the term parameter inference, because we interpret by its meaning from the fields of logic and reasoning, where it describes the process from moving from premises to logical consequences.

Statistical inference uses mathematics to draw conclusions in the presence of uncertainty. This generalizes deterministic reasoning, with the absence of uncertainty as a special case. Statistical inference uses quantitative or qualitative (categorical) data which may be subject to random variations.

Supported inference backends

Backend	Supported Algorithms	Inference	Hierarchical Models
numpyro	Markov Chain Monte Carlo (MCMC), Stochastic Variational Inference (SVI)	✅	✅
pymoo	(Global) Multi-objective optimization	✅	plan
pyabc	Approximate Bayes	✅	plan
scipy	Local optimization (minimize)	dev	plan
pymc	MCMC	plan	plan
sbi	Simulation Based Inference (in planning)	hold	hold
interactive	interactive backend in jupyter notebookswith parameter sliders	✅	plan

Why use pymob.inferer?

The goal of pymob is to reduce the overhead for exploring, parameterizing and comparing models so that the focus can remain on model development. Because after all, we would like to use models to anwer questions.

pymob.inferer only enforces minimal conventions, which allows the interoperability between different frameworks. At the same time it allows you to take the models you have already written and just plug them into the framework and use it to do inference on it.

• Observations/data have to be provided as one xarray.Dataset.

• Simulations have to be wrapped inside a function, which converts the simulation output to a dictionary of results, so that results can be converted to xarray.Datasets

• It is recommended that priors are specified following the scipy syntax, because functionality is implemented to translate scipy priors to different probabilistic programming languages (PPL). This is no strict requirement, but using PPL specific notation, sacrifices the possibility of interoperability between frameworks

Parameter inference

The components of parameter inference or optimization are usually the same or very similar. We need:

• priors

• a deterministic function (which can be an identity)

• an error model (or error distribution)

If you are thinking I don’t need any priors for an optimization task, you are somehow correct, but we could also say that the specified bounds (lb, ub) on a parameter act as a prior (with the corresponding distribution being \(Uniform(lb, ub)\)). Even an unconstrained parameter is a Uniform between \([-\infty, \infty]\). Despite pymob sticks to the convention and let’s you specify bounds or priors as you like.

Example: The hierarchical Lotka-Volterra case-study

The hierarchical Lotka-Volterra case study assumes that observations were taken from two rabbit populations (species: Cottontail and Jackrabbit), observed in 2010, 2011 and 2012 in different valleys, which define the number of in-treatment replicates.

The case study is designed to study the possibilities of recovering species level parameters from nested error structures.

We set up the case study like this:

from pymob import Config

config = Config()
config.case_study.name = "lotka_volterra_case_study"
config.case_study.scenario = "test_hierarchical"
config.case_study.simulation = "HierarchicalSimulation"
config.import_casestudy_modules(reset_path=True)
Simulation = config.import_simulation_from_case_study()

sim = Simulation(config)
sim.setup()

Creating parameters and priors

A parameter in pymob for instance for the hierarchical Lotka-Volterra case study can be specified like this:

from pymob.sim.parameters import Param
alpha_species = Param(
    value=0.5, free=True, hyper=True, min=None, max=None,
    prior="norm(loc=[[1],[3]],scale=0.1,dims=('rabbit_species','experiment'))"
)

The prior assumes that we have a lot of prior knowledge of the alpha parameter of the different species.

Distribution shapes and dimensionality

Take good care to specify the shape of your priors correctly. Dimensions are broadcasted following the normal rules of numpy. The above means, in the 0-th dimension (axis), we have two different assumptions loc=1, and loc=3. The 0-th dimension is the dimension of the rabbit species. And the parameters are simply broadcasted to match the shapes of any other dimension.

The specification loc=[1,3] would be understood as [[1,3]] and be interpreted as the experiment dimension. If the dimensions have different shapes (i.e. different number of unique elements), this will result in an error. If the dimensions have the same size, this will lead to serious problems, because the above will work, but lead to incorrect results.

Ideally, the dimensionality is so low that you can be specific about the priors. I.e.: loc=[[1,1,1],[3,3,3]]. This of course expects you know about the dimensional sizes of the prior (i.e. the unique elements of the dimensions).

We can use the hyperprior and define a prior alpha which actually parameterizes our model.

from pymob.sim.parameters import Param
alpha = Param(
    value=0.5, free=True, hyper=False,
    prior="lognorm(s=0.1,scale=alpha_species[rabbit_species_index, experiment_index],dims=('id',))"
)

This prior samples from a Lognormal distribution, but the scale (log-mean) is parameterized by the hyper prior alpha_species. In addition, the indices rabbit_species_index and experiment_index are responsible for selecting the samples from the hyper-prior alpha_species so that they correctly map to the ids of the different observations. Pymob automatically creats these indices for you, if you provide the experimental layout.

Parsing priors as strings

The prior is currently passed as a string, which is somewhat inconvenient, but the problem is that in more complex distributions, some variables are only known when they are sampled.

This means we have to postpone evaluation of these priors to a later point. Pymob handles this by parsing the string into an Abstract Syntax Tree, where only general syntax problems are detected. Any problems with the indices, dimensions will only occurr, when samples are taken from the prior.

Generating the experimental layout

sim.define_observations_replicated_multi_experiment(n=12) 
y0 = sim.parse_input("y0", drop_dims=["time"])
sim.model_parameters["y0"] = y0

Sampling from the generated distribution

from pymob.inference.scipy_backend import ScipyBackend

inferer = ScipyBackend(simulation=sim)
theta = inferer.sample_distribution()

TODO

• Create IO tools for importing existing models and data

• Write documentation with examples