internals.fitting

internals.fitting

Fitting sims to data: likelihoods, particle filters, samplers, and what sim.fit() runs

These exist so that a model can be fitted, with its uncertainty, in the same way whichever solver runs it. Parameters are chosen by giving them bounds (ss.par), in declared modules or in the parameters of modules written as code. Some methods are for some solvers only: gradients come from JAX, and particle filters need a stochastic compartmental solver.

Imported when it is first used (by ss.fitting.X, sim.fit(), or sim.sensitivity()), not by import starsim.

Used with: agents and compartments; declared and coded modules

Classes

Name	Description
Fit	The result of fitting a sim to its data; see `sim.fit()`
ModelComparison	The result of comparing models by how well they predict held-out data; see `ss.fitting.compare_models()`
Optimization	The result of allocating spending to maximize (or minimize) a result; see `sim.optimize()`
ParticleFilter	An estimate of the likelihood of a stochastic sim’s data, from a bootstrap particle filter (see `sim.fit()` with `method='pmcmc'`)
Swarm	The parameter values of each particle of IF2, iterated filtering, and the rates of the run that they give (see `sim.fit()` with method=‘if2’)
maximize	The goal of `sim.optimize()`: a result to make as large as possible
method	A fitting method and its settings, for `sim.fit(method=...)`
minimize	The goal of `sim.optimize()`: a result to make as small as possible

Fit

internals.fitting.Fit(
    sim,
    before,
    pars,
    table,
    trials,
    method,
    objective,
    objective_before,
    measure,
    samples,
    diagnostics,
    converged,
    why,
    scores,
    study=None,
    trajectories=None,
)

The result of fitting a sim to its data; see sim.fit()

Attributes

Name	Type	Description
sim	`ss.Sim`	the calibrated sim: a copy of the original with the best-fit values, already run
pars	`dict`	the best-fit value of each fitted parameter; can be passed to `ss.Sim(pars=...)`
table	`df`	parameter, initial value, best fit, and bounds (and with `method='if2'`, the lowest and highest of the chains’ estimates, as `spread_lower` and `spread_upper`)
before	`ss.Sim`	the sim with the initial values, already run
objective	`float`	what the fit minimized, for the calibrated sim: its loss, or its negative log-likelihood (or log-posterior) if the fit used likelihoods; for a stochastic sim, the particle filter’s estimate of it (with Optuna, where the filter can run the sim), or else the mean over several runs of each run’s
measure	`str`	what `objective` measures, e.g. ‘negative log-likelihood’, or for a stochastic sim ‘negative log-likelihood (estimated by a particle filter of 200 particles)’ or ‘negative mean log-likelihood of 4 runs’ (which is not the likelihood: the mean of the logs, not the log of the mean)
objective_before	`float`	the same, before the fit
converged	`bool`	whether the searches from the best few points of the scan that are apart end at the same value (for SciPy; with a single search, SciPy’s own stopping rule), stopped improving (for Optuna: by less than 1% over the last quarter of the trials), or has chains that agree (for NUTS and particle MCMC: an R-hat of at most 1.01, and no divergent draws; for IF2: chains from different starts that end within 2 of the best log-likelihood); None for ABC-SMC, whose posterior is approximate, by an amount that depends on the tolerance it reached, with no rule for when it is close enough
why	`str`	what `converged` is based on
scores	`list`	the score of the best-fit values in each run (one for a deterministic sim, and one, the particle filter’s, where it scored the points); `objective` is their mean, so it is the score that was minimized
trials	`df`	every evaluation: the values, and the `objective` that the method minimized (what `measure` says; with `method='if2'`, also each pass of each chain: the middle of the swarm’s values, and minus the log-likelihood of the pass, whose values were perturbed)
samples	`df`	with `method='nuts'`, `'pmcmc'`, or `'abc-smc'`, the draws from the posterior: one row per draw, with its chain, its number, and a column for each parameter
trajectories	`list`	for a stochastic compartmental sim with `method='pmcmc'` or `'abc-smc'`, the run of the sim that goes with each draw of `samples` (one per row), so that each pair is a draw from the posterior of the values and the epidemic together: with `'pmcmc'`, a run drawn from the particle filter, up to the last data point (see `ss.fitting.ParticleFilter`); with `'abc-smc'`, the run that was accepted. Each is, by the module that runs the compartments, the time point it ends at (`ti`), the compartments there (`y`), and by disease how many were in each state (`stocks`) and the flow along each transition (`flows`) at each time point before it; a draw that the chain stayed at shares its run with the draw before. `fit.project()` continues them. (For a deterministic sim, the run at each draw is its trajectory, so none are kept)
study	`optuna.Study`	with `method='optuna'`, Optuna’s study, e.g. for its plots (`optuna.visualization.plot_param_importances(fit.study)`)
diagnostics	`dict`	what else the method reports: `divergences` (with `method='nuts'`, the number of draws that followed a divergent trajectory, which suggests the sampler missed part of the posterior), `acceptance` (with `method='pmcmc'`, the fraction of proposals that each chain accepted after its warm-up; with `'abc-smc'`, that of its last generation), `tolerance` and `generations` (with `'abc-smc'`, the tolerance of its last generation, which sets how close its posterior is, and the number of generations run), and `band` (with a posterior, for each data series, the 2.5th and 97.5th percentiles of the model’s value over draws from it, at each step, or for a flow with data at a longer period than the step, of its totals over the data’s periods; with `'pmcmc'`, each from a run with a seed of its own, so the band includes chance)

Methods

