Grey-seal integrated population model
This is the capstone case study: a full integrated population model (IPM) for the Baltic grey seal, ported from the n-kall/sealIPM Stan program. An IPM fits one latent population process to several independent data streams at once — here eight of them, from aerial pup counts to hunting-bag totals to reproductive-tract signs — so the shared demographic parameters are informed by everything observed about the population.
It exercises nearly the whole StanBlocks surface in one model: a year-recursive state process (a scan), a numerical ODE solved inside that scan, a custom Dirichlet-multinomial fate allocation, six user-defined distribution families, two observation-stream submodels, and named-tuple state carriers. Because the pieces are many and interdependent, the model is assembled from a small library of cards with compile_slic_bundle rather than one inline @slic block — 31 @deffun helper cards, two anonymous observation-stream submodels, and one parent body. Everything below is evaluated at documentation-build time, so the displayed Julia is the exact source that produced the ~52 KB Stan program beside it. The arrays are a tiny build fixture (3 age classes × 2 sexes × 3 years), not real survey data.
Companion port in BayesianRegressionModels.jl
The same grey-seal IPM is also ported in BRM, onto its StanBlocks backend: Grey-seal IPM. Both port the identical n-kall/sealIPM source verbatim through compile_slic_bundle.
The shape of the model
The parent body is short and readable — it declares the demographic and observation parameters, computes a handful of transformed quantities, runs the state process once, and then attaches each observation stream as a single line.
The state process is a scan.
run_state_processwalks the years forward: each year updates the density-dependent birth rate, the age/sex population composition, the Sweden/Finland hunting pressure (whose expected catch is anode_rk45_tolintegral of a hunting-hazard ODE,dH_dt), and a stochastic split of every demographic class into survived / bycaught / hunted-SE / hunted-FI viamultinomial_allocation. A recurrence like this cannot live in@slic(no control flow) — it is a@deffun, called loop-free from the body.It returns a named tuple of state carriers.
run_state_processreturns(; birth_rate, pregnancy_rate, population_total, hunted_sweden, …, reproductive_probs), and the body reads them by name —state.population_total,state.bycatch_expected. (StanBlocks named tuples are authoring sugar that lower to Stan's positional tuple access; see the wastewater study's tuple note.)Eight observation streams, each one removable line. Every stream is a Form-A observation (
data ~ family(...)ordata ~ submodel(...)): the real observed vector is on the left, so deleting the line drops that stream (and a submodel node drops its own parameters too). Two streams that carry their own parameters — the aerial count's over-dispersion, the bycatch composition's selectivity bias — are submodels; the other six call custom families directly.Six custom distribution families.
aerial_count(negative-binomial),harvest_bags(normal with a shared CV),hunting_comp/bycatch_comp/reproductive_signs(multinomial compositions), andpregnancy(binomial) are each defined as a@deffun_lpmf/_lpmfs/_rngtriad (the density card carries@lhs @lpxf), so the model gets posterior-predictive draws and pointwise log-likelihoods for every stream for free.
julia
# ---- STATE-PROCESS parameters (base — the scan needs them) ----
phi_a_sc ~ uniform(0.0, 1.0)
phi_sc ~ uniform(0.0, 1.0)
survival_shape ~ uniform(0.0, 1.0)
male_pup_survival_offset ~ cauchy(0.0, 1.0)
male_adult_survival_offset ~ cauchy(0.0, 1.0)
carrying_capacity ~ lognormal(11.0, 1.0)
max_baseline_birth_rate ~ uniform(0.0, 1.0)
min_baseline_birth_rate_sc ~ uniform(0.0, 1.0)
herring_intercept_scaled ~ normal(0.0, 1.0)
herring_slope ~ normal(0.0, 1.0)
herring_weight ~ uniform(0.0, 1.0)
hunting_selectivity_sweden :: vector[n_demo] ~ normal(0.0, 1.0)
hunting_selectivity_finland :: vector[n_demo] ~ normal(0.0, 1.0)
hunting_effort_sd_sweden ~ cauchy(0.0, 1.0; lower=0.0)
hunting_effort_sd_finland ~ cauchy(0.0, 1.0; lower=0.0)
population_init_size ~ lognormal(11.0, 1.0)
# harvest-bag CV: an OBSERVATION param, but shared by SW+FI, so it can't live
# in a single Form-A stream (two obs would share it) — it stays in the base.
harvest_bag_cv ~ lognormal(0.0, 1.0; lower=0.0)
epsilon_birth :: vector[n_state_years] ~ std_normal()
epsilon_sex :: vector[n_state_years] ~ std_normal()
epsilon_h_sw :: vector[n_state_years] ~ std_normal()
epsilon_h_fi :: vector[n_state_years] ~ std_normal()
transition_noise_vec :: vector[3 * n_demo * n_state_years] ~ std_normal()
# reproductive-signs reporting params (feed pi_s/pi_c into the scan — FLAG h)
report_ca_mean ~ uniform(0.0, 1.0)
report_placental_mean ~ uniform(0.0, 1.0)
prob_of_ca ~ uniform(0.0, 1.0)
report_placental_sd ~ normal(0.0, 0.1; lower=0.0)
report_ca_sd ~ normal(0.0, 0.1; lower=0.0)
epsilon_ca :: vector[n_state_years] ~ std_normal(; lower=0.0)
epsilon_placental :: vector[n_state_years] ~ std_normal(; lower=0.0)
# ---- transformed quantities feeding the scan ----
phi_a = phi_a_sc
phi_pup = phi_sc * phi_a
mu_m = mortality_rates(phi_pup, phi_a, survival_shape, n_age,
male_pup_survival_offset, male_adult_survival_offset)
S_diag = exp(-mu_m)
aging = create_aging_matrix(n_demo, n_age)
baseline_bbr = compute_baseline_birth_rate(
max_baseline_birth_rate * min_baseline_birth_rate_sc, max_baseline_birth_rate,
herring_intercept_scaled, herring_slope, herring_weight,
herring_index_1, herring_index_2)
dd_scaled = birth_rate_at_carrying_capacity(phi_a, mu_m, n_age)
dd_intercept = compute_density_dependence_intercept(max_baseline_birth_rate, dd_scaled)
dd_slope = -log(carrying_capacity)
pop_first = initialize_population(population_init, population_init_size,
population_burn_in, baseline_bbr[1], aging, S_diag, n_age)
pi_s = report_placental_mean * exp(-epsilon_placental * report_placental_sd)
pi_c = report_ca_mean * exp(-epsilon_ca * report_ca_sd)
transition_noise_raw = to_matrix(
transition_noise_vec, 3 * n_demo, n_state_years)
# ---- THE state process: computed ONCE (state carriers) ----
state = run_state_process(
n_state_years, n_age, pop_first, baseline_bbr[1], sum(pop_first),
baseline_bbr, dd_intercept, dd_slope, aging, S_diag, mu_m,
hunting_selectivity_sweden, hunting_selectivity_finland,
hunting_quota_sweden, hunting_quota_finland,
hunting_effort_sd_sweden, hunting_effort_sd_finland,
epsilon_h_sw, epsilon_h_fi, t_mate_to_preg, t_birth_to_end_hunt,
epsilon_birth, epsilon_sex, transition_noise_raw,
pi_s, pi_c, prob_of_ca, ode_init_state, ode_times)
# =========================================================================
# OBSERVATION NODES — Form A: the REAL observed data is on the LHS. DELETE a
# line to deactivate that stream (a submodel node also drops its own params).
# The 2 submodels take EVERY name they read as a KWARG (no scope-flow —
# StanBlocks snag data-submodel-li-75dc835a); the 6 direct family calls read
# the model's own `state`/data/params in scope. Embed only via `~`.
# =========================================================================
obs_aerial_count ~ aerial_stream(;
aerial_year, population_total = state.population_total)
obs_bycatch_comp ~ bycatch_stream(;
bycatch_comp_year, bycatch_comp_sample_size, n_demo,
bycatch_expected = state.bycatch_expected)
obs_hunting_bag_sweden ~ harvest_bags(hunting_bag_year_sweden,
state.hunting_bag_total_sweden, harvest_bag_cv)
obs_hunting_bag_finland ~ harvest_bags(hunting_bag_year_finland,
state.hunting_bag_total_finland, harvest_bag_cv)
obs_hunting_comp_sweden ~ hunting_comp(hunting_comp_year_sweden,
state.hunted_sweden, state.hunting_bag_total_sweden,
hunting_comp_sample_size_sweden)
obs_hunting_comp_finland ~ hunting_comp(hunting_comp_year_finland,
state.hunted_finland, state.hunting_bag_total_finland,
hunting_comp_sample_size_finland)
obs_pregnancy_count ~ pregnancy(pregnancy_count_year,
pregnancy_sample_size, state.pregnancy_rate)
obs_reproductive_signs_finland ~ reproductive_signs(reproductive_signs_year,
state.reproductive_probs, reproductive_signs_sample_size)julia
run_state_process(n_state_years::int, n_age::int,
pop_first::vector[n_demo], birth_rate_first::real, pop_total_first::real,
baseline_bbr::vector[Tb], dd_intercept::real, dd_slope::real,
aging::matrix[n_demo, n_demo], S_diag::vector[n_demo], mu_m::vector[n_demo],
hs_sw::vector[n_demo], hs_fi::vector[n_demo],
hq_sw::int[n_state_years], hq_fi::int[n_state_years],
he_sd_sw::real, he_sd_fi::real,
eps_h_sw::vector[n_state_years], eps_h_fi::vector[n_state_years],
t_mate_to_preg::real, t_birth_to_end_hunt::real,
eps_birth::vector[n_state_years], eps_sex::vector[n_state_years],
transition_noise_raw::matrix[Tn, n_state_years],
pi_s::vector[n_state_years], pi_c::vector[n_state_years], prob_of_ca::real,
ode_init_state::vector[1], ode_times::vector[1]) = begin
n_demo = 2 * n_age
ode_ts = to_array_1d(ode_times)
birth_rate::vector[n_state_years]
pregnancy_rate::vector[n_state_years]
population_total::vector[n_state_years]
hunted_sweden::matrix[n_demo, n_state_years]
hunted_finland::matrix[n_demo, n_state_years]
bycatch_expected::matrix[n_demo, n_state_years]
hunting_bag_total_sweden::vector[n_state_years]
hunting_bag_total_finland::vector[n_state_years]
reproductive_probs::matrix[4, n_state_years]
population_comp::matrix[n_demo, n_state_years]
survivors::matrix[n_demo, n_state_years]
for year in 1:n_state_years
if year == 1
birth_rate[year] = birth_rate_first
population_comp[:, year] = pop_first
population_total[year] = pop_total_first
else
birth_rate[year] = update_birth_rate(
baseline_bbr[year], dd_intercept, dd_slope, population_total[year - 1])
population_comp[:, year] = update_population_from_survivors(
survivors[:, year - 1], aging, birth_rate[year],
eps_birth[year], eps_sex[year], n_age)
population_total[year] = sum(population_comp[:, year])
end
pregnancy_rate[year] = update_pregnancy_rate(
baseline_bbr[year + 1], dd_intercept, dd_slope,
population_total[year], t_mate_to_preg)
hp_sw::vector[n_demo]
hp_fi::vector[n_demo]
log_N = log(population_comp[:, year])
log_denom_sw = log_sum_exp(hs_sw + log_N)
log_denom_fi = log_sum_exp(hs_fi + log_N)
if hq_sw[year] == 0
hp_sw = rep_vector(0.0, n_demo)
else
hp_sw = exp(hs_sw + log(hq_sw[year]) + log(2.0)
- 2.0 * log(t_birth_to_end_hunt)
- eps_h_sw[year] * he_sd_sw - log_denom_sw)
end
if hq_fi[year] == 0
hp_fi = rep_vector(0.0, n_demo)
else
hp_fi = exp(hs_fi + log(hq_fi[year]) + log(2.0)
- 2.0 * log(t_birth_to_end_hunt)
- eps_h_fi[year] * he_sd_fi - log_denom_fi)
end
exp_hunted_sw::vector[n_demo]
exp_hunted_fi::vector[n_demo]
for demo in 1:n_demo
# Reference package defaults are literal here because Stan requires
# solver controls to be data-only and @deffun has no such qualifier yet.
sol_sw = ode_rk45_tol(dH_dt, ode_init_state, 0.0, ode_ts, 1.0e-6, 1.0e-6, 1000,
population_comp[demo, year], t_birth_to_end_hunt,
hp_sw[demo], hp_fi[demo], mu_m[demo])
exp_hunted_sw[demo] = sol_sw[1][1]
sol_fi = ode_rk45_tol(dH_dt, ode_init_state, 0.0, ode_ts, 1.0e-6, 1.0e-6, 1000,
population_comp[demo, year], t_birth_to_end_hunt,
hp_fi[demo], hp_sw[demo], mu_m[demo])
exp_hunted_fi[demo] = sol_fi[1][1]
end
transition_matrix = create_transition_matrix(
exp_hunted_sw, exp_hunted_fi, population_comp[:, year],
hp_sw, hp_fi, t_birth_to_end_hunt, S_diag)
expected_fate = to_matrix(transition_matrix * population_comp[:, year], n_demo, 4)
noise_year = to_matrix(transition_noise_raw[:, year], n_demo, 3)
realized_fate::matrix[n_demo, 4]
for demo in 1:n_demo
realized_fate[demo, :] = multinomial_allocation(
expected_fate[demo, :], noise_year[demo, :], population_comp[demo, year])
end
survivors[:, year] = realized_fate[:, 1]
bycatch_expected[:, year] = realized_fate[:, 2]
hunted_sweden[:, year] = realized_fate[:, 3]
hunted_finland[:, year] = realized_fate[:, 4]
hunting_bag_total_sweden[year] = sum(hunted_sweden[:, year])
hunting_bag_total_finland[year] = sum(hunted_finland[:, year])
reproductive_probs[2, year] = birth_rate[year] * pi_s[year] * (1.0 - pi_c[year])
reproductive_probs[3, year] = birth_rate[year] * (1.0 - pi_s[year]) * pi_c[year] +
(1.0 - birth_rate[year]) * prob_of_ca * pi_c[year]
reproductive_probs[4, year] = birth_rate[year] * pi_s[year] * pi_c[year]
reproductive_probs[1, year] = 1.0 - sum(reproductive_probs[2:4, year])
end
(; birth_rate, pregnancy_rate, population_total,
hunted_sweden, hunted_finland, bycatch_expected,
hunting_bag_total_sweden, hunting_bag_total_finland, reproductive_probs)
endjulia
dH_dt(tau::real, H::vector[ny], n0::real, k::real,
E_1::real, E_2::real, mu::real)::vector[ny] = begin
surv = exp(-(E_1 + E_2) * (k * tau - tau * tau / 2) - mu * tau)
rep_vector(n0 * E_1 * surv * (k - tau), 1)
endjulia
multinomial_allocation(eta_row::row_vector[4], u_row::row_vector[3], N::real)::row_vector[4] = begin
eta = eta_row'
eta_adj = eta * (1.0 + 1.0 / min(eta))
mean_logratio = (digamma(eta_adj[2:4]) - digamma(eta_adj[1]))'
Sigma = rep_matrix(trigamma(eta_adj[1]), 3, 3) + diag_matrix(trigamma(eta_adj[2:4]))
L = cholesky_decompose(Sigma)
logits = append_col(rep_row_vector(0.0, 1), mean_logratio + u_row * L')
allocation = softmax(logits')'
allocation * N
endjulia
# @slic markers: stanonly, lhs, lpxf
aerial_count_lpmf(obs::int[n_obs], year::int[n_obs],
population_total::vector[T], mu::real, phi::real)::real = begin
neg_binomial_2_lpmf(obs, mu * population_total[year], phi)
endstan
functions {
vector mortality_rates(
real phi_pup,
real phi_adult,
real c,
int n_age,
real male_pup_offset,
real male_adult_offset
) {
int n_demo = (2 * n_age);
vector[n_demo] mu_m;
real mu_pup_f = (-log(phi_pup));
real mu_ad_f = (-log(phi_adult));
real mu_pup_m = exp((log(mu_pup_f) + male_pup_offset));
real mu_ad_m = exp((log(mu_ad_f) + male_adult_offset));
mu_m[1] = mu_pup_f;
mu_m[n_age] = mu_ad_f;
for(j in 2:(n_age - 1)) {
real w = exp((c * log(((j - 1.0) / (n_age - 1.0)))));
mu_m[j] = exp((log(mu_pup_f) + (w * (log(mu_ad_f) - log(mu_pup_f)))));
}
mu_m[(n_age + 1)] = mu_pup_m;
mu_m[n_demo] = mu_ad_m;
for(j in 2:(n_age - 1)) {
real w_male = exp((c * log(((j - 1.0) / (n_age - 1.0)))));
mu_m[(n_age + j)] = exp((log(mu_pup_m) + (w_male * (log(mu_ad_m) - log(mu_pup_m)))));
}
return mu_m;
}
matrix create_aging_matrix(
int n_demo,
int n_age
) {
matrix[n_demo, n_demo] A = rep_matrix(0.0, n_demo, n_demo);
for(i in 2:n_age) {
A[i, (i - 1)] = 1.0;
}
A[n_age, n_age] = 1.0;
for(i in (n_age + 2):n_demo) {
A[i, (i - 1)] = 1.0;
}
A[n_demo, n_demo] = 1.0;
return A;
}
vector compute_baseline_birth_rate(
real min_bbr,
real max_bbr,
real h_int,
real h_slope,
real h_weight,
vector h1,
vector h2
) {
vector[dims(h1)[1]] weighted_h = ((h_weight * h1) + ((1.0 - h_weight) * h2));
return (min_bbr + ((max_bbr - min_bbr) * inv_logit((h_slope * (h_int + weighted_h)))));
}
real birth_rate_at_carrying_capacity(
real phi_a,
vector mu_m,
int n_age
) {
return ((2.0 * (1.0 - phi_a)) / exp(sum((-mu_m[1:(n_age - 1)]))));
}
real compute_density_dependence_intercept(
real max_bbr,
real dd_scaled
) {
return (-log((max_bbr + ((1.0 - max_bbr) * dd_scaled))));
}
vector initialize_population(
vector pop_init,
real pop_init_size,
int burn_in,
real birth_rate_year,
matrix aging,
vector S_diag,
int n_age
) {
int n_demo = dims(pop_init)[1];
if (dims(aging)[1] != n_demo) reject("initialize_population: dim mismatch — `aging` dim 1 (= ", dims(aging)[1], ") does not match `n_demo` (= ", n_demo, "), inferred from `pop_init` dim 1. `n_demo` sizes: `pop_init` dim 1 (= ", dims(pop_init)[1], "), `aging` dim 1 (= ", dims(aging)[1], "), `aging` dim 2 (= ", dims(aging)[2], "), `S_diag` dim 1 (= ", dims(S_diag)[1], ").");
if (dims(aging)[2] != n_demo) reject("initialize_population: dim mismatch — `aging` dim 2 (= ", dims(aging)[2], ") does not match `n_demo` (= ", n_demo, "), inferred from `pop_init` dim 1. `n_demo` sizes: `pop_init` dim 1 (= ", dims(pop_init)[1], "), `aging` dim 1 (= ", dims(aging)[1], "), `aging` dim 2 (= ", dims(aging)[2], "), `S_diag` dim 1 (= ", dims(S_diag)[1], ").");
if (dims(S_diag)[1] != n_demo) reject("initialize_population: dim mismatch — `S_diag` dim 1 (= ", dims(S_diag)[1], ") does not match `n_demo` (= ", n_demo, "), inferred from `pop_init` dim 1. `n_demo` sizes: `pop_init` dim 1 (= ", dims(pop_init)[1], "), `aging` dim 1 (= ", dims(aging)[1], "), `aging` dim 2 (= ", dims(aging)[2], "), `S_diag` dim 1 (= ", dims(S_diag)[1], ").");
vector[dims(pop_init)[1]] pop = pop_init;
for(k in 1:burn_in) {
pop = (aging * diag_matrix(S_diag) * pop);
pop[1] = ((birth_rate_year / 2.0) * pop[n_age]);
pop[(n_age + 1)] = ((birth_rate_year / 2.0) * pop[n_age]);
}
return ((pop * pop_init_size) / sum(pop));
}
tuple(vector, vector, vector, matrix, matrix, matrix, vector, vector, matrix) run_state_process(
int n_state_years,
int n_age,
vector pop_first,
real birth_rate_first,
real pop_total_first,
vector baseline_bbr,
real dd_intercept,
real dd_slope,
matrix aging,
vector S_diag,
vector mu_m,
vector hs_sw,
vector hs_fi,
array[] int hq_sw,
array[] int hq_fi,
real he_sd_sw,
real he_sd_fi,
vector eps_h_sw,
vector eps_h_fi,
real t_mate_to_preg,
real t_birth_to_end_hunt,
vector eps_birth,
vector eps_sex,
matrix transition_noise_raw,
vector pi_s,
vector pi_c,
real prob_of_ca,
vector ode_init_state,
vector ode_times
) {
int n_demo = dims(pop_first)[1];
if (dims(aging)[1] != n_demo) reject("run_state_process: dim mismatch — `aging` dim 1 (= ", dims(aging)[1], ") does not match `n_demo` (= ", n_demo, "), inferred from `pop_first` dim 1. `n_demo` sizes: `pop_first` dim 1 (= ", dims(pop_first)[1], "), `aging` dim 1 (= ", dims(aging)[1], "), `aging` dim 2 (= ", dims(aging)[2], "), `S_diag` dim 1 (= ", dims(S_diag)[1], "), `mu_m` dim 1 (= ", dims(mu_m)[1], "), `hs_sw` dim 1 (= ", dims(hs_sw)[1], "), `hs_fi` dim 1 (= ", dims(hs_fi)[1], ").");
if (dims(aging)[2] != n_demo) reject("run_state_process: dim mismatch — `aging` dim 2 (= ", dims(aging)[2], ") does not match `n_demo` (= ", n_demo, "), inferred from `pop_first` dim 1. `n_demo` sizes: `pop_first` dim 1 (= ", dims(pop_first)[1], "), `aging` dim 1 (= ", dims(aging)[1], "), `aging` dim 2 (= ", dims(aging)[2], "), `S_diag` dim 1 (= ", dims(S_diag)[1], "), `mu_m` dim 1 (= ", dims(mu_m)[1], "), `hs_sw` dim 1 (= ", dims(hs_sw)[1], "), `hs_fi` dim 1 (= ", dims(hs_fi)[1], ").");
if (dims(S_diag)[1] != n_demo) reject("run_state_process: dim mismatch — `S_diag` dim 1 (= ", dims(S_diag)[1], ") does not match `n_demo` (= ", n_demo, "), inferred from `pop_first` dim 1. `n_demo` sizes: `pop_first` dim 1 (= ", dims(pop_first)[1], "), `aging` dim 1 (= ", dims(aging)[1], "), `aging` dim 2 (= ", dims(aging)[2], "), `S_diag` dim 1 (= ", dims(S_diag)[1], "), `mu_m` dim 1 (= ", dims(mu_m)[1], "), `hs_sw` dim 1 (= ", dims(hs_sw)[1], "), `hs_fi` dim 1 (= ", dims(hs_fi)[1], ").");
if (dims(mu_m)[1] != n_demo) reject("run_state_process: dim mismatch — `mu_m` dim 1 (= ", dims(mu_m)[1], ") does not match `n_demo` (= ", n_demo, "), inferred from `pop_first` dim 1. `n_demo` sizes: `pop_first` dim 1 (= ", dims(pop_first)[1], "), `aging` dim 1 (= ", dims(aging)[1], "), `aging` dim 2 (= ", dims(aging)[2], "), `S_diag` dim 1 (= ", dims(S_diag)[1], "), `mu_m` dim 1 (= ", dims(mu_m)[1], "), `hs_sw` dim 1 (= ", dims(hs_sw)[1], "), `hs_fi` dim 1 (= ", dims(hs_fi)[1], ").");
if (dims(hs_sw)[1] != n_demo) reject("run_state_process: dim mismatch — `hs_sw` dim 1 (= ", dims(hs_sw)[1], ") does not match `n_demo` (= ", n_demo, "), inferred from `pop_first` dim 1. `n_demo` sizes: `pop_first` dim 1 (= ", dims(pop_first)[1], "), `aging` dim 1 (= ", dims(aging)[1], "), `aging` dim 2 (= ", dims(aging)[2], "), `S_diag` dim 1 (= ", dims(S_diag)[1], "), `mu_m` dim 1 (= ", dims(mu_m)[1], "), `hs_sw` dim 1 (= ", dims(hs_sw)[1], "), `hs_fi` dim 1 (= ", dims(hs_fi)[1], ").");
if (dims(hs_fi)[1] != n_demo) reject("run_state_process: dim mismatch — `hs_fi` dim 1 (= ", dims(hs_fi)[1], ") does not match `n_demo` (= ", n_demo, "), inferred from `pop_first` dim 1. `n_demo` sizes: `pop_first` dim 1 (= ", dims(pop_first)[1], "), `aging` dim 1 (= ", dims(aging)[1], "), `aging` dim 2 (= ", dims(aging)[2], "), `S_diag` dim 1 (= ", dims(S_diag)[1], "), `mu_m` dim 1 (= ", dims(mu_m)[1], "), `hs_sw` dim 1 (= ", dims(hs_sw)[1], "), `hs_fi` dim 1 (= ", dims(hs_fi)[1], ").");
n_demo = (2 * n_age);
array[dims(ode_times)[1]] real ode_ts = to_array_1d(ode_times);
vector[n_state_years] birth_rate;
vector[n_state_years] pregnancy_rate;
vector[n_state_years] population_total;
matrix[n_demo, n_state_years] hunted_sweden;
matrix[n_demo, n_state_years] hunted_finland;
matrix[n_demo, n_state_years] bycatch_expected;
vector[n_state_years] hunting_bag_total_sweden;
vector[n_state_years] hunting_bag_total_finland;
matrix[4, n_state_years] reproductive_probs;
matrix[n_demo, n_state_years] population_comp;
matrix[n_demo, n_state_years] survivors;
for(year in 1:n_state_years) {
if((year == 1)) {
birth_rate[year] = birth_rate_first;
population_comp[:, year] = pop_first;
population_total[year] = pop_total_first;
} else {
birth_rate[year] = update_birth_rate(baseline_bbr[year], dd_intercept, dd_slope, population_total[(year - 1)]);
population_comp[:, year] = update_population_from_survivors(
survivors[:, (year - 1)],
aging,
birth_rate[year],
eps_birth[year],
eps_sex[year],
n_age
);
population_total[year] = sum(population_comp[:, year]);
}
pregnancy_rate[year] = update_pregnancy_rate(
baseline_bbr[(year + 1)],
dd_intercept,
dd_slope,
population_total[year],
t_mate_to_preg
);
vector[n_demo] hp_sw;
vector[n_demo] hp_fi;
vector[n_demo] log_N = log(population_comp[:, year]);
real log_denom_sw = log_sum_exp((hs_sw + log_N));
real log_denom_fi = log_sum_exp((hs_fi + log_N));
if((hq_sw[year] == 0)) {
hp_sw = rep_vector(0.0, n_demo);
} else {
hp_sw = exp(
(
(
((hs_sw + log(hq_sw[year]) + log(2.0)) - (2.0 * log(t_birth_to_end_hunt))) -
(eps_h_sw[year] * he_sd_sw)
) -
log_denom_sw
)
);
}
if((hq_fi[year] == 0)) {
hp_fi = rep_vector(0.0, n_demo);
} else {
hp_fi = exp(
(
(
((hs_fi + log(hq_fi[year]) + log(2.0)) - (2.0 * log(t_birth_to_end_hunt))) -
(eps_h_fi[year] * he_sd_fi)
) -
log_denom_fi
)
);
}
vector[n_demo] exp_hunted_sw;
vector[n_demo] exp_hunted_fi;
for(demo in 1:n_demo) {
array[dims(ode_times)[1]] vector[dims(ode_init_state)[1]] sol_sw = ode_rk45_tol(
dH_dt,
ode_init_state,
0.0,
ode_ts,
1.0e-6,
1.0e-6,
1000,
population_comp[demo, year],
t_birth_to_end_hunt,
hp_sw[demo],
hp_fi[demo],
mu_m[demo]
);
exp_hunted_sw[demo] = sol_sw[1][1];
array[dims(ode_times)[1]] vector[dims(ode_init_state)[1]] sol_fi = ode_rk45_tol(
dH_dt,
ode_init_state,
0.0,
ode_ts,
1.0e-6,
1.0e-6,
1000,
population_comp[demo, year],
t_birth_to_end_hunt,
hp_fi[demo],
hp_sw[demo],
mu_m[demo]
);
exp_hunted_fi[demo] = sol_fi[1][1];
}
matrix[(4 * dims(S_diag)[1]), dims(S_diag)[1]] transition_matrix = create_transition_matrix(
exp_hunted_sw,
exp_hunted_fi,
population_comp[:, year],
hp_sw,
hp_fi,
t_birth_to_end_hunt,
S_diag
);
matrix[n_demo, 4] expected_fate = to_matrix((transition_matrix * population_comp[:, year]), n_demo, 4);
matrix[n_demo, 3] noise_year = to_matrix(transition_noise_raw[:, year], n_demo, 3);
matrix[n_demo, 4] realized_fate;
for(demo in 1:n_demo) {
realized_fate[demo, :] = multinomial_allocation(expected_fate[demo, :], noise_year[demo, :], population_comp[demo, year]);
}
survivors[:, year] = realized_fate[:, 1];
bycatch_expected[:, year] = realized_fate[:, 2];
hunted_sweden[:, year] = realized_fate[:, 3];
hunted_finland[:, year] = realized_fate[:, 4];
hunting_bag_total_sweden[year] = sum(hunted_sweden[:, year]);
hunting_bag_total_finland[year] = sum(hunted_finland[:, year]);
reproductive_probs[2, year] = (birth_rate[year] * pi_s[year] * (1.0 - pi_c[year]));
reproductive_probs[3, year] = (
(birth_rate[year] * (1.0 - pi_s[year]) * pi_c[year]) +
((1.0 - birth_rate[year]) * prob_of_ca * pi_c[year])
);
reproductive_probs[4, year] = (birth_rate[year] * pi_s[year] * pi_c[year]);
reproductive_probs[1, year] = (1.0 - sum(reproductive_probs[2:4, year]));
}
return (
birth_rate,
pregnancy_rate,
population_total,
hunted_sweden,
hunted_finland,
bycatch_expected,
hunting_bag_total_sweden,
hunting_bag_total_finland,
reproductive_probs
);
}
real update_birth_rate(
real b0,
real theta0,
real theta1,
real N_prev
) {
return (b0 * exp(((-theta0) * (exp((theta1 * N_prev)) - 1.0))));
}
vector update_population_from_survivors(
vector prev_surv,
matrix aging,
real birth_rate,
real eps_birth,
real eps_sex,
int n_age
) {
int n_demo = dims(prev_surv)[1];
if (dims(aging)[1] != n_demo) reject("update_population_from_survivors: dim mismatch — `aging` dim 1 (= ", dims(aging)[1], ") does not match `n_demo` (= ", n_demo, "), inferred from `prev_surv` dim 1. `n_demo` sizes: `prev_surv` dim 1 (= ", dims(prev_surv)[1], "), `aging` dim 1 (= ", dims(aging)[1], "), `aging` dim 2 (= ", dims(aging)[2], ").");
if (dims(aging)[2] != n_demo) reject("update_population_from_survivors: dim mismatch — `aging` dim 2 (= ", dims(aging)[2], ") does not match `n_demo` (= ", n_demo, "), inferred from `prev_surv` dim 1. `n_demo` sizes: `prev_surv` dim 1 (= ", dims(prev_surv)[1], "), `aging` dim 1 (= ", dims(aging)[1], "), `aging` dim 2 (= ", dims(aging)[2], ").");
vector[dims(aging)[1]] N = (aging * prev_surv);
real N_temp = ((N[n_age] * birth_rate) + (sqrt((N[n_age] * birth_rate * (1.0 - birth_rate))) * eps_birth));
N[1] = ((N_temp / 2.0) + (sqrt((N_temp / 4.0)) * eps_sex));
N[(n_age + 1)] = (N_temp - N[1]);
return N;
}
real update_pregnancy_rate(
real baseline_birth_rate,
real theta0,
real theta1,
real N_tot,
real tau_s
) {
return (baseline_birth_rate * exp((theta0 * (1.0 - (tau_s * exp((theta1 * N_tot)))))));
}
vector dH_dt(
real tau,
vector H,
real n0,
real k,
real E_1,
real E_2,
real mu
) {
int ny = dims(H)[1];
real surv = exp((((-(E_1 + E_2)) * ((k * tau) - ((tau * tau) / 2))) - (mu * tau)));
return rep_vector((n0 * E_1 * surv * (k - tau)), 1);
}
matrix create_transition_matrix(
vector hc_sw,
vector hc_fi,
vector pop,
vector hp_sw,
vector hp_fi,
real tau_h,
vector S_diag
) {
int n = dims(hc_sw)[1];
if (dims(hc_fi)[1] != n) reject("create_transition_matrix: dim mismatch — `hc_fi` dim 1 (= ", dims(hc_fi)[1], ") does not match `n` (= ", n, "), inferred from `hc_sw` dim 1. `n` sizes: `hc_sw` dim 1 (= ", dims(hc_sw)[1], "), `hc_fi` dim 1 (= ", dims(hc_fi)[1], "), `pop` dim 1 (= ", dims(pop)[1], "), `hp_sw` dim 1 (= ", dims(hp_sw)[1], "), `hp_fi` dim 1 (= ", dims(hp_fi)[1], "), `S_diag` dim 1 (= ", dims(S_diag)[1], ").");
if (dims(pop)[1] != n) reject("create_transition_matrix: dim mismatch — `pop` dim 1 (= ", dims(pop)[1], ") does not match `n` (= ", n, "), inferred from `hc_sw` dim 1. `n` sizes: `hc_sw` dim 1 (= ", dims(hc_sw)[1], "), `hc_fi` dim 1 (= ", dims(hc_fi)[1], "), `pop` dim 1 (= ", dims(pop)[1], "), `hp_sw` dim 1 (= ", dims(hp_sw)[1], "), `hp_fi` dim 1 (= ", dims(hp_fi)[1], "), `S_diag` dim 1 (= ", dims(S_diag)[1], ").");
if (dims(hp_sw)[1] != n) reject("create_transition_matrix: dim mismatch — `hp_sw` dim 1 (= ", dims(hp_sw)[1], ") does not match `n` (= ", n, "), inferred from `hc_sw` dim 1. `n` sizes: `hc_sw` dim 1 (= ", dims(hc_sw)[1], "), `hc_fi` dim 1 (= ", dims(hc_fi)[1], "), `pop` dim 1 (= ", dims(pop)[1], "), `hp_sw` dim 1 (= ", dims(hp_sw)[1], "), `hp_fi` dim 1 (= ", dims(hp_fi)[1], "), `S_diag` dim 1 (= ", dims(S_diag)[1], ").");
if (dims(hp_fi)[1] != n) reject("create_transition_matrix: dim mismatch — `hp_fi` dim 1 (= ", dims(hp_fi)[1], ") does not match `n` (= ", n, "), inferred from `hc_sw` dim 1. `n` sizes: `hc_sw` dim 1 (= ", dims(hc_sw)[1], "), `hc_fi` dim 1 (= ", dims(hc_fi)[1], "), `pop` dim 1 (= ", dims(pop)[1], "), `hp_sw` dim 1 (= ", dims(hp_sw)[1], "), `hp_fi` dim 1 (= ", dims(hp_fi)[1], "), `S_diag` dim 1 (= ", dims(S_diag)[1], ").");
if (dims(S_diag)[1] != n) reject("create_transition_matrix: dim mismatch — `S_diag` dim 1 (= ", dims(S_diag)[1], ") does not match `n` (= ", n, "), inferred from `hc_sw` dim 1. `n` sizes: `hc_sw` dim 1 (= ", dims(hc_sw)[1], "), `hc_fi` dim 1 (= ", dims(hc_fi)[1], "), `pop` dim 1 (= ", dims(pop)[1], "), `hp_sw` dim 1 (= ", dims(hp_sw)[1], "), `hp_fi` dim 1 (= ", dims(hp_fi)[1], "), `S_diag` dim 1 (= ", dims(S_diag)[1], ").");
vector[dims(pop)[1]] M_hunted_sw = (hc_sw ./ pop);
vector[dims(pop)[1]] M_hunted_fi = (hc_fi ./ pop);
vector[dims(S_diag)[1]] M_survived = (exp((((-(hp_sw + hp_fi)) * (tau_h * tau_h)) / 2.0)) .* S_diag);
vector[dims(S_diag)[1]] M_died = (((1.0 - M_survived) - M_hunted_sw) - M_hunted_fi);
return append_row(
diag_matrix(M_survived),
append_row(diag_matrix(M_died), append_row(diag_matrix(M_hunted_sw), diag_matrix(M_hunted_fi)))
);
}
row_vector multinomial_allocation(
row_vector eta_row,
row_vector u_row,
real N
) {
vector[dims(eta_row)[1]] eta = (eta_row');
vector[dims(eta_row)[1]] eta_adj = (eta * (1.0 + (1.0 / min(eta))));
row_vector[(1 + (4 - 2))] mean_logratio = ((digamma(eta_adj[2:4]) - digamma(eta_adj[1]))');
matrix[(1 + (4 - 2)), (1 + (4 - 2))] Sigma = (rep_matrix(trigamma(eta_adj[1]), 3, 3) + diag_matrix(trigamma(eta_adj[2:4])));
matrix[(1 + (4 - 2)), (1 + (4 - 2))] L = cholesky_decompose(Sigma);
row_vector[(1 + (1 + (4 - 2)))] logits = append_col(rep_row_vector(0.0, 1), (mean_logratio + (u_row * (L'))));
row_vector[(1 + (1 + (4 - 2)))] allocation = (softmax((logits'))');
return (allocation * N);
}
real aerial_count_lpmf(
array[] int obs,
array[] int year,
vector population_total,
real mu,
real phi
) {
int n_obs = dims(obs)[1];
if (dims(year)[1] != n_obs) reject("aerial_count_lpmf: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], ").");
return neg_binomial_2_lpmf(obs | (mu * population_total[year]), phi);
}
vector aerial_count_lpmfs(
array[] int obs,
array[] int year,
vector population_total,
real mu,
real phi
) {
int n_obs = dims(obs)[1];
if (dims(year)[1] != n_obs) reject("aerial_count_lpmfs: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], ").");
return neg_binomial_2_lpmfs(obs, (mu * population_total[year]), phi);
}
vector neg_binomial_2_lpmfs(
array[] int obs,
vector mu,
real phi
) {
return jbroadcasted_neg_binomial_2_lpmfs(obs, mu, phi);
}
vector jbroadcasted_neg_binomial_2_lpmfs(
array[] int x1,
vector x2,
real x3
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = neg_binomial_2_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i), x3);
}
return rv;
}
real neg_binomial_2_lpmfs(
int args1,
real args2,
real args3
) {
return neg_binomial_2_lpmf(args1 | args2, args3);
}
int broadcasted_getindex(array[] int x, int i) {
return x[i];
}
real broadcasted_getindex(vector x, int i) {
return x[i];
}
array[] int aerial_count_int_rng(
int anontok__1,
array[] int year,
vector population_total,
real mu,
real phi
) {
int n_obs = anontok__1;
if (dims(year)[1] != n_obs) reject("aerial_count_rng: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `anontok__1` dim 1. `n_obs` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `year` dim 1 (= ", dims(year)[1], ").");
return neg_binomial_2_rng((mu * population_total[year]), phi);
}
real bycatch_comp_lpmf(
array[, ] int obs,
array[] int year,
matrix bycatch_expected,
vector bias,
array[] int row_N
) {
int n_obs = dims(obs)[1];
int K = dims(obs)[2];
if (dims(year)[1] != n_obs) reject("bycatch_comp_lpmf: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(row_N)[1] != n_obs) reject("bycatch_comp_lpmf: dim mismatch — `row_N` dim 1 (= ", dims(row_N)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(bycatch_expected)[1] != K) reject("bycatch_comp_lpmf: dim mismatch — `bycatch_expected` dim 1 (= ", dims(bycatch_expected)[1], ") does not match `K` (= ", K, "), inferred from `obs` dim 2. `K` sizes: `obs` dim 2 (= ", dims(obs)[2], "), `bycatch_expected` dim 1 (= ", dims(bycatch_expected)[1], "), `bias` dim 1 (= ", dims(bias)[1], ").");
if (dims(bias)[1] != K) reject("bycatch_comp_lpmf: dim mismatch — `bias` dim 1 (= ", dims(bias)[1], ") does not match `K` (= ", K, "), inferred from `obs` dim 2. `K` sizes: `obs` dim 2 (= ", dims(obs)[2], "), `bycatch_expected` dim 1 (= ", dims(bycatch_expected)[1], "), `bias` dim 1 (= ", dims(bias)[1], ").");
vector[dims(bias)[1]] w = exp(bias);
real lp = 0.0;
for(i in 1:n_obs) {
int t = year[i];
lp += multinomial_lpmf(obs[i, :] |
((w .* bycatch_expected[:, t]) ./ dot_product(w, bycatch_expected[:, t]))
);
}
return lp;
}
vector bycatch_comp_lpmfs(
array[, ] int obs,
array[] int year,
matrix bycatch_expected,
vector bias,
array[] int row_N
) {
int n_obs = dims(obs)[1];
int K = dims(obs)[2];
if (dims(year)[1] != n_obs) reject("bycatch_comp_lpmfs: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(row_N)[1] != n_obs) reject("bycatch_comp_lpmfs: dim mismatch — `row_N` dim 1 (= ", dims(row_N)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(bycatch_expected)[1] != K) reject("bycatch_comp_lpmfs: dim mismatch — `bycatch_expected` dim 1 (= ", dims(bycatch_expected)[1], ") does not match `K` (= ", K, "), inferred from `obs` dim 2. `K` sizes: `obs` dim 2 (= ", dims(obs)[2], "), `bycatch_expected` dim 1 (= ", dims(bycatch_expected)[1], "), `bias` dim 1 (= ", dims(bias)[1], ").");
if (dims(bias)[1] != K) reject("bycatch_comp_lpmfs: dim mismatch — `bias` dim 1 (= ", dims(bias)[1], ") does not match `K` (= ", K, "), inferred from `obs` dim 2. `K` sizes: `obs` dim 2 (= ", dims(obs)[2], "), `bycatch_expected` dim 1 (= ", dims(bycatch_expected)[1], "), `bias` dim 1 (= ", dims(bias)[1], ").");
vector[dims(bias)[1]] w = exp(bias);
vector[n_obs] ll;
for(i in 1:n_obs) {
int t = year[i];
ll[i] = multinomial_lpmf(obs[i, :] |
((w .* bycatch_expected[:, t]) ./ dot_product(w, bycatch_expected[:, t]))
);
}
return ll;
}
array[, ] int bycatch_comp_int_rng(
tuple(int, int) anontok__1,
array[] int year,
matrix bycatch_expected,
vector bias,
array[] int row_N
) {
int n_obs = anontok__1.1;
int K = anontok__1.2;
if (dims(year)[1] != n_obs) reject("bycatch_comp_rng: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `anontok__1` dim 1. `n_obs` sizes: `anontok__1` dim 1 (= ", anontok__1.1, "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(row_N)[1] != n_obs) reject("bycatch_comp_rng: dim mismatch — `row_N` dim 1 (= ", dims(row_N)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `anontok__1` dim 1. `n_obs` sizes: `anontok__1` dim 1 (= ", anontok__1.1, "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(bycatch_expected)[1] != K) reject("bycatch_comp_rng: dim mismatch — `bycatch_expected` dim 1 (= ", dims(bycatch_expected)[1], ") does not match `K` (= ", K, "), inferred from `anontok__1` dim 2. `K` sizes: `anontok__1` dim 2 (= ", anontok__1.2, "), `bycatch_expected` dim 1 (= ", dims(bycatch_expected)[1], "), `bias` dim 1 (= ", dims(bias)[1], ").");
if (dims(bias)[1] != K) reject("bycatch_comp_rng: dim mismatch — `bias` dim 1 (= ", dims(bias)[1], ") does not match `K` (= ", K, "), inferred from `anontok__1` dim 2. `K` sizes: `anontok__1` dim 2 (= ", anontok__1.2, "), `bycatch_expected` dim 1 (= ", dims(bycatch_expected)[1], "), `bias` dim 1 (= ", dims(bias)[1], ").");
vector[dims(bias)[1]] w = exp(bias);
array[n_obs, K] int y;
for(i in 1:n_obs) {
int t = year[i];
y[i, :] = multinomial_rng(((w .* bycatch_expected[:, t]) ./ dot_product(w, bycatch_expected[:, t])), row_N[i]);
}
return y;
}
real harvest_bags_lpdf(
vector obs,
array[] int year,
vector hunted_total,
real cv
) {
int n_obs = dims(obs)[1];
if (dims(year)[1] != n_obs) reject("harvest_bags_lpdf: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], ").");
vector[dims(year)[1]] expected = hunted_total[year];
return normal_lpdf(obs | expected, (cv * expected));
}
vector harvest_bags_lpdfs(
vector obs,
array[] int year,
vector hunted_total,
real cv
) {
int n_obs = dims(obs)[1];
if (dims(year)[1] != n_obs) reject("harvest_bags_lpdfs: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], ").");
vector[dims(year)[1]] expected = hunted_total[year];
return normal_lpdfs(obs, expected, (cv * expected));
}
vector normal_lpdfs(
vector obs,
vector loc,
vector scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
vector x2,
vector x3
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = normal_lpdfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i), broadcasted_getindex(x3, i));
}
return rv;
}
real normal_lpdfs(
real args1,
real args2,
real args3
) {
return normal_lpdf(args1 | args2, args3);
}
vector harvest_bags_vector_rng(
int anontok__1,
array[] int year,
vector hunted_total,
real cv
) {
int n_obs = anontok__1;
if (dims(year)[1] != n_obs) reject("harvest_bags_rng: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `anontok__1` dim 1. `n_obs` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `year` dim 1 (= ", dims(year)[1], ").");
vector[dims(year)[1]] expected = hunted_total[year];
return to_vector(normal_rng(expected, (cv * expected)));
}
real hunting_comp_lpmf(
array[, ] int obs,
array[] int year,
matrix hunted,
vector hunted_total,
array[] int row_N
) {
int n_obs = dims(obs)[1];
int K = dims(obs)[2];
int T = dims(hunted)[2];
if (dims(year)[1] != n_obs) reject("hunting_comp_lpmf: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(row_N)[1] != n_obs) reject("hunting_comp_lpmf: dim mismatch — `row_N` dim 1 (= ", dims(row_N)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(hunted)[1] != K) reject("hunting_comp_lpmf: dim mismatch — `hunted` dim 1 (= ", dims(hunted)[1], ") does not match `K` (= ", K, "), inferred from `obs` dim 2. `K` sizes: `obs` dim 2 (= ", dims(obs)[2], "), `hunted` dim 1 (= ", dims(hunted)[1], ").");
if (dims(hunted_total)[1] != T) reject("hunting_comp_lpmf: dim mismatch — `hunted_total` dim 1 (= ", dims(hunted_total)[1], ") does not match `T` (= ", T, "), inferred from `hunted` dim 2. `T` sizes: `hunted` dim 2 (= ", dims(hunted)[2], "), `hunted_total` dim 1 (= ", dims(hunted_total)[1], ").");
real lp = 0.0;
for(i in 1:n_obs) {
int t = year[i];
lp += multinomial_lpmf(obs[i, :] | (hunted[:, t] ./ hunted_total[t]));
}
return lp;
}
vector hunting_comp_lpmfs(
array[, ] int obs,
array[] int year,
matrix hunted,
vector hunted_total,
array[] int row_N
) {
int n_obs = dims(obs)[1];
int K = dims(obs)[2];
int T = dims(hunted)[2];
if (dims(year)[1] != n_obs) reject("hunting_comp_lpmfs: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(row_N)[1] != n_obs) reject("hunting_comp_lpmfs: dim mismatch — `row_N` dim 1 (= ", dims(row_N)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(hunted)[1] != K) reject("hunting_comp_lpmfs: dim mismatch — `hunted` dim 1 (= ", dims(hunted)[1], ") does not match `K` (= ", K, "), inferred from `obs` dim 2. `K` sizes: `obs` dim 2 (= ", dims(obs)[2], "), `hunted` dim 1 (= ", dims(hunted)[1], ").");
if (dims(hunted_total)[1] != T) reject("hunting_comp_lpmfs: dim mismatch — `hunted_total` dim 1 (= ", dims(hunted_total)[1], ") does not match `T` (= ", T, "), inferred from `hunted` dim 2. `T` sizes: `hunted` dim 2 (= ", dims(hunted)[2], "), `hunted_total` dim 1 (= ", dims(hunted_total)[1], ").");
vector[n_obs] ll;
for(i in 1:n_obs) {
int t = year[i];
ll[i] = multinomial_lpmf(obs[i, :] | (hunted[:, t] ./ hunted_total[t]));
}
return ll;
}
array[, ] int hunting_comp_int_rng(
tuple(int, int) anontok__1,
array[] int year,
matrix hunted,
vector hunted_total,
array[] int row_N
) {
int n_obs = anontok__1.1;
int K = anontok__1.2;
int T = dims(hunted)[2];
if (dims(year)[1] != n_obs) reject("hunting_comp_rng: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `anontok__1` dim 1. `n_obs` sizes: `anontok__1` dim 1 (= ", anontok__1.1, "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(row_N)[1] != n_obs) reject("hunting_comp_rng: dim mismatch — `row_N` dim 1 (= ", dims(row_N)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `anontok__1` dim 1. `n_obs` sizes: `anontok__1` dim 1 (= ", anontok__1.1, "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(hunted)[1] != K) reject("hunting_comp_rng: dim mismatch — `hunted` dim 1 (= ", dims(hunted)[1], ") does not match `K` (= ", K, "), inferred from `anontok__1` dim 2. `K` sizes: `anontok__1` dim 2 (= ", anontok__1.2, "), `hunted` dim 1 (= ", dims(hunted)[1], ").");
if (dims(hunted_total)[1] != T) reject("hunting_comp_rng: dim mismatch — `hunted_total` dim 1 (= ", dims(hunted_total)[1], ") does not match `T` (= ", T, "), inferred from `hunted` dim 2. `T` sizes: `hunted` dim 2 (= ", dims(hunted)[2], "), `hunted_total` dim 1 (= ", dims(hunted_total)[1], ").");
array[n_obs, K] int y;
for(i in 1:n_obs) {
int t = year[i];
y[i, :] = multinomial_rng((hunted[:, t] ./ hunted_total[t]), row_N[i]);
}
return y;
}
real pregnancy_lpmf(
array[] int obs,
array[] int year,
array[] int sample_size,
vector pregnancy_rate
) {
int n_obs = dims(obs)[1];
if (dims(year)[1] != n_obs) reject("pregnancy_lpmf: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `sample_size` dim 1 (= ", dims(sample_size)[1], ").");
if (dims(sample_size)[1] != n_obs) reject("pregnancy_lpmf: dim mismatch — `sample_size` dim 1 (= ", dims(sample_size)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `sample_size` dim 1 (= ", dims(sample_size)[1], ").");
return binomial_lpmf(obs | sample_size, pregnancy_rate[year]);
}
vector pregnancy_lpmfs(
array[] int obs,
array[] int year,
array[] int sample_size,
vector pregnancy_rate
) {
int n_obs = dims(obs)[1];
if (dims(year)[1] != n_obs) reject("pregnancy_lpmfs: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `sample_size` dim 1 (= ", dims(sample_size)[1], ").");
if (dims(sample_size)[1] != n_obs) reject("pregnancy_lpmfs: dim mismatch — `sample_size` dim 1 (= ", dims(sample_size)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `sample_size` dim 1 (= ", dims(sample_size)[1], ").");
return binomial_lpmfs(obs, sample_size, pregnancy_rate[year]);
}
vector binomial_lpmfs(
array[] int y,
array[] int args1,
vector args2
) {
return jbroadcasted_binomial_lpmfs(y, args1, args2);
}
vector jbroadcasted_binomial_lpmfs(
array[] int x1,
array[] int x2,
vector x3
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = binomial_lpmfs(
broadcasted_getindex(x1, i),
broadcasted_getindex(x2, i),
broadcasted_getindex(x3, i)
);
}
return rv;
}
real binomial_lpmfs(
int args1,
int args2,
real args3
) {
return binomial_lpmf(args1 | args2, args3);
}
array[] int pregnancy_int_rng(
int anontok__1,
array[] int year,
array[] int sample_size,
vector pregnancy_rate
) {
int n_obs = anontok__1;
if (dims(year)[1] != n_obs) reject("pregnancy_rng: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `anontok__1` dim 1. `n_obs` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `year` dim 1 (= ", dims(year)[1], "), `sample_size` dim 1 (= ", dims(sample_size)[1], ").");
if (dims(sample_size)[1] != n_obs) reject("pregnancy_rng: dim mismatch — `sample_size` dim 1 (= ", dims(sample_size)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `anontok__1` dim 1. `n_obs` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `year` dim 1 (= ", dims(year)[1], "), `sample_size` dim 1 (= ", dims(sample_size)[1], ").");
return binomial_rng(sample_size, pregnancy_rate[year]);
}
real reproductive_signs_lpmf(
array[, ] int obs,
array[] int year,
matrix reproductive_probs,
array[] int row_N
) {
int n_obs = dims(obs)[1];
if (dims(year)[1] != n_obs) reject("reproductive_signs_lpmf: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(row_N)[1] != n_obs) reject("reproductive_signs_lpmf: dim mismatch — `row_N` dim 1 (= ", dims(row_N)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
real lp = 0.0;
for(i in 1:n_obs) {
lp += multinomial_lpmf(obs[i, :] | reproductive_probs[:, year[i]]);
}
return lp;
}
vector reproductive_signs_lpmfs(
array[, ] int obs,
array[] int year,
matrix reproductive_probs,
array[] int row_N
) {
int n_obs = dims(obs)[1];
if (dims(year)[1] != n_obs) reject("reproductive_signs_lpmfs: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(row_N)[1] != n_obs) reject("reproductive_signs_lpmfs: dim mismatch — `row_N` dim 1 (= ", dims(row_N)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `obs` dim 1. `n_obs` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
vector[n_obs] ll;
for(i in 1:n_obs) {
ll[i] = multinomial_lpmf(obs[i, :] | reproductive_probs[:, year[i]]);
}
return ll;
}
array[, ] int reproductive_signs_int_rng(
tuple(int, int) anontok__1,
array[] int year,
matrix reproductive_probs,
array[] int row_N
) {
int n_obs = anontok__1.1;
if (dims(year)[1] != n_obs) reject("reproductive_signs_rng: dim mismatch — `year` dim 1 (= ", dims(year)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `anontok__1` dim 1. `n_obs` sizes: `anontok__1` dim 1 (= ", anontok__1.1, "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
if (dims(row_N)[1] != n_obs) reject("reproductive_signs_rng: dim mismatch — `row_N` dim 1 (= ", dims(row_N)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `anontok__1` dim 1. `n_obs` sizes: `anontok__1` dim 1 (= ", anontok__1.1, "), `year` dim 1 (= ", dims(year)[1], "), `row_N` dim 1 (= ", dims(row_N)[1], ").");
array[n_obs, 4] int y;
for(i in 1:n_obs) {
y[i, :] = multinomial_rng(reproductive_probs[:, year[i]], row_N[i]);
}
return y;
}
}
data {
int n_demo;
int n_state_years;
int n_age;
int herring_index_1_n;
vector[herring_index_1_n] herring_index_1;
int herring_index_2_n;
vector[herring_index_2_n] herring_index_2;
int population_init_n;
vector[population_init_n] population_init;
int population_burn_in;
int hunting_quota_sweden_n;
array[hunting_quota_sweden_n] int hunting_quota_sweden;
int hunting_quota_finland_n;
array[hunting_quota_finland_n] int hunting_quota_finland;
real t_mate_to_preg;
real t_birth_to_end_hunt;
int ode_init_state_n;
vector[ode_init_state_n] ode_init_state;
int ode_times_n;
vector[ode_times_n] ode_times;
int obs_aerial_count_n;
array[obs_aerial_count_n] int obs_aerial_count;
int aerial_year_n;
array[aerial_year_n] int aerial_year;
int obs_bycatch_comp_m;
int obs_bycatch_comp_n;
array[obs_bycatch_comp_m, obs_bycatch_comp_n] int obs_bycatch_comp;
int bycatch_comp_sample_size_n;
int bycatch_comp_year_n;
array[bycatch_comp_year_n] int bycatch_comp_year;
array[bycatch_comp_sample_size_n] int bycatch_comp_sample_size;
int obs_hunting_bag_sweden_n;
vector[obs_hunting_bag_sweden_n] obs_hunting_bag_sweden;
int hunting_bag_year_sweden_n;
array[hunting_bag_year_sweden_n] int hunting_bag_year_sweden;
int obs_hunting_bag_finland_n;
vector[obs_hunting_bag_finland_n] obs_hunting_bag_finland;
int hunting_bag_year_finland_n;
array[hunting_bag_year_finland_n] int hunting_bag_year_finland;
int obs_hunting_comp_sweden_m;
int obs_hunting_comp_sweden_n;
array[obs_hunting_comp_sweden_m, obs_hunting_comp_sweden_n] int obs_hunting_comp_sweden;
int hunting_comp_sample_size_sweden_n;
int hunting_comp_year_sweden_n;
array[hunting_comp_year_sweden_n] int hunting_comp_year_sweden;
array[hunting_comp_sample_size_sweden_n] int hunting_comp_sample_size_sweden;
int obs_hunting_comp_finland_m;
int obs_hunting_comp_finland_n;
array[obs_hunting_comp_finland_m, obs_hunting_comp_finland_n] int obs_hunting_comp_finland;
int hunting_comp_sample_size_finland_n;
int hunting_comp_year_finland_n;
array[hunting_comp_year_finland_n] int hunting_comp_year_finland;
array[hunting_comp_sample_size_finland_n] int hunting_comp_sample_size_finland;
int obs_pregnancy_count_n;
array[obs_pregnancy_count_n] int obs_pregnancy_count;
int pregnancy_sample_size_n;
int pregnancy_count_year_n;
array[pregnancy_count_year_n] int pregnancy_count_year;
array[pregnancy_sample_size_n] int pregnancy_sample_size;
int obs_reproductive_signs_finland_m;
int obs_reproductive_signs_finland_n;
array[obs_reproductive_signs_finland_m, obs_reproductive_signs_finland_n] int obs_reproductive_signs_finland;
int reproductive_signs_sample_size_n;
int reproductive_signs_year_n;
array[reproductive_signs_year_n] int reproductive_signs_year;
array[reproductive_signs_sample_size_n] int reproductive_signs_sample_size;
}
transformed data {
matrix[n_demo, n_demo] aging = create_aging_matrix(n_demo, n_age);
}
parameters {
real<lower=0.0, upper=1.0> phi_a_sc;
real<lower=0.0, upper=1.0> phi_sc;
real<lower=0.0, upper=1.0> survival_shape;
real male_pup_survival_offset;
real male_adult_survival_offset;
real<lower=0.0> carrying_capacity;
real<lower=0.0, upper=1.0> max_baseline_birth_rate;
real<lower=0.0, upper=1.0> min_baseline_birth_rate_sc;
real herring_intercept_scaled;
real herring_slope;
real<lower=0.0, upper=1.0> herring_weight;
vector[n_demo] hunting_selectivity_sweden;
vector[n_demo] hunting_selectivity_finland;
real<lower=0.0> hunting_effort_sd_sweden;
real<lower=0.0> hunting_effort_sd_finland;
real<lower=0.0> population_init_size;
real<lower=0.0> harvest_bag_cv;
vector[n_state_years] epsilon_birth;
vector[n_state_years] epsilon_sex;
vector[n_state_years] epsilon_h_sw;
vector[n_state_years] epsilon_h_fi;
vector[(3 * n_demo * n_state_years)] transition_noise_vec;
real<lower=0.0, upper=1.0> report_ca_mean;
real<lower=0.0, upper=1.0> report_placental_mean;
real<lower=0.0, upper=1.0> prob_of_ca;
real<lower=0.0> report_placental_sd;
real<lower=0.0> report_ca_sd;
vector<lower=0.0>[n_state_years] epsilon_ca;
vector<lower=0.0>[n_state_years] epsilon_placental;
real<lower=0, upper=1> obs_aerial_count_aerial_count_mu;
real<lower=0.0> obs_aerial_count_aerial_count_overdispersion;
vector[n_demo] obs_bycatch_comp_bycatch_bias;
}
transformed parameters {
real<lower=0.0, upper=1.0> phi_a = phi_a_sc;
real phi_pup = (phi_sc * phi_a);
vector[(2 * n_age)] mu_m = mortality_rates(
phi_pup,
phi_a,
survival_shape,
n_age,
male_pup_survival_offset,
male_adult_survival_offset
);
vector[(2 * n_age)] S_diag = exp((-mu_m));
vector[herring_index_1_n] baseline_bbr = compute_baseline_birth_rate(
(max_baseline_birth_rate * min_baseline_birth_rate_sc),
max_baseline_birth_rate,
herring_intercept_scaled,
herring_slope,
herring_weight,
herring_index_1,
herring_index_2
);
real dd_scaled = birth_rate_at_carrying_capacity(phi_a, mu_m, n_age);
real dd_intercept = compute_density_dependence_intercept(max_baseline_birth_rate, dd_scaled);
real dd_slope = (-log(carrying_capacity));
vector[population_init_n] pop_first = initialize_population(
population_init,
population_init_size,
population_burn_in,
baseline_bbr[1],
aging,
S_diag,
n_age
);
vector[n_state_years] pi_s = (report_placental_mean * exp(((-epsilon_placental) * report_placental_sd)));
vector[n_state_years] pi_c = (report_ca_mean * exp(((-epsilon_ca) * report_ca_sd)));
matrix[(3 * n_demo), n_state_years] transition_noise_raw = to_matrix(transition_noise_vec, (3 * n_demo), n_state_years);
tuple(
vector[n_state_years],
vector[n_state_years],
vector[n_state_years],
matrix[n_demo, n_state_years],
matrix[n_demo, n_state_years],
matrix[n_demo, n_state_years],
vector[n_state_years],
vector[n_state_years],
matrix[4, n_state_years]
) state = run_state_process(
n_state_years,
n_age,
pop_first,
baseline_bbr[1],
sum(pop_first),
baseline_bbr,
dd_intercept,
dd_slope,
aging,
S_diag,
mu_m,
hunting_selectivity_sweden,
hunting_selectivity_finland,
hunting_quota_sweden,
hunting_quota_finland,
hunting_effort_sd_sweden,
hunting_effort_sd_finland,
epsilon_h_sw,
epsilon_h_fi,
t_mate_to_preg,
t_birth_to_end_hunt,
epsilon_birth,
epsilon_sex,
transition_noise_raw,
pi_s,
pi_c,
prob_of_ca,
ode_init_state,
ode_times
);
}
model {
phi_a_sc ~ uniform(0.0, 1.0);
phi_sc ~ uniform(0.0, 1.0);
survival_shape ~ uniform(0.0, 1.0);
male_pup_survival_offset ~ cauchy(0.0, 1.0);
male_adult_survival_offset ~ cauchy(0.0, 1.0);
carrying_capacity ~ lognormal(11.0, 1.0);
max_baseline_birth_rate ~ uniform(0.0, 1.0);
min_baseline_birth_rate_sc ~ uniform(0.0, 1.0);
herring_intercept_scaled ~ normal(0.0, 1.0);
herring_slope ~ normal(0.0, 1.0);
herring_weight ~ uniform(0.0, 1.0);
hunting_selectivity_sweden ~ normal(0.0, 1.0);
hunting_selectivity_finland ~ normal(0.0, 1.0);
hunting_effort_sd_sweden ~ cauchy(0.0, 1.0);
hunting_effort_sd_finland ~ cauchy(0.0, 1.0);
population_init_size ~ lognormal(11.0, 1.0);
harvest_bag_cv ~ lognormal(0.0, 1.0);
epsilon_birth ~ std_normal();
epsilon_sex ~ std_normal();
epsilon_h_sw ~ std_normal();
epsilon_h_fi ~ std_normal();
transition_noise_vec ~ std_normal();
report_ca_mean ~ uniform(0.0, 1.0);
report_placental_mean ~ uniform(0.0, 1.0);
prob_of_ca ~ uniform(0.0, 1.0);
report_placental_sd ~ normal(0.0, 0.1);
report_ca_sd ~ normal(0.0, 0.1);
epsilon_ca ~ std_normal();
epsilon_placental ~ std_normal();
obs_aerial_count_aerial_count_mu ~ beta(2.0, 2.0);
obs_aerial_count_aerial_count_overdispersion ~ lognormal(0.0, 1.0);
obs_aerial_count ~ aerial_count(
aerial_year,
state.3,
obs_aerial_count_aerial_count_mu,
obs_aerial_count_aerial_count_overdispersion
);
obs_bycatch_comp_bycatch_bias ~ normal(0.0, 1.0);
obs_bycatch_comp ~ bycatch_comp(bycatch_comp_year, state.6, obs_bycatch_comp_bycatch_bias, bycatch_comp_sample_size);
obs_hunting_bag_sweden ~ harvest_bags(hunting_bag_year_sweden, state.7, harvest_bag_cv);
obs_hunting_bag_finland ~ harvest_bags(hunting_bag_year_finland, state.8, harvest_bag_cv);
obs_hunting_comp_sweden ~ hunting_comp(hunting_comp_year_sweden, state.4, state.7, hunting_comp_sample_size_sweden);
obs_hunting_comp_finland ~ hunting_comp(hunting_comp_year_finland, state.5, state.8, hunting_comp_sample_size_finland);
obs_pregnancy_count ~ pregnancy(pregnancy_count_year, pregnancy_sample_size, state.2);
obs_reproductive_signs_finland ~ reproductive_signs(reproductive_signs_year, state.9, reproductive_signs_sample_size);
}
generated quantities {
vector[aerial_year_n] obs_aerial_count_likelihood = aerial_count_lpmfs(
obs_aerial_count,
aerial_year,
state.3,
obs_aerial_count_aerial_count_mu,
obs_aerial_count_aerial_count_overdispersion
);
array[aerial_year_n] int obs_aerial_count_gen = aerial_count_int_rng(
obs_aerial_count_n,
aerial_year,
state.3,
obs_aerial_count_aerial_count_mu,
obs_aerial_count_aerial_count_overdispersion
);
vector[bycatch_comp_sample_size_n] obs_bycatch_comp_likelihood = bycatch_comp_lpmfs(
obs_bycatch_comp,
bycatch_comp_year,
state.6,
obs_bycatch_comp_bycatch_bias,
bycatch_comp_sample_size
);
array[bycatch_comp_sample_size_n, n_demo] int obs_bycatch_comp_gen = bycatch_comp_int_rng(
(obs_bycatch_comp_m, obs_bycatch_comp_n),
bycatch_comp_year,
state.6,
obs_bycatch_comp_bycatch_bias,
bycatch_comp_sample_size
);
vector[hunting_bag_year_sweden_n] obs_hunting_bag_sweden_likelihood = harvest_bags_lpdfs(obs_hunting_bag_sweden, hunting_bag_year_sweden, state.7, harvest_bag_cv);
vector[hunting_bag_year_sweden_n] obs_hunting_bag_sweden_gen = harvest_bags_vector_rng(obs_hunting_bag_sweden_n, hunting_bag_year_sweden, state.7, harvest_bag_cv);
vector[hunting_bag_year_finland_n] obs_hunting_bag_finland_likelihood = harvest_bags_lpdfs(obs_hunting_bag_finland, hunting_bag_year_finland, state.8, harvest_bag_cv);
vector[hunting_bag_year_finland_n] obs_hunting_bag_finland_gen = harvest_bags_vector_rng(
obs_hunting_bag_finland_n,
hunting_bag_year_finland,
state.8,
harvest_bag_cv
);
vector[hunting_comp_sample_size_sweden_n] obs_hunting_comp_sweden_likelihood = hunting_comp_lpmfs(
obs_hunting_comp_sweden,
hunting_comp_year_sweden,
state.4,
state.7,
hunting_comp_sample_size_sweden
);
array[hunting_comp_sample_size_sweden_n, n_demo] int obs_hunting_comp_sweden_gen = hunting_comp_int_rng(
(obs_hunting_comp_sweden_m, obs_hunting_comp_sweden_n),
hunting_comp_year_sweden,
state.4,
state.7,
hunting_comp_sample_size_sweden
);
vector[hunting_comp_sample_size_finland_n] obs_hunting_comp_finland_likelihood = hunting_comp_lpmfs(
obs_hunting_comp_finland,
hunting_comp_year_finland,
state.5,
state.8,
hunting_comp_sample_size_finland
);
array[hunting_comp_sample_size_finland_n, n_demo] int obs_hunting_comp_finland_gen = hunting_comp_int_rng(
(obs_hunting_comp_finland_m, obs_hunting_comp_finland_n),
hunting_comp_year_finland,
state.5,
state.8,
hunting_comp_sample_size_finland
);
vector[pregnancy_sample_size_n] obs_pregnancy_count_likelihood = pregnancy_lpmfs(obs_pregnancy_count, pregnancy_count_year, pregnancy_sample_size, state.2);
array[pregnancy_sample_size_n] int obs_pregnancy_count_gen = pregnancy_int_rng(obs_pregnancy_count_n, pregnancy_count_year, pregnancy_sample_size, state.2);
vector[reproductive_signs_sample_size_n] obs_reproductive_signs_finland_likelihood = reproductive_signs_lpmfs(
obs_reproductive_signs_finland,
reproductive_signs_year,
state.9,
reproductive_signs_sample_size
);
array[reproductive_signs_sample_size_n, 4] int obs_reproductive_signs_finland_gen = reproductive_signs_int_rng(
(obs_reproductive_signs_finland_m, obs_reproductive_signs_finland_n),
reproductive_signs_year,
state.9,
reproductive_signs_sample_size
);
}The Julia panes above are, in order: the parent @slic body, then four representative @deffun cards — the state-process scan, the ODE right-hand side, the fate-allocation helper, and one custom-family density head. The full library (all 31 UDF cards + the two observation-stream submodels) is vendored verbatim in docs/grey_seal_ipm.jl; they all appear in the generated Stan's functions {} block. The Stan pane is the complete emitted program.
What this exercises
compile_slic_bundle— a multi-source workspace (31 UDF cards + two anonymous submodels + a parent body) assembled and traced in one call, the natural shape for a model too large for a single inline block.A year-recursive scan (
run_state_process) in a@deffun, with named-tuple state carriers read by field (state.population_total).A numerical ODE inside the scan —
ode_rk45_tolover a@deffunhunting-hazard right-hand side, per demographic class per year.A custom Dirichlet-multinomial fate allocation (
multinomial_allocation) usingdigamma/trigamma/cholesky_decompose/softmax.Six user-defined distribution families (
_lpmf/_lpmfs/_rngtriads with@lhs @lpxf) driving negative-binomial, normal, binomial, and multinomial observation likelihoods — each with automatic posterior-predictive and pointwise-log-likelihood twins.Observation submodels (
data ~ submodel(...)) and direct-family observations (data ~ family(...)), eight streams in total, each added or dropped as one line.A full executable descriptor — the assembled model offers
fit/predict/pointwise_loglik.