Skip to main content

SimInf

🧠0% 🤖100%

This page is for people who use SimInf in R: stochastic compartmental models in each of thousands of herds or farms, linked by a register of births, deaths, and movements between holdings, written with mparse() or SimInf’s predefined models, and fitted with pfilter(), pmcmc(), and abc(). It shows a SimInf model next to the same model in Starsim, how their results agree, how long each takes, and what works and what does not yet, including what fails.

The same model in both

SimInf’s events_SIR() example: an SIR in each of 1,600 cattle herds over four years, linked only by 466,692 scheduled events (animals entering, leaving, and moving between herds):

#| eval: false
set.seed(123)
set_num_threads(1)

u0 <- u0_SIR()
u0$I[1] <- 1
model <- SIR(
    u0     = u0,
    tspan  = seq(from = 1, to = 4*365, by = 1),
    events = events_SIR(),
    beta   = 0.16,
    gamma  = 0.01
)

result <- run(model)
plot(result)

In Starsim, SimInf’s events become a movement register, a table with columns date, from, to, n, and states, in place of event codes and the select and shift matrices. Here u0.csv and events.csv.gz are SimInf’s u0_SIR() (with the one infected animal) and events_SIR(), saved from R. Each herd is a level of a stratification, and the model runs by Gillespie’s method in each herd, as SimInf does:

import pandas as pd
import starsim as ss

u0 = pd.read_csv('u0.csv')
events = pd.read_csv('events.csv.gz')

# SimInf's events as a register: entering animals are susceptible, and the others are chosen from every state
register = pd.DataFrame({
    'date':   events.time,
    'from':   events.node.where(events.event != 'enter'),
    'to':     events.node.where(events.event == 'enter', events.dest.where(events.event == 'extTrans')),
    'n':      events.n,
    'states': (events.event == 'enter').map({True:'S', False:None}),
})

S, I, R = ss.states('S I R')
sir = ss.Disease(
    infection = ss.Transition(S, I, rate=ss.perday(0.16) * I),
    recovery  = ss.Transition(I, R, rate=ss.perday(0.01)),
    init_n    = dict(I=ss.by(herd=u0.I.values)),
)
size = u0[['S', 'I', 'R']].sum(axis=1)
sir.stratify(herd=ss.levels(dict(zip(list(u0.node.astype(str)), size/size.sum())), kind='space', mobility=ss.geo.register(register)))
sim = ss.Sim(sir, solver='gillespie', unit='day', start=0, stop=1460, total_pop=size.sum())
sim.run()
sim.plot()

This cell is not run when the docs are built, since it needs SimInf’s data. Seven of its lines turn SimInf’s events into a register and four read the data, which SimInf gets from the package; that is why it is longer than the original.

How the results agree

Model Agreement
events_SIR() (above) Agree: the mean of 20 runs is within about two standard errors of SimInf’s mean of 200 at every day checked, for the totals, the number of herds infected, and the first herd (susceptible on day 100: 199,916 ± 11 against 199,906)
SISe3 (three age groups, environmental infectious pressure phi with seasonal decay, 783,773 events) Agree, within a herd, between herds, and in the demography: every series within about one standard error of SimInf’s mean on every day checked (on day 1,457, 184 ± 16 herds infected against 166)
A rate for each herd (mparse() with ldata) Agree: infected on day 1,460, 7,267 ± 688 (20 runs) against 7,613 ± 363 (50 runs); herds infected 1,111 ± 105 against 1,178 ± 56

In SISe3, what SimInf does in C (the environmental pressure, its decay, and its daily update) is declared: phi is a state that is not people, and shedding and decay are two transitions.

Timings

Lines are lines of model code (not counting blank lines, comments, imports, or plotting, and counting a line over 100 characters twice); times are for one run, the first in a new process and a rerun of the same model. Both were measured as described on the comparison page. SimInf 11.1.0’s times were recorded on the same machine. Starsim fails the north star’s test of lines on both models, and passes on time, being less than 1 s slower.

Model Lines: SimInf Lines: Starsim One run: SimInf One run: Starsim
events_SIR(), 1,600 herds 12 21 first 0.17 s, rerun 0.078 s first 0.27 s, rerun 0.14 s Fail (lines)
SISe3, 1,600 herds 27 39 first 0.31 s, rerun 0.094 s first 0.92 s, rerun 0.55 s Fail (lines)
events_SIR() copied ten times: 16,000 herds and 4.7 million events 1.9 s to run (2.0 s on 8 threads), 2.1 s to build 3.1 s to run, 2.0 s to build Pass