Name	Description
plot	Plot the results that have data against it, before and after fitting; a flow with data at a longer period than the step (e.g. weekly cases) is shown as its totals over the data’s periods, as the fit compares them
project	Run the fitted model forward under each scenario, with the uncertainty in the fitted values: one run per draw from the posterior, each with a seed of its own (or runs of the best fit, each with its own seed, for a method that gives none), with the same draws and seeds in every scenario
summarize	Print the goodness of fit before and after, and return the estimates (a dataframe), which are printed too in a script (see `ss.utils.return_table()`)
to_df	Every evaluation: the parameter values tried, and the objective at each

plot

internals.fitting.Fit.plot(residuals=False, **kwargs)

Plot the results that have data against it, before and after fitting; a flow with data at a longer period than the step (e.g. weekly cases) is shown as its totals over the data’s periods, as the fit compares them

Parameters

Name	Type	Description	Default
residuals	`bool`	whether to add the residuals of the fitted model (the data minus the model, as the fit compares them) in axes below each result	`False`
kwargs	`dict`	see `ss.plot_args()`, e.g. `figsize=(8, 4)`

Example

fit = sim.fit()
fit.plot(residuals=True) # A pattern in the residuals points to the structure of the model

project

internals.fitting.Fit.project(
    scenarios=None,
    n_runs=None,
    rand_seed=None,
    restart=None,
    **kwargs,
)

Run the fitted model forward under each scenario, with the uncertainty in the fitted values: one run per draw from the posterior, each with a seed of its own (or runs of the best fit, each with its own seed, for a method that gives none), with the same draws and seeds in every scenario

With a run kept for each draw (fit.trajectories), each run continues the draw’s own, from where it ends (with a seed of its own for what follows): up to there it is the epidemic that the draw’s values give with the data, and the scenarios differ only after it. With restart=True, each draw is run from the start instead, as for a fit that keeps no runs, so its course is not conditioned on the data, only its values are; this is also what happens, with a warning saying why, when the runs can’t be continued (a scenario the kernels don’t run whole, e.g. an intervention they do not take, or different time points). Either way, for a deterministic sim, the runs are the same.

Parameters

Name	Type	Description	Default
scenarios	`dict`	for each scenario, what changes: None, an intervention or a list of them, or a dict of arguments for `sim.new()` (default: the fitted sim alone)	`None`
n_runs	`int`	the number of runs of each scenario (default: 100; but one for the best fit of a deterministic sim (an ODE), whose runs would all be the same)	`None`
rand_seed	`int`	the seed of the draws	`None`
restart	`bool`	whether to run each draw from the start, rather than continue its run in `fit.trajectories` (if the fit kept them); by default, continue where the runs can be continued	`None`
kwargs	`dict`	passed to `sim.new()` for every scenario, e.g. a later `stop=`

Returns an ss.MultiSim, whose results have a run dimension labeled by scenario, with the draws as runs of each scenario (see msim.to_df(quantiles=...), whose quantiles are over the draws). Runs that continue those of the draws are kept as results only (with no results by group), since a sim cannot start partway through.

Examples

fit = sim.fit(method='abc-smc')
proj = fit.project(dict(baseline=None, vaccine=ss.vaccinate(...)), stop='2021-01-01')
df = proj.to_df(quantiles=[0.05, 0.5, 0.95])

summarize

internals.fitting.Fit.summarize()

Print the goodness of fit before and after, and return the estimates (a dataframe), which are printed too in a script (see ss.utils.return_table())

to_df

internals.fitting.Fit.to_df()

Every evaluation: the parameter values tried, and the objective at each

ModelComparison

internals.fitting.ModelComparison(
    fits,
    table,
    holdout,
    measure,
    splits,
    clear=True,
)

The result of comparing models by how well they predict held-out data; see ss.fitting.compare_models()

Attributes

Name	Type	Description
fits	`dict`	each model’s `ss.fitting.Fit` (a list of them, one per fold, if there are several folds)
table	`df`	the fitted parameters, the error on the data fitted to, the prediction error, and the rank of each model
holdout	`str`	a description of what was held out
measure	`str`	how the fit error and the prediction error are measured (the same for every model)
splits	`list`	the data fitted to and the data held out, for each fold
clear	`bool`	whether the best model is clearly better, given the spread across folds

Methods

Name	Description
plot	Plot each model’s fit and prediction against the data, with the held-out data hollow; kwargs are as for `ss.plot_args()`
summarize	Print what was held out, and return the table, which is printed too in a script (see `ss.utils.return_table()`)
to_df	The table, as a dataframe

plot

internals.fitting.ModelComparison.plot(**kwargs)

Plot each model’s fit and prediction against the data, with the held-out data hollow; kwargs are as for ss.plot_args()

Parameters

Name	Type	Description	Default
kwargs	`dict`	passed to `ss.plot_args()`

summarize

internals.fitting.ModelComparison.summarize()

Print what was held out, and return the table, which is printed too in a script (see ss.utils.return_table())

to_df

internals.fitting.ModelComparison.to_df()

The table, as a dataframe

Optimization

internals.fitting.Optimization(**kwargs)

The result of allocating spending to maximize (or minimize) a result; see sim.optimize()

Attributes

Name	Type	Description
alloc	`dict`	the best allocation, by program
value	`float`	the result with it
total	`float`	its total (with `target=`, the least that reaches the target)
alloc_before	`dict`	the allocation before
value_before	`float`	the result with it
sim	`ss.Sim`	the sim run with the best allocation
n_runs	`int`	the number of runs the search took

Methods

Name	Description
summarize	Print the result before and after, and return the allocation before and after (a dataframe), which is printed too in a script (see `ss.utils.return_table()`)

summarize

internals.fitting.Optimization.summarize()

Print the result before and after, and return the allocation before and after (a dataframe), which is printed too in a script (see ss.utils.return_table())

ParticleFilter

