Skip to content

Commit

Permalink
sir model update
Browse files Browse the repository at this point in the history
  • Loading branch information
stevencarlislewalker committed Dec 13, 2023
1 parent 9f603bb commit afaf0fb
Show file tree
Hide file tree
Showing 2 changed files with 22 additions and 20 deletions.
19 changes: 15 additions & 4 deletions inst/starter_models/sir/example.R
Original file line number Diff line number Diff line change
@@ -1,11 +1,22 @@
source("inst/starter_models/sir/model.R")

## -------------------------
## parameterize model
## -------------------------

tmb_simulator$update$transformations(Log("beta"))
tmb_simulator$replace$params(log(init_mats$get("beta")), "log_beta")
tmb_simulator ## note the new expression before the simulation loop

## -------------------------
## simulate fake data
## -------------------------

time_steps = 100L
true_beta = 0.4
observed_data = tmb_simulator$report(true_beta)

tmb_simulator$replace$time_steps(100L)
observed_data = tmb_simulator$report(log(true_beta))
observed_data$value = rpois(100, observed_data$value)

if (interactive()) {
Expand All @@ -28,13 +39,13 @@ tmb_simulator$update$matrices(
## -------------------------

if (interactive()) {
betas = seq(from = 0.01, to = 1, length = 100)
log_betas = seq(from = log(0.1), to = log(1), length = 100)
ll = vapply(
betas
log_betas
, tmb_simulator$objective
, numeric(1L)
)
plot(betas, ll, type = "l")
plot(exp(log_betas), ll, type = "l")
abline(v = true_beta)
}

Expand Down
23 changes: 7 additions & 16 deletions inst/starter_models/sir/model.R
Original file line number Diff line number Diff line change
Expand Up @@ -11,23 +11,20 @@ computations = list(
## absolute flow rates (per time only)
flow_rates = list(
infection ~ S * I * beta / N
, recovered ~ gamma * I
, recovery ~ gamma * I
)

## state updates
state_updates = list(
S ~ S - infection
, I ~ I + infection - recovered
, R ~ R + recovered
, I ~ I + infection - recovery
, R ~ R + recovery
)

model_evaluation = list(
log_likelihood ~ dpois(I_obs, rbind_time(I, I_obs_times))
)

## data model comparison
obj_fn = ~ -sum(log_likelihood)

## simple unstructured scalar expression
expr_list = ExprList(
before = computations
Expand All @@ -44,27 +41,21 @@ init_mats = MatsList(
, gamma = 0.1
, N = 100
, infection = empty_matrix
, recovered = empty_matrix
, recovery = empty_matrix
, log_likelihood = empty_matrix
, I_obs = empty_matrix
, I_obs_times = empty_matrix
, .mats_to_save = "I"
, .mats_to_return = "I"
)

obj_fn = ObjectiveFunction(~ -sum(log_likelihood))

## simulator
###################################################

tmb_simulator = TMBModel(
init_mats = init_mats
, expr_list = expr_list
, time_steps = Time(100L)
, params = OptParamsList(
init_mats$get("beta")
, par_id = c(0L)
, mat = c("beta")
, row_id = c(0L)
, col_id = c(0L)
)
, obj_fn = ObjectiveFunction(obj_fn)
, obj_fn = obj_fn
)$simulator()

0 comments on commit afaf0fb

Please sign in to comment.