Skip to content

Latest commit

 

History

History
161 lines (99 loc) · 7.62 KB

lotkavolterra_inference.rst

File metadata and controls

161 lines (99 loc) · 7.62 KB
Raw

Lotka-Volterra system inference

pyselfi

The Lotka-Volterra model

The Lotka-Volterra model (found in src/pyselfi/lotkavolterra/) defines a prior (lotkavolterra.prior) and a selfi child class (lotkavolterra.selfi). The lotkavolterra.selfi class allows optimization of prior hyperparameters (see sections II.C and II.E in Leclercq et al. 2019_, respectively).

The prior for time-dependent smooth functions

The class lotkavolterra.prior introduces a prior for time-dependent smooth functions. Let x(ti) be the number function of prey and y(ti) be the number function of predators at timesteps ti, with initially x0 ≡ x(t0) and y0 ≡ y(t0).

The prior on {{x(ti)},{y(ti)}} is a Gaussian with mean θ0 (a fiducial population number function) and covariance matrix S = αnorm2V ∘ K (where ∘ is the Hadamard product), such that: $\textbf{K} = \begin{pmatrix} \textbf{K}_x & \textbf{0} \\ \textbf{0} & \textbf{K}_y \end{pmatrix}$, $(\textbf{K}_z)_{ij} = -\dfrac{1}{2} \left( \dfrac{t_i-t_j}{t_\mathrm{smooth}} \right)^2$ for z ∈ {x,y}, $\textbf{V} = \begin{pmatrix} x_0 \textbf{u}\textbf{u}^\intercal & \textbf{0} \\ \textbf{0} & y_0 \textbf{u}\textbf{u}^\intercal \end{pmatrix}$, $(\textbf{u})_i = 1 + \dfrac{t_i}{t_\mathrm{chaos}}$.

The prior is characterized by 3 hyperparameters:

  • αnorm is the overall amplitude,
  • tsmooth is a typical time characterizing the smoothness of the population functions,
  • tchaos is a typical time characterizing the chaotic behaviour of the system.

It is used like this:

from pyselfi.lotkavolterra.prior import lotkavolterra_prior
alpha_norm = 0.02
t_smooth = 1.6
t_chaos = 8.2
prior=lotkavolterra_prior(t_s,X0,Y0,theta_0,alpha_norm,t_smooth,t_chaos)

t_s is a vector containing the support wavenumbers, X0 and Y0 are the initial conditions, and theta_0 is the desired prior mean (which should correspond to the expansion point for the effective likelihood). t_s and theta_0 shall be both of size S.

The SELFI model

The model is used like this:

from pyselfi.lotkavolterra.selfi import lotkavolterra_selfi
selfi=lotkavolterra_selfi(fname,pool_prefix,pool_suffix,prior,blackbox,
                          theta_0,Ne,Ns,Delta_theta,phi_obs,LV)

where LV is a LVsimulator (introduced by lotkavolterra_simulator).

The main computations (calculation of the effective likelihood) are done using:

selfi.run_simulations()
selfi.compute_likelihood()
selfi.save_likelihood()

Once the effective likelihood has been characterized, the prior can be later modified without having to rerun any simulation, for example:

alpha_norm = 0.05
t_smooth = 2.1
selfi.prior.alpha_norm=alpha_norm
selfi.prior.t_smooth=t_smooth
selfi.compute_prior()
selfi.save_prior()
selfi.compute_posterior()
selfi.save_posterior()

The prior can be optimized using a set of fiducial simulations theta_fid:

selfi.optimize_prior(theta_fid, x0,
                     alpha_norm_min, alpha_norm_max,
                     t_smooth_min, t_smooth_max,
                     t_chaos_min, t_chaos_max)
selfi.save_prior()

Check for model misspecification and score compression

The Mahalanobis distance between the posterior mean and the prior can be used to check for model misspecification (see Leclercq 2022_). It is obtained using

Mahalanobis = selfi.Mahalanobis_ms_check()

A score compressor for top-level parameters can be obtained without having to rerun any simulation. It is obtained using

selfi.compute_grad_f_omega(t, t_s, tres, Delta_omega)
selfi.compute_Fisher_matrix()
omega_tilde = selfi.score_compression(phi_obs)

For more details see the example notebook (selfi_LV.ipynb) and check the API documentation (pyselfi.lotkavolterra.prior, pyselfi.lotkavolterra.selfi).

Inference of top-level parameters

The Fisher-Rao distance between two points in top-level parameter space can be obtained using:

Fdist = selfi.FisherRao_distance(omega_tilde1, omega_tilde2)

It can be used for example, within likelihood-free rejection sampling (Approximate Bayesian Computation, ABC). See the example notebook ABC_LV.ipynb.

Blackbox using lotkavolterra_simulator

pyselfi_examples.lotkavolterra.model

A blackbox using the lotkavolterra_simulator code is available in the package pyselfi_examples installed with the code. It is defined as an instance of the class blackbox_LV.blackbox (see Leclercq 2022_ for details):

from pyselfi_examples.lotkavolterra.model.blackbox_LV import blackbox
# correct data model
bbA=blackbox(P, X0=X0, Y0=Y0, seed=seed, fixnoise=fixnoise, model=0,
             Xefficiency=Xefficiency0, Yefficiency=Yefficiency0,
             threshold=threshold, mask=mask, obsP=obsP,
             obsQ=obsQ, obsR=obsR0, obsS=obsS, obsT=obsT)
# misspecified data model
bbB=blackbox(P, X0=X0, Y0=Y0, seed=seed, fixnoise=fixnoise, model=1,
             mask=mask, obsR=obsR1)

See the API documentation for the description of arguments, and the file src/pyselfi_examples/lotkavolterra/model/setup_LV.py for an example.

pyselfi_examples.lotkavolterra.model.blackbox_LV.blackbox

Synthetic observations to be used in the inference can be generated using the method make_data:

phi_obs = bbA.make_data(theta_true, seed=seed)

The blackbox object also contains the method compute_pool with the necessary prototype (see this page), as well as a method evaluate to generate individual realizations from any vector of input parameters θ.

A full worked-out example using the Lotka-Volterra blackbox is available in this notebook.