internals.fitting.ParticleFilter(sim, likelihood=None, weights=None)

An estimate of the likelihood of a stochastic sim’s data, from a bootstrap particle filter (see sim.fit() with method='pmcmc')

Many runs of the sim (the particles, ss.solvers.Particles) go forward together. At each data point, each is weighted by the likelihood of the data given its results, and when the weights become uneven, the runs are drawn again in proportion to them, so that runs that have strayed from the data are dropped and those that follow it are copied. The product over the data of the mean weight is an unbiased estimate of the likelihood (Andrieu, Doucet & Holenstein 2010), which is what particle MCMC needs.

Each declared disease is filtered on its own, since the diseases of a sim that the kernels run whole do not affect each other, except those tracked jointly (see ss.solvers.JointStates), which are filtered together, as the runs of their program. An outcome is compared with the data as its expected number (rather than drawn, as a run of the sim draws it), as the likelihood would compare it.

After each call, path is one run drawn from the filter: a particle at the last data point, drawn by its weight, and the particles it was copied from back to the start (Andrieu, Doucet & Holenstein 2010), as particle MCMC keeps with each draw (see Fit.trajectories). For each disease filtered (or diseases tracked jointly, by the first’s name), it has the time point after the last step that its data read, ti; the compartments there, y; and, by disease, how many were in each state (stocks) and the flow along each transition (flows) at each time point before it (totals over the groups), as ss.solvers.Compartmental.run_many() takes them to continue the run.

Each particle has random numbers of its own, so the particles are shared between threads (see ss.options.threads; one in the workers of ss.parallel() and of a calibration), and the estimate does not depend on how many there are.

Parameters

Name	Type	Description	Default
sim	`ss.Sim`	the sim, initialized; it is restarted with each set of values before the filter is called	required
likelihood	`dict`	the likelihood of each data series, if not declared on its outcome	`None`
weights	`dict`	a multiplier on each series’ log-likelihood (default 1)	`None`

Methods

Name	Description
derive	The values of a result worked out from the others, at the steps given: an outcome from the flows (and just the outcomes it needs), or anything else as a run works it out
filter	The estimate of the log-likelihood of the data of one disease (or of the diseases of joint states, run by its program)
keep_path	Keep the run of particle `b` in `path`: its stocks and flows, by tracing the particles it was copied from at each
reach	For a result worked out from the others: how many steps before a time point its value there depends on the run, or None if on the whole run before it
reader	Where a result is in the stocks or flows of a run, if it is read straight from them: (‘stocks’ or ‘flows’, indices to sum)
reason	Why the particle filter cannot run a sim (apart from any data series without a likelihood), if it cannot
sources	An outcome, and the outcomes it is worked out from

derive

internals.fitting.ParticleFilter.derive(s, stocks, flows)

The values of a result worked out from the others, at the steps given: an outcome from the flows (and just the outcomes it needs), or anything else as a run works it out

filter

internals.fitting.ParticleFilter.filter(js, n, rng, joint=None, swarm=None)

The estimate of the log-likelihood of the data of one disease (or of the diseases of joint states, run by its program)

keep_path

internals.fitting.ParticleFilter.keep_path(
    js,
    runs,
    b,
    resampled,
    joint=None,
    by_group=False,
)

Keep the run of particle b in path: its stocks and flows, by tracing the particles it was copied from at each resampling back to the start, and its compartments now. At a resampling, the runs keep the steps from keep on (which only grows), so a particle’s own steps are those from the keep of the last resampling before it, to that of the next

reach

internals.fitting.ParticleFilter.reach(engine, name)

For a result worked out from the others: how many steps before a time point its value there depends on the run, or None if on the whole run before it

reader

internals.fitting.ParticleFilter.reader(engine, name)

Where a result is in the stocks or flows of a run, if it is read straight from them: (‘stocks’ or ‘flows’, indices to sum)

reason

internals.fitting.ParticleFilter.reason(sim)

Why the particle filter cannot run a sim (apart from any data series without a likelihood), if it cannot

sources

internals.fitting.ParticleFilter.sources(engine, name)

An outcome, and the outcomes it is worked out from

Swarm

internals.fitting.Swarm(
    sim,
    fitted,
    set_pars,
    to_x,
    z0,
    sd=0.05,
    cooling=0.5,
    n_points=1,
)

The parameter values of each particle of IF2, iterated filtering, and the rates of the run that they give (see sim.fit() with method=‘if2’)

Each particle carries values of its own, on the unbounded scale that the fit searches on, which a random walk perturbs at the start of each pass of the filter and before each data point, by steps that shrink geometrically, to cooling of their first size after 50 passes (as pomp’s mif2() and camdl do; Ionides et al. 2015). The rates are worked out from each particle’s values by the model’s own expressions, so the kernels run each particle with its own, and they are checked against the sim’s own.

Parameters

Name	Type	Description	Default
sim	`ss.Sim`	the sim, initialized, with the data of one declared disease	required
fitted	`dict`	the fitted parameters, by name	required
set_pars	`func`	restarts the sim with an array of values	required
to_x	`func`	gives the values (an array, or particles by values) from the search scale, and their places between the bounds	required
z0	`array`	the values the sim is restarted with now, on the search scale	required
sd	`float`	the first size of a step of the random walk, on the search scale	`0.05`
cooling	`float`	the fraction of the first size that a step has after 50 passes	`0.5`
n_points	`int`	the data points in a pass	`1`

Methods