The last row was measured on 18 October 2026 on a machine shared with other jobs, and is the switching demo’s test_siminf_scale(), with the results kept for every herd on every day, as SimInf keeps them. Starsim’s Gillespie solver steps only the herds that can have an event (a herd with nobody infected is skipped until the register moves animals into or out of it), and is always exact, however many events there are.

What works and what does not yet

These are from Starsim’s switching demo for SimInf users (tests/switching/test_siminf.py in this repo), which checks each task against SimInf’s numbers, and from the review it was built from. “Fails” means the demo’s check does not pass today.

Task Today
Herds linked by a register of births, deaths, and movements Works: on SimInf’s 1,600 herds for 100 days, the susceptible total and the number of herds infected are inside SimInf’s 5% to 95% range
Selective movements by state (states=), and movements that change state (to_states=) Works: columns of the register
A rate for each herd, written ss.perday(ss.by(herd=b)) Works: b is a list or a NumPy array, one value per herd
A model with a continuous state per herd (SISe’s phi) Works: declared, with no C
Fit by particle MCMC (pmcmc()) Works: sim.fit(method='pmcmc'), for a stochastic model with a likelihood, runs four chains, and reports R-hat; on SimInf’s ?pmcmc example the 95% intervals (beta 0.123 to 0.186, gamma 0.066 to 0.089) cover the true values (0.16 and 0.077)
The speed of the particle filter Works: 0.57 ms per particle run on one thread, the same as SimInf’s 0.57 ms (the demo asks for within 3 times), and 0.23 ms on four threads
Fit by ABC-SMC (abc()) Works: method='abc-smc', below
Speed on a national register Works: see the timings above
Build and run a model from R Works: through reticulate, with ss$Transition(), states named with ss$states(), and name=
The same model as animals Works: sim.new(solver='abm', n_agents=...) runs one agent per animal (196,169 animals in 108 s for events_SIR())

SimInf’s ?pmcmc example in miniature, fitted by ABC-SMC: one herd of 1,000 with an SIR at beta = 0.16 and gamma = 0.077, observed through the number infected every five days:

import numpy as np
import pandas as pd
import starsim as ss

data = pd.DataFrame(dict(t=np.arange(5, 101, 5), I=[9, 9, 22, 25, 25, 56, 56, 83, 93, 113, 131, 152, 171, 162, 118, 128, 128, 107, 71, 76]))
beta  = ss.par(ss.perday(0.2), bounds=[0.01, 1])
gamma = ss.par(ss.perday(0.1), bounds=[0.01, 1])
S, I, R = ss.states('S I R')
sir = ss.Disease(infection=ss.Transition(S, I, rate=beta * I), recovery=ss.Transition(I, R, rate=gamma), init_n=dict(I=5))
sim = ss.Sim(sir, total_pop=1000, solver='gillespie', unit='day', start=1, stop=100, data=data, verbose=0)

fit = sim.fit(method=ss.fitting.method('abc-smc', samples=100, generations=7), rand_seed=1, verbose=0)
fit.table # The posterior median of each parameter and its 95% interval (lower, upper)
par initial best low high median lower upper ess
0 beta 0.2 0.157636 0.01 1.0 0.179391 0.062257 0.291824 75.852305
1 gamma 0.1 0.078069 0.01 1.0 0.085115 0.024455 0.195542 75.852305

Both intervals cover the true values. They are wider than those of particle MCMC on the same data (in the table above), but that is not a rule: ABC-SMC compares the runs with the data by a distance, with no likelihood, and its result approximates the posterior as closely as its final tolerance allows, so its width depends on the number of generations as well as on the data. With a likelihood, particle MCMC gives the posterior itself.

Bringing your model in

ss.import_from('model.R', format='siminf') reads a SimInf model: transitions written with mparse(), or one of SimInf’s predefined models (SIR(), SIS(), SEIR(), SISe(), SISe3()). The R is read, not run, so what the script computes (such as u0_SIR() or a data file) can be passed in, e.g. ss.import_from('model.R', format='siminf', u0='u0.csv', events='events.csv'). The nodes become a stratification of places, the events a movement register, and the continuous state of the SISe models a state that is not people. C code (pts_fun) is listed as not translated, and sim.imported.prompt is a prompt for an AI assistant with Starsim’s convert-model skill to finish the job. sim.describe() shows how each part was read.

See Import/export for importing in general: checking an import, keeping it as code, and finishing one that is incomplete.

Further reading