Name	Description
convert	A fitted value as the engine uses it: in the sim’s unit, as `ss.gradients.Density.convert()`
estimate	The values at the middle of the swarm (its mean on the search scale, as pomp’s `coef()` of a `mif2()` fit)
fracs	The time fraction between places of every particle, as `ss.kernels.advance()` takes it (particles by entries; default, the swarm’s)
perturb	A step of the random walk for every particle
resample	The values of the particles drawn again
start	Start every one of n particles at the same values (on the search scale)
tables	The parts of the program that vary in time, as `ss.kernels.advance()` takes them: for every particle, with its own values (particles by values; default, the swarm’s), at some rows (of those that vary in time)

convert

internals.fitting.Swarm.convert(i, val)

A fitted value as the engine uses it: in the sim’s unit, as ss.gradients.Density.convert()

estimate

internals.fitting.Swarm.estimate()

The values at the middle of the swarm (its mean on the search scale, as pomp’s coef() of a mif2() fit)

fracs

internals.fitting.Swarm.fracs(x=None)

The time fraction between places of every particle, as ss.kernels.advance() takes it (particles by entries; default, the swarm’s)

perturb

internals.fitting.Swarm.perturb(rng)

A step of the random walk for every particle

resample

internals.fitting.Swarm.resample(picks)

The values of the particles drawn again

start

internals.fitting.Swarm.start(z, n)

Start every one of n particles at the same values (on the search scale)

tables

internals.fitting.Swarm.tables(x=None, rows=slice(None))

The parts of the program that vary in time, as ss.kernels.advance() takes them: for every particle, with its own values (particles by values; default, the swarm’s), at some rows (of those that vary in time)

maximize

internals.fitting.maximize(result, at=None)

The goal of sim.optimize(): a result to make as large as possible

Parameters

Name	Type	Description	Default
result	`str` / `func`	a result of a module (e.g. `'all_vs'`, or `'hiv.n_vs'`), or a function of the run sim that returns a number	required
at	`float` / date	the time of the result (default: the end), read as any time is (see `ss.connectors.to_time()`: a bare number is a year if the sim has a calendar); without it, a function of the sim decides	`None`

Example

best = sim.optimize(ss.fitting.maximize('all_vs', at=2030), spend=['testing', 'art'])

method

internals.fitting.method(name, **settings)

A fitting method and its settings, for sim.fit(method=...)

sim.fit(method='pmcmc') uses a method with its default settings; ss.fitting.method() gives them, e.g. sim.fit(method=ss.fitting.method('pmcmc', particles=500, chains=2)). A setting the method does not take is an error that lists those it does.

Parameters

Name Type Description Default

name str the method (see sim.fit()): ‘nelder-mead’, ‘powell’, ‘l-bfgs-b’, ‘slsqp’ (SciPy), ‘optuna’, ‘nuts’, ‘pmcmc’, ‘if2’, or ‘abc-smc’ required

settings

Name	Type	Description	Default
name	`str`	the method (see `sim.fit()`): ‘nelder-mead’, ‘powell’, ‘l-bfgs-b’, ‘slsqp’ (SciPy), ‘optuna’, ‘nuts’, ‘pmcmc’, ‘if2’, or ‘abc-smc’	required
settings	`dict`	the method’s settings; those of each method, with their defaults: SciPy methods: `scan` (the points spread over the bounds to try first: 10 per parameter, at most 100), `restarts` (the searches from the best of them that are apart: 3, or 1 for one parameter of a deterministic sim), and SciPy’s own options, e.g. `maxiter=500` (see `scipy.optimize.show_options()`) `'optuna'`: `n_trials` (default: 50 per fitted parameter, at least 100), `particles` (for a stochastic sim that the particle filter can run, with a likelihood for every data series: the particles of the filter that scores each point, 200; see `ss.fitting.ParticleFilter`), `workers` (processes that run trials at once, sharing what they learn: 1), `storage` (a file name, a database URL, or an Optuna storage, where the trials are kept so that calling the fit again resumes it), and Optuna’s own keywords for `study.optimize()`, e.g. `timeout=600` `'nuts'`: `draws='auto'` (until the chains agree, with an R-hat of at most 1.01, and every value has 400 effective draws; or a number), `warmup=300`, `chains=4`, `scan` (each chain starts from one of the best points of a scan that are apart), and for `draws='auto'`: `ess=400`, `rhat=1.01`, `max_draws=5000`, `batch=100` `'pmcmc'`: `draws='auto'` (as for NUTS), `warmup=1000`, `chains=4`, `max_draws=20000`, `batch=250`, `particles` (of the particle filter: by default chosen by each chain in its warm-up; see `ss.fitting.pmmh()`; none for a deterministic sim, whose likelihood is exact), `scan` (each chain starts from one of the best points of a scan that are apart), `parallel=True` (to run the chains in parallel processes) `'if2'`: `particles=1000` (of the particle filter), `iterations=50` (passes of the filter in each chain), `cooling=0.5` (the fraction of the first size of the random walk’s steps that is left after 50 passes, as pomp’s `cooling.fraction.50`), `sd=0.05` (that first size, on the scale searched on), `chains=4` `'abc-smc'`: `samples=200` (the parameter values kept in each generation, which are the draws of the posterior), `generations=8`, `quantile=0.5` (of the distances, which sets the next tolerance), `tolerance` (the distance to stop at) for a stochastic sim, how each point is scored (SciPy methods and Optuna; `n_runs` also for particle MCMC’s start): `n_runs` (runs per point: 4 for a stochastic compartmental sim, 1 for agents and for a deterministic one; every point is scored by the mean over the same seeds of each run’s score; for Optuna, giving it scores by runs in place of the particle filter), and `reseed` (whether each point takes fresh seeds: by default only for Optuna, whose search then averages over chance) for the methods with gradients (‘l-bfgs-b’, ‘slsqp’, ‘nuts’): `require_jax=True` (whether the gradients have to be taken with the log-likelihood as one function in JAX; with `False`, something JAX does not cover yet is taken with Python in the loop instead, which is much slower)

dict

the method’s settings; those of each method, with their defaults:

SciPy methods: scan (the points spread over the bounds to try first: 10 per parameter, at most 100), restarts (the searches from the best of them that are apart: 3, or 1 for one parameter of a deterministic sim), and SciPy’s own options, e.g. maxiter=500 (see scipy.optimize.show_options())
'optuna': n_trials (default: 50 per fitted parameter, at least 100), particles (for a stochastic sim that the particle filter can run, with a likelihood for every data series: the particles of the filter that scores each point, 200; see ss.fitting.ParticleFilter), workers (processes that run trials at once, sharing what they learn: 1), storage (a file name, a database URL, or an Optuna storage, where the trials are kept so that calling the fit again resumes it), and Optuna’s own keywords for study.optimize(), e.g. timeout=600
'nuts': draws='auto' (until the chains agree, with an R-hat of at most 1.01, and every value has 400 effective draws; or a number), warmup=300, chains=4, scan (each chain starts from one of the best points of a scan that are apart), and for draws='auto': ess=400, rhat=1.01, max_draws=5000, batch=100
'pmcmc': draws='auto' (as for NUTS), warmup=1000, chains=4, max_draws=20000, batch=250, particles (of the particle filter: by default chosen by each chain in its warm-up; see ss.fitting.pmmh(); none for a deterministic sim, whose likelihood is exact), scan (each chain starts from one of the best points of a scan that are apart), parallel=True (to run the chains in parallel processes)
'if2': particles=1000 (of the particle filter), iterations=50 (passes of the filter in each chain), cooling=0.5 (the fraction of the first size of the random walk’s steps that is left after 50 passes, as pomp’s cooling.fraction.50), sd=0.05 (that first size, on the scale searched on), chains=4
'abc-smc': samples=200 (the parameter values kept in each generation, which are the draws of the posterior), generations=8, quantile=0.5 (of the distances, which sets the next tolerance), tolerance (the distance to stop at)
for a stochastic sim, how each point is scored (SciPy methods and Optuna; n_runs also for particle MCMC’s start): n_runs (runs per point: 4 for a stochastic compartmental sim, 1 for agents and for a deterministic one; every point is scored by the mean over the same seeds of each run’s score; for Optuna, giving it scores by runs in place of the particle filter), and reseed (whether each point takes fresh seeds: by default only for Optuna, whose search then averages over chance)
for the methods with gradients (‘l-bfgs-b’, ‘slsqp’, ‘nuts’): require_jax=True (whether the gradients have to be taken with the log-likelihood as one function in JAX; with False, something JAX does not cover yet is taken with Python in the loop instead, which is much slower)

Examples

import starsim as ss

beta = ss.par(0.3, bounds=[0.1, 1])
S, I, R = ss.states(3)
infection = S >> I @ beta*I
sir = ss.Disease(infection, I >> R @ ss.days(8), ss.outcome(0.5*infection, name='cases', likelihood=ss.poisson()), init=dict(I=0.01), name='sir')
pars = dict(total_pop=10_000, unit='day', stop=60, solver='binomial')
data = ss.Sim(sir, pars=dict(beta=0.5), rand_seed=2, **pars).run().to_synthetic_data('cases', resample='week')
sim = ss.Sim(sir, data=data, **pars)
fit = sim.fit(method=ss.fitting.method('optuna', n_trials=50, n_runs=4))
post = sim.fit(method=ss.fitting.method('pmcmc', draws=200, warmup=200, chains=2, particles=100))

Attributes

Name	Description
options	The settings for the method’s own function
scoring	The settings of how points are scored, for the Problem

minimize

internals.fitting.minimize(result, at=None)

The goal of sim.optimize(): a result to make as small as possible

Parameters

Name	Type	Description	Default
result	`str` / `func`	a result of a module (e.g. `'sir.R'`), or a function of the run sim that returns a number	required
at	`float` / date	the time of the result (default: the end); without it, a function of the sim decides	`None`

Example

best = sim.optimize(ss.fitting.minimize('sir.R'), spend=['tests', 'masks'])

Functions

Name	Description
abc_smc	Approximate Bayesian computation by sequential Monte Carlo (Beaumont et al. 2009): a posterior without a likelihood
compare_models	Compare models by how well they predict data they were not fitted to
diagnose	The R-hat (near 1 if the chains agree) and the effective sample size of each value, from draws (chains × draws × values)
fit_stages	Staged fitting: run each method in turn, each starting from the best fit of the one before; see `sim.fit()` with `method=` a list
loglik	The log-likelihood of a sim’s data (larger is better): the sum over the series, which are treated as independent
mismatch	The loss of a sim against its data: the mismatch between them (smaller is better)
nuts	Sample from a posterior with the No-U-Turn sampler (NUTS, from BlackJAX), with Stan’s warm-up of the step size and the mass matrix
pmmh	Sample from a posterior by particle marginal Metropolis-Hastings: a random-walk Metropolis sampler on a log-density that is estimated by a particle filter (or is exact)
prior_draws	Values of the model’s parameters drawn from their priors (`ss.par(prior=...)`, truncated to the bounds if it has them, as the samplers use it; or uniform over the bounds of one with bounds and no prior), for `ss.parallel(sim, priors=True)`
profile	The profile likelihood of one parameter: at each of its values, the others refitted by maximum likelihood (`sim.fit(pars=...)`), and the confidence interval it gives
sensitivity	Which parameters matter: Latin hypercube samples over the bounds of the parameters, with the partial rank correlation of each with a summary of each result
split_data	Split data into the part to fit to and the part to predict, for each fold

abc_smc

internals.fitting.abc_smc(
    distance,
    logprior,
    z0,
    samples=200,
    generations=8,
    quantile=0.5,
    tolerance=None,
    min_acceptance=0.02,
    rand_seed=None,
    verbose=True,
)

Approximate Bayesian computation by sequential Monte Carlo (Beaumont et al. 2009): a posterior without a likelihood

A population of samples (parameter values) is drawn from the prior, and in each generation a new one is drawn from the last (each sample perturbed by a Gaussian kernel twice as wide as the population) and kept only if its run is within the tolerance of the data; the tolerance is the median distance of the last population, so it shrinks as the population moves toward the data. Weights correct for the samples not being drawn from the prior.

Parameters

Name	Type	Description	Default
distance	`func`	the distance of a run from the data, given a point (on the sampler’s scale) and a seed; or that and something to keep with the point if it is accepted (e.g. the run)	required
logprior	`func`	the log of the prior density of a point, on the sampler’s scale	required
z0	`array`	a point on the sampler’s scale (only its length is used: the first generation is drawn uniformly on that scale and weighted by the prior)	required
samples	`int`	the size of the population	`200`
generations	`int`	the most generations; fewer if the tolerance or the acceptance rate says to stop	`8`
quantile	`float`	the fraction of the last population’s distances below the next tolerance	`0.5`
tolerance	`float`	stop once the tolerance is this small	`None`
min_acceptance	`float`	stop once a generation accepts fewer of its proposals than this	`0.02`
rand_seed	`int`	the random seed	`None`
verbose	`bool`	whether to say what is being done	`True`

Returns the population (samples × values, on the sampler’s scale), its normalized weights and distances, the tolerance and acceptance rate of each generation, and the number of runs; and if distance gives something to keep, what it gave with each sample of the population (or None for the first generation, whose are not kept).

compare_models

internals.fitting.compare_models(
    sims,
    data=None,
    holdout=0.25,
    folds=None,
    n_runs=None,
    **kwargs,
)

Compare models by how well they predict data they were not fitted to

Each model is fitted (with sim.fit()) to all but the last part of each data series, run forward over that part, and scored on it: by the default goodness of fit, or by the log predictive density if every series of every model has a likelihood; the same measure is used for every model. A model with more parameters usually fits better; the prediction error needs no separate penalty for them.

A stochastic model is run n_runs times at its best-fit values, and scored by the log of the mean over the runs of the likelihood (or by the mean loss), since its predictive probability is the average over the epidemics it could give, and a single run scores worse than that average. The runs start at the beginning and are not conditioned on the data fitted to (as a particle filter would condition them), so a run that drifts early is still scored on the held-out part, which is harsher on a stochastic model than on a deterministic one; and they are at the best-fit values, not averaged over a posterior.

Parameters

Name	Type	Description	Default
sims	`list` / `dict`	the sims to compare; labeled by their `label`, else by position (or by the keys of a dict)	required
data	`df` / `str`	the data, attached to every sim (default: each sim’s own data, which must be the same for all)	`None`
holdout	`float` / date	the fraction of each series to hold out (default: the last 25%), or the date (or time) from which to hold it out	`0.25`
folds	`int`	if given, repeat the comparison with this many cut points over the held-out part (fit to the first part, predict the next block), and report the mean and spread of the prediction error	`None`
n_runs	`int`	the number of runs of a stochastic model at its best-fit values that it is scored by (default: 100, as for `fit.project()`; a deterministic one is run once)	`None`
kwargs	`dict`	passed to `sim.fit()`, e.g. `weights=`

Returns

Name	Type	Description
		An `ss.fitting.ModelComparison`; `comp.fits` holds each model’s fit

Examples

comp = ss.fitting.compare_models([sim_sir, sim_seir], data=df)
comp.summarize()
comp.plot()

diagnose

internals.fitting.diagnose(draws)

The R-hat (near 1 if the chains agree) and the effective sample size of each value, from draws (chains × draws × values)

Both are rank-normalized and from split chains (Vehtari et al. 2021, as in Stan, ArviZ and BlackJAX, whose results these match): the R-hat is the larger of that of the draws and that of their distance from the median, and the effective sample size is of the bulk. Both are NaN for chains of fewer than 4 draws, too short to tell. In NumPy, since in JAX they took seconds to compile.

Parameters

Name	Type	Description	Default
draws	`array`	the draws, with shape (chains, draws, values)	required

fit_stages

internals.fitting.fit_stages(sim, stages, verbose=None, **kwargs)

Staged fitting: run each method in turn, each starting from the best fit of the one before; see sim.fit() with method= a list

Each stage is a method’s name, or ss.fitting.method() with that method’s settings. A stage that does not converge stops the sequence (with a warning), and the fit returned is that of the last stage that did, or the first stage’s best attempt if it was the one that failed.

Parameters

Name	Type	Description	Default
sim	Sim	the sim to fit	required
stages	`list`	the stages, each a method’s name, or an `ss.fitting.method()`	required
verbose	`bool`	whether to say what is being done (default: `ss.options.verbose`)	`None`
kwargs	`dict`	passed to `sim.fit()` for every stage

Examples

import starsim as ss

beta = ss.par(0.3, bounds=[0.1, 1])
S, I, R = ss.states(3)
infection = S >> I @ beta*I
sir = ss.Disease(infection, I >> R @ ss.days(8), init=dict(I=0.01), name='sir')
pars = dict(total_pop=10_000, unit='day', stop=60)
data = ss.Sim(sir, pars=dict(beta=0.5), **pars).run().to_synthetic_data('I')
fit = ss.Sim(sir, data=data, **pars).fit(method=[ss.fitting.method('optuna', n_trials=30), ss.fitting.method('nelder-mead', maxiter=100)]) # A search, then a polish
print(fit.pars)

loglik

internals.fitting.loglik(
    sim,
    data=None,
    likelihood=None,
    weights=None,
    by_series=False,
)

The log-likelihood of a sim’s data (larger is better): the sum over the series, which are treated as independent

Parameters

Name	Type	Description	Default
sim	`ss.Sim`	a sim that has been run	required
data	`df`	the data, if not the sim’s	`None`
likelihood	`ss.Dist`/dict	the likelihood of every series (e.g. `ss.poisson()`), or of each by name (e.g. `dict(cases=ss.poisson())`), if not declared on its outcome	`None`
weights	`dict`	a multiplier on each series’ log-likelihood (default 1)	`None`
by_series	`bool`	return a dict of the log-likelihood of each series instead of the total	`False`

mismatch

internals.fitting.mismatch(sim, data=None, weights=None, loss=None)

The loss of a sim against its data: the mismatch between them (smaller is better)

By default, for each data series, the normalized absolute error: the mean of |model − data| divided by the mean of the data; the series are then averaged. A series whose data are all zero is divided by the mean of the model instead, so it adds at most 1 (not raw counts, which would swamp the other series).

Parameters

Name	Type	Description	Default
sim	`ss.Sim`	a sim that has been run	required
data	`df`	the data, if not the sim’s	`None`
weights	`dict`	the relative weight of each data series (default: equal)	`None`
loss	`str` / `func`	‘nae’ (default), ‘rmse’, or a function of the sim returning a number	`None`

nuts

internals.fitting.nuts(
    logdensity,
    z0,
    draws='auto',
    warmup=300,
    chains=4,
    traced=False,
    forward=True,
    rand_seed=None,
    verbose=True,
    ess=ESS_TARGET,
    rhat=RHAT_TARGET,
    max_draws=5000,
    batch=100,
)

Sample from a posterior with the No-U-Turn sampler (NUTS, from BlackJAX), with Stan’s warm-up of the step size and the mass matrix

The log-density is given in one of two ways. By default, as a function of NumPy arrays that returns its value and gradient, so it can be anything (here, a run of the sim, differentiated by ss.gradients.gradient()); it is called back from JAX, and the chains are run one after another. With traced=True, as a JAX function (see ss.gradients.Density), so that the whole of the sampling is compiled, with the chains run together (each on cores of its own); its gradient is taken in forward mode (one direction per parameter), which for a few parameters is much cheaper than reverse mode through the steps of an ODE, or with forward=False in reverse mode (one pass whatever their number, for many parameters). The mass matrix is dense (it learns the correlations between values, which a fit’s often has) for up to DENSE_MAX values, and diagonal for more. With draws='auto' (the default), the chains are sampled in batches of batch draws until they agree, with enough effective draws (rhat, ess; see enough_draws()), or reach max_draws.

Parameters

Name	Type	Description	Default
logdensity	`func`	called with a point (an array, on an unbounded scale); returns the log of the posterior density there and its gradient (or, with `traced=True`, the log-density alone, as a JAX function)	required
z0	`array`	the point to start each chain near, or a point for each chain (chains × values)	required
draws	`int` / `str`	the draws to keep from each chain, or ‘auto’ to stop once they are enough	`'auto'`
warmup	`int`	the draws for the warm-up of each chain, which are not kept	`300`
chains	`int`	the number of chains	`4`
traced	`bool`	whether `logdensity` is a JAX function	`False`
forward	`bool`	for a JAX function, whether its gradient is taken in forward mode (else reverse mode, for which it has to be differentiable)	`True`
rand_seed	`int`	the random seed	`None`
verbose	`bool`	whether to say what is being done	`True`
ess	`float`	with `draws='auto'`, the effective draws of each value to stop at	`ESS_TARGET`
rhat	`float`	with `draws='auto'`, the R-hat to stop at	`RHAT_TARGET`
max_draws	`int`	with `draws='auto'`, the most draws to take from each chain	`5000`
batch	`int`	with `draws='auto'`, the draws each chain takes between checks	`100`

Returns

Name	Type	Description
		The draws (chains × draws × values), the log-density at each (chains × draws), and whether each followed a divergent trajectory (chains × draws)

pmmh

internals.fitting.pmmh(
    logpost,
    z0,
    draws='auto',
    warmup=1000,
    chains=4,
    particles=None,
    rand_seed=None,
    parallel=True,
    verbose=True,
    ess=ESS_TARGET,
    rhat=RHAT_TARGET,
    max_draws=20000,
    batch=250,
)

Sample from a posterior by particle marginal Metropolis-Hastings: a random-walk Metropolis sampler on a log-density that is estimated by a particle filter (or is exact)

Since the estimate of the likelihood is unbiased, the chains sample from the exact posterior (Andrieu, Doucet & Holenstein 2010), however noisy it is; more noise only makes them move more slowly. In the warm-up, the proposal’s covariance is learned from the chain so far (Haario et al. 2001), and its scale is tuned toward accepting 15% of proposals (about the best when the estimate is noisy; Sherlock et al. 2015; 23.4% for an exact log-density). In the first half of the warm-up, the estimate where the chain is is made again at each step, so that a chain far from the posterior, where the estimate is noisiest, cannot be held by a lucky one. With particles=0, the log-density is exact (e.g. of an ODE), so it is worked out once at each point, and this is an adaptive random-walk Metropolis sampler. With draws='auto' (the default), the chains are taken further in rounds of batch draws until they agree, with enough effective draws (rhat, ess; see enough_draws()), or reach max_draws.

Parameters

Name	Type	Description	Default
logpost	`func`	called with a point (an array, on an unbounded scale), a random number generator (`np.random.Generator`) to draw any seeds from, and the number of particles; returns an estimate of the log of the posterior density there, or that and something to keep with each draw at the point (e.g. a run of the sim drawn from the filter)	required
z0	`array`	the point to start each chain from, or a point for each chain (chains × values)	required
draws	`int` / `str`	the draws to keep from each chain, or ‘auto’ to stop once they are enough	`'auto'`
warmup	`int`	the draws for the warm-up of each chain, which are not kept	`1000`
chains	`int`	the number of chains	`4`
particles	`int`	the number of particles; by default, each chain uses 100 for the first half of its warm-up, and then enough that the variance of the estimate at the middle of where it has been lately is about 1 (Doucet et al. 2015); 0 for a log-density that is exact	`None`
rand_seed	`int`	the random seed	`None`
parallel	`bool`	whether to run the chains in parallel processes	`True`
verbose	`bool`	whether to say what is being done	`True`
ess	`float`	with `draws='auto'`, the effective draws of each value to stop at	`ESS_TARGET`
rhat	`float`	with `draws='auto'`, the R-hat to stop at	`RHAT_TARGET`
max_draws	`int`	with `draws='auto'`, the most draws to take from each chain	`20000`
batch	`int`	with `draws='auto'`, the draws each chain takes between checks	`250`

Returns

Name	Type	Description
		The draws (chains × draws × values), the estimate of the log-density at each (chains × draws), the fraction of proposals accepted after the warm-up in each chain, and the number of particles that each chain used;
		and if `logpost` gives something to keep, what it gave at each draw (a list for each chain, in which a draw that stayed where the chain was holds the same object as the draw before)

prior_draws

internals.fitting.prior_draws(sim, n, rand_seed=None)

Values of the model’s parameters drawn from their priors (ss.par(prior=...), truncated to the bounds if it has them, as the samplers use it; or uniform over the bounds of one with bounds and no prior), for ss.parallel(sim, priors=True)

The draws are a Latin hypercube, as in sim.sensitivity(): each parameter’s range is cut into n slices of equal probability with one draw in each, which covers it evenly in few runs.

Returns a dict with a list of n values for each parameter that has a prior or bounds; a parameter with a value for each level (ss.par(ss.by(age=[...]), prior=...)) has a draw for each level, as “beta[a]”, “beta[b]”, …, as fit() fits it.

profile

internals.fitting.profile(sim, par, values=None, level=0.95, **kwargs)

The profile likelihood of one parameter: at each of its values, the others refitted by maximum likelihood (sim.fit(pars=...)), and the confidence interval it gives

The interval is the values whose negative log-likelihood is within half the chi-squared quantile with one degree of freedom (1.92 for 95%) of the lowest. Unlike the spread of a point estimate, it also shows whether the data determine the parameter at all: a profile that stays flat to a bound does not end within it.

For a deterministic sim (an ODE) whose data series all have a likelihood, and parameters with no priors (so each refit is the maximum likelihood); a stochastic sim’s estimates of the likelihood are noisy, and need a smoothed profile, which is not yet supported.

Parameters

Name	Type	Description	Default
sim	`ss.Sim`	the sim, with data, which is not changed	required
par	`str`	the name of the parameter to profile, one with a single value	required
values	`list`	its values (default: the best fit of every parameter, 10 points spread over the bounds of this one, on the scale it is searched on, then 10 more at a time over the interval and the points next to it, until 5 are within it, at most 4 times)	`None`
level	`float`	the confidence level of the interval	`0.95`
kwargs	`dict`	passed to `sim.fit()` for each fit, e.g. `likelihood=` or `pars=` (the parameters to refit; default: all that have bounds)

Returns an objdict with table (a dataframe with a row for each value: the value, the refitted values of the others, the negative log-likelihood nll, and whether the value is in the interval), best (the value with the lowest), interval (its low and high ends, where the profile crosses the threshold, interpolated between the values; an end that the values do not reach is the outermost of them, with a warning), and threshold.

Examples

prof = ss.fitting.profile(sim, 'beta', likelihood=dict(cases=ss.poisson()))
prof.interval # The 95% confidence interval of beta
prof.table.plot('beta', 'nll')

sensitivity

internals.fitting.sensitivity(
    sim,
    pars=None,
    results=None,
    n=100,
    rand_seed=None,
    verbose=None,
)

Which parameters matter: Latin hypercube samples over the bounds of the parameters, with the partial rank correlation of each with a summary of each result

Needs no data. A result is summarized by its peak and its final value (a flow by its peak and its total). A partial rank correlation near 1 or -1 means the result goes up or down with the parameter, with the others held fixed; near 0, that it does not matter, so it cannot be fitted either.

Parameters

Name	Type	Description	Default
sim	`ss.Sim`	the sim, which is not changed	required
pars	`list`	the names of the parameters (default: all that have bounds)	`None`
results	`list`	the results, e.g. `['I']` (default: every count of a state, and every flow)	`None`
n	`int`	the number of samples	`100`
rand_seed	`int`	the seed of the samples	`None`
verbose	`bool`	whether to print the table	`None`

Returns a dataframe of the partial rank correlations, with a row for each parameter and a column for each summary of a result.

split_data

internals.fitting.split_data(data, holdout=0.25, folds=None)

Split data into the part to fit to and the part to predict, for each fold

Each series is split on its own times. With a fraction, the last holdout of each series’ points is held out; with folds=k, that part is cut into k blocks, and fold j fits to everything before block j and predicts block j. With a date or time, everything from it on is held out.

Returns a list of (fitting data, held-out data) pairs, one per fold.

Parameters

Name	Type	Description	Default
data	`dataframe`	the data	required
holdout	`float` / date	the fraction of each series to hold out, or the date or time from which to hold out everything	`0.25`
folds	`int`	the number of folds	`None`