Grey-seal IPM on the @brm formula surface
This page ports the full Baltic grey-seal integrated population model — a joint state-space model of an age/sex-structured population fit to eight heterogeneous observation streams — onto the @brm formula surface. It is a faithful port of the verbatim SlicTranspiler research model (the same one rendered as a StanBlocks-native case study): every regression / random-effect / covariate piece is hoisted onto a brms-style formula, while the mechanistic core stays a @deffun scan — exactly as the CDC ww-inference model keeps its renewal scan in a @deffun. It is rendered in the standard four-pane view and verified (transpile + stanc + compiles) by test/seal_brm.jl.
What is a brms formula here
Hoisted from the upstream @slic parent onto the formula surface:
Birth-rate covariate regression —
bbr_logit ~ 1 + herring_index_1 + herring_index_2(the upstreamcompute_baseline_birth_ratelogistic-on-herring), bounded inside the scan.Per-year process/effort random effects, with estimated group SDs —
eps_birth ~ 0 + (1|year),eps_sex,eps_h_sw,eps_h_fi,eps_ca,eps_placental. The upstreamepsilon_*fixed-unitstd_normalinnovations become proper(1|year)random effects; the estimated group SD is the brms upgrade (and the reason the@brmform has a few more parameters than the@slicone).Per-demographic-class hunting selectivity + bycatch bias —
hs_sw ~ 0 + (1|demo),hs_fi,bycatch_bias.Per-cell transition process noise —
tnoise ~ 0 + (1|noise_cell), a crossed cell random effect reshaped into the transition-noise matrix in the scan.
These live on four different frames (year, demo, noise_cell, plus the herring covariate frame) simultaneously — multi-frame is fine on the formula surface.
What stays a @deffun scan
The carried-state recurrence — the right home for it, as CDC's renewal is: the verbatim run_state_process (Leslie aging, size-structured mortality, density-dependent births, a within-year hunting ODE ode_rk45_tol on dH_dt, and a multinomial_allocation logistic-normal survived/bycatch/hunted split) and its twelve verbatim UDFs. A seal_state wrapper runs the parent's transformed quantities + the scan and returns a NamedTuple of carriers, each observation reading a named field (state.population_total, state.hunted_sweden, …) via @brm body field access — one scan, exactly as the @slic seal reads its run_state_process struct.
The eight observation streams
Each is a @brm family likelihood reading a state.<field> indexed to that stream's observation years, on its own frame:
aerial pup counts —
NegativeBinomial2onstate.population_total;Swedish / Finnish harvest bags —
Normal(cv-scaled sd) on the hunting-bag totals;Swedish / Finnish hunting age-composition, bycatch composition, reproductive-signs composition — per-row
Multinomialon the per-year simplices derived from the hunted / bycatch / reproductive carriers;pregnancy counts —
Binomialonstate.pregnancy_rate.
The hidden setup below evaluates the @deffun scan + helpers and the fixture from the checked-in source before the declaration is rendered.
The declaration below is extracted verbatim from grey_seal_brm_model in the reviewed source, followed by the StanBlocks model it emits, the generated Stan (whose functions {} block carries the verbatim scan + UDFs + the composition densities), and the Turing pane.
brm-comparison
Full grey-seal IPM (@brm formula surface)julia
function grey_seal_brm_model(df = grey_seal_brm_fixture())
@brm df begin
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)
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)
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)
aerial_mu ~ Uniform(0.0, 1.0); phi_aerial ~ LogNormal(0.0, 1.0)
harvest_bag_cv ~ LogNormal(0.0, 1.0; lower = 0.0)
# brms formulas
bbr_logit ~ 1 + herring_index_1 + herring_index_2
eps_birth ~ 0 + (1 | year); eps_sex ~ 0 + (1 | year)
eps_h_sw ~ 0 + (1 | year); eps_h_fi ~ 0 + (1 | year)
eps_ca ~ 0 + (1 | year); eps_placental ~ 0 + (1 | year)
hs_sw ~ 0 + (1 | demo); hs_fi ~ 0 + (1 | demo)
bycatch_bias ~ 0 + (1 | demo)
tnoise ~ 0 + (1 | noise_cell)
state = seal_state(bbr_logit, eps_birth, eps_sex, eps_h_sw, eps_h_fi, eps_ca, eps_placental,
hs_sw, hs_fi, tnoise, phi_a_sc, phi_sc, survival_shape,
male_pup_survival_offset, male_adult_survival_offset, carrying_capacity,
max_baseline_birth_rate, min_baseline_birth_rate_sc,
hunting_effort_sd_sweden, hunting_effort_sd_finland, population_init_size,
report_ca_mean, report_placental_mean, report_ca_sd, report_placental_sd,
prob_of_ca, population_init, hunting_quota_sweden, hunting_quota_finland,
t_mate_to_preg, t_birth_to_end_hunt, population_burn_in, n_age, ode_init_state, ode_times)
obs_aerial_count ~ NegativeBinomial2(aerial_mean(aerial_mu, state.population_total, aerial_year), phi_aerial)
obs_hunting_bag_sweden ~ Normal(harvest_mean(state.hunting_bag_total_sweden, hunting_bag_year_sweden), harvest_sd(state.hunting_bag_total_sweden, hunting_bag_year_sweden, harvest_bag_cv))
obs_hunting_bag_finland ~ Normal(harvest_mean(state.hunting_bag_total_finland, hunting_bag_year_finland), harvest_sd(state.hunting_bag_total_finland, hunting_bag_year_finland, harvest_bag_cv))
obs_hunting_comp_sweden ~ Multinomial(hunting_comp_sample_size_sweden, comp_hunted(state.hunted_sweden, state.hunting_bag_total_sweden, hunting_comp_year_sweden))
obs_hunting_comp_finland ~ Multinomial(hunting_comp_sample_size_finland, comp_hunted(state.hunted_finland, state.hunting_bag_total_finland, hunting_comp_year_finland))
obs_bycatch_comp ~ Multinomial(bycatch_comp_sample_size, comp_bycatch(state.bycatch_expected, bycatch_bias, bycatch_comp_year))
obs_pregnancy_count ~ Binomial(pregnancy_sample_size, state.pregnancy_rate[pregnancy_count_year])
obs_reproductive_signs_finland ~ Multinomial(reproductive_signs_sample_size, comp_repro(state.reproductive_probs, reproductive_signs_year))
end
endjulia
BRMI:
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)
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)
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)
aerial_mu ~ Uniform(0.0, 1.0)
phi_aerial ~ LogNormal(0.0, 1.0)
harvest_bag_cv ~ LogNormal(0.0, 1.0; lower=0.0)
herring_index_1: data (eltype=Float64, n=4)
herring_index_2: data (eltype=Float64, n=4)
bbr_logit ~ 1 + herring_index_1 + herring_index_2
year: data (eltype=Int64, n=3)
eps_birth ~ 0 + (1 | year)
eps_sex ~ 0 + (1 | year)
eps_h_sw ~ 0 + (1 | year)
eps_h_fi ~ 0 + (1 | year)
eps_ca ~ 0 + (1 | year)
eps_placental ~ 0 + (1 | year)
demo: data (eltype=Int64, n=6)
hs_sw ~ 0 + (1 | demo)
hs_fi ~ 0 + (1 | demo)
bycatch_bias ~ 0 + (1 | demo)
noise_cell: data (eltype=Int64, n=54)
tnoise ~ 0 + (1 | noise_cell)
population_init: data (eltype=Float64, n=6)
hunting_quota_sweden: data (eltype=Int64, n=3)
hunting_quota_finland: data (eltype=Int64, n=3)
t_mate_to_preg: data (eltype=Float64, n=1)
t_birth_to_end_hunt: data (eltype=Float64, n=1)
population_burn_in: data (eltype=Int64, n=1)
n_age: data (eltype=Int64, n=1)
ode_init_state: data (eltype=Float64, n=1)
ode_times: data (eltype=Float64, n=1)
:state = seal_state(bbr_logit, eps_birth, eps_sex, eps_h_sw, eps_h_fi, eps_ca, eps_placental, hs_sw, hs_fi, tnoise, phi_a_sc, phi_sc, survival_shape, male_pup_survival_offset, male_adult_survival_offset, carrying_capacity, max_baseline_birth_rate, min_baseline_birth_rate_sc, hunting_effort_sd_sweden, hunting_effort_sd_finland, population_init_size, report_ca_mean, report_placental_mean, report_ca_sd, report_placental_sd, prob_of_ca, population_init, hunting_quota_sweden, hunting_quota_finland, t_mate_to_preg, t_birth_to_end_hunt, population_burn_in, n_age, ode_init_state, ode_times)
aerial_year: data (eltype=Int64, n=3)
obs_aerial_count ~ NegativeBinomial2(aerial_mean(aerial_mu, getproperty(state, :population_total), aerial_year), phi_aerial)
hunting_bag_year_sweden: data (eltype=Int64, n=3)
obs_hunting_bag_sweden ~ Normal(harvest_mean(getproperty(state, :hunting_bag_total_sweden), hunting_bag_year_sweden), harvest_sd(getproperty(state, :hunting_bag_total_sweden), hunting_bag_year_sweden, harvest_bag_cv))
hunting_bag_year_finland: data (eltype=Int64, n=3)
obs_hunting_bag_finland ~ Normal(harvest_mean(getproperty(state, :hunting_bag_total_finland), hunting_bag_year_finland), harvest_sd(getproperty(state, :hunting_bag_total_finland), hunting_bag_year_finland, harvest_bag_cv))
hunting_comp_sample_size_sweden: data (eltype=Int64, n=3)
hunting_comp_year_sweden: data (eltype=Int64, n=3)
obs_hunting_comp_sweden ~ Multinomial(hunting_comp_sample_size_sweden, comp_hunted(getproperty(state, :hunted_sweden), getproperty(state, :hunting_bag_total_sweden), hunting_comp_year_sweden))
hunting_comp_sample_size_finland: data (eltype=Int64, n=3)
hunting_comp_year_finland: data (eltype=Int64, n=3)
obs_hunting_comp_finland ~ Multinomial(hunting_comp_sample_size_finland, comp_hunted(getproperty(state, :hunted_finland), getproperty(state, :hunting_bag_total_finland), hunting_comp_year_finland))
bycatch_comp_sample_size: data (eltype=Int64, n=3)
bycatch_comp_year: data (eltype=Int64, n=3)
obs_bycatch_comp ~ Multinomial(bycatch_comp_sample_size, comp_bycatch(getproperty(state, :bycatch_expected), bycatch_bias, bycatch_comp_year))
pregnancy_sample_size: data (eltype=Int64, n=3)
pregnancy_count_year: data (eltype=Int64, n=3)
obs_pregnancy_count ~ Binomial(pregnancy_sample_size, getindex(getproperty(state, :pregnancy_rate), pregnancy_count_year))
reproductive_signs_sample_size: data (eltype=Int64, n=3)
reproductive_signs_year: data (eltype=Int64, n=3)
obs_reproductive_signs_finland ~ Multinomial(reproductive_signs_sample_size, comp_repro(getproperty(state, :reproductive_probs), reproductive_signs_year))julia
SBBRMI with data keys = [:aerial_year, :bycatch_comp_sample_size, :bycatch_comp_year, :demo, :demo_idx, :herring_index_1, :herring_index_2, :hunting_bag_year_finland, :hunting_bag_year_sweden, :hunting_comp_sample_size_finland, :hunting_comp_sample_size_sweden, :hunting_comp_year_finland, :hunting_comp_year_sweden, :hunting_quota_finland, :hunting_quota_sweden, :n_age, :n_demo, :n_noise_cell, :n_year, :noise_cell, :noise_cell_idx, :obs_aerial_count, :obs_bycatch_comp, :obs_hunting_bag_finland, :obs_hunting_bag_sweden, :obs_hunting_comp_finland, :obs_hunting_comp_sweden, :obs_pregnancy_count, :obs_reproductive_signs_finland, :ode_init_state, :ode_times, :population_burn_in, :population_init, :pregnancy_count_year, :pregnancy_sample_size, :reproductive_signs_sample_size, :reproductive_signs_year, :t_birth_to_end_hunt, :t_mate_to_preg, :year, :year_idx]
emitted @slic body:
begin
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)
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)
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)
aerial_mu ~ uniform(0.0, 1.0)
phi_aerial ~ lognormal(0.0, 1.0)
harvest_bag_cv ~ lognormal(0.0, 1.0; lower = 0.0)
X_bbr_logit = hcat(rep_vector(1.0, num_elements(herring_index_1)), herring_index_1, herring_index_2)
pop_bbr_logit ~ popefs(; X = X_bbr_logit)
bbr_logit = pop_bbr_logit
r_eps_birth_year ~ ranef_intercept(; group_idx = year_idx, n_groups = n_year)
eps_birth = r_eps_birth_year
r_eps_sex_year ~ ranef_intercept(; group_idx = year_idx, n_groups = n_year)
eps_sex = r_eps_sex_year
r_eps_h_sw_year ~ ranef_intercept(; group_idx = year_idx, n_groups = n_year)
eps_h_sw = r_eps_h_sw_year
r_eps_h_fi_year ~ ranef_intercept(; group_idx = year_idx, n_groups = n_year)
eps_h_fi = r_eps_h_fi_year
r_eps_ca_year ~ ranef_intercept(; group_idx = year_idx, n_groups = n_year)
eps_ca = r_eps_ca_year
r_eps_placental_year ~ ranef_intercept(; group_idx = year_idx, n_groups = n_year)
eps_placental = r_eps_placental_year
r_hs_sw_demo ~ ranef_intercept(; group_idx = demo_idx, n_groups = n_demo)
hs_sw = r_hs_sw_demo
r_hs_fi_demo ~ ranef_intercept(; group_idx = demo_idx, n_groups = n_demo)
hs_fi = r_hs_fi_demo
r_bycatch_bias_demo ~ ranef_intercept(; group_idx = demo_idx, n_groups = n_demo)
bycatch_bias = r_bycatch_bias_demo
r_tnoise_noise_cell ~ ranef_intercept(; group_idx = noise_cell_idx, n_groups = n_noise_cell)
tnoise = r_tnoise_noise_cell
state = (Main.seal_brm.seal_state)(bbr_logit, eps_birth, eps_sex, eps_h_sw, eps_h_fi, eps_ca, eps_placental, hs_sw, hs_fi, tnoise, phi_a_sc, phi_sc, survival_shape, male_pup_survival_offset, male_adult_survival_offset, carrying_capacity, max_baseline_birth_rate, min_baseline_birth_rate_sc, hunting_effort_sd_sweden, hunting_effort_sd_finland, population_init_size, report_ca_mean, report_placental_mean, report_ca_sd, report_placental_sd, prob_of_ca, population_init, hunting_quota_sweden, hunting_quota_finland, t_mate_to_preg, t_birth_to_end_hunt, population_burn_in, n_age, ode_init_state, ode_times)
obs_aerial_count ~ neg_binomial_2((Main.seal_brm.aerial_mean)(aerial_mu, state.population_total, aerial_year), phi_aerial)
obs_hunting_bag_sweden ~ normal((Main.seal_brm.harvest_mean)(state.hunting_bag_total_sweden, hunting_bag_year_sweden), (Main.seal_brm.harvest_sd)(state.hunting_bag_total_sweden, hunting_bag_year_sweden, harvest_bag_cv))
obs_hunting_bag_finland ~ normal((Main.seal_brm.harvest_mean)(state.hunting_bag_total_finland, hunting_bag_year_finland), (Main.seal_brm.harvest_sd)(state.hunting_bag_total_finland, hunting_bag_year_finland, harvest_bag_cv))
obs_hunting_comp_sweden ~ brm_multinomial((Main.seal_brm.comp_hunted)(state.hunted_sweden, state.hunting_bag_total_sweden, hunting_comp_year_sweden), hunting_comp_sample_size_sweden)
obs_hunting_comp_finland ~ brm_multinomial((Main.seal_brm.comp_hunted)(state.hunted_finland, state.hunting_bag_total_finland, hunting_comp_year_finland), hunting_comp_sample_size_finland)
obs_bycatch_comp ~ brm_multinomial((Main.seal_brm.comp_bycatch)(state.bycatch_expected, bycatch_bias, bycatch_comp_year), bycatch_comp_sample_size)
obs_pregnancy_count ~ binomial(pregnancy_sample_size, state.pregnancy_rate[pregnancy_count_year])
obs_reproductive_signs_finland ~ brm_multinomial((Main.seal_brm.comp_repro)(state.reproductive_probs, reproductive_signs_year), reproductive_signs_sample_size)
endstan
functions {
matrix hcat(vector x, vector y, vector z) {
return hcat(hcat(x, y), z);
}
matrix hcat(
matrix x,
vector y
) {
int m = dims(x)[1];
int n = dims(x)[2];
if (dims(y)[1] != m) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `m` (= ", m, "), inferred from `x` dim 1. `m` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
matrix hcat(
vector x,
vector y
) {
int n = dims(x)[1];
if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
return append_col(x, y);
}
tuple(vector, vector, vector, matrix, matrix, matrix, vector, vector, matrix) seal_state(
vector bbr_logit,
vector eps_birth,
vector eps_sex,
vector eps_h_sw,
vector eps_h_fi,
vector eps_ca,
vector eps_placental,
vector hs_sw,
vector hs_fi,
vector tnoise,
real phi_a_sc,
real phi_sc,
real survival_shape,
real male_pup_offset,
real male_adult_offset,
real carrying_capacity,
real max_bbr,
real min_bbr_sc,
real he_sd_sw,
real he_sd_fi,
real pop_init_size,
real report_ca_mean,
real report_placental_mean,
real report_ca_sd,
real report_placental_sd,
real prob_of_ca,
vector population_init,
array[] int hq_sw,
array[] int hq_fi,
real t_mate,
real t_hunt,
int burn_in,
int n_age,
vector ode_init,
vector ode_times
) {
int T = dims(eps_birth)[1];
int nd = dims(hs_sw)[1];
if (dims(eps_sex)[1] != T) reject("seal_state: dim mismatch — `eps_sex` dim 1 (= ", dims(eps_sex)[1], ") does not match `T` (= ", T, "), inferred from `eps_birth` dim 1. `T` sizes: `eps_birth` dim 1 (= ", dims(eps_birth)[1], "), `eps_sex` dim 1 (= ", dims(eps_sex)[1], "), `eps_h_sw` dim 1 (= ", dims(eps_h_sw)[1], "), `eps_h_fi` dim 1 (= ", dims(eps_h_fi)[1], "), `eps_ca` dim 1 (= ", dims(eps_ca)[1], "), `eps_placental` dim 1 (= ", dims(eps_placental)[1], "), `hq_sw` dim 1 (= ", dims(hq_sw)[1], "), `hq_fi` dim 1 (= ", dims(hq_fi)[1], ").");
if (dims(eps_h_sw)[1] != T) reject("seal_state: dim mismatch — `eps_h_sw` dim 1 (= ", dims(eps_h_sw)[1], ") does not match `T` (= ", T, "), inferred from `eps_birth` dim 1. `T` sizes: `eps_birth` dim 1 (= ", dims(eps_birth)[1], "), `eps_sex` dim 1 (= ", dims(eps_sex)[1], "), `eps_h_sw` dim 1 (= ", dims(eps_h_sw)[1], "), `eps_h_fi` dim 1 (= ", dims(eps_h_fi)[1], "), `eps_ca` dim 1 (= ", dims(eps_ca)[1], "), `eps_placental` dim 1 (= ", dims(eps_placental)[1], "), `hq_sw` dim 1 (= ", dims(hq_sw)[1], "), `hq_fi` dim 1 (= ", dims(hq_fi)[1], ").");
if (dims(eps_h_fi)[1] != T) reject("seal_state: dim mismatch — `eps_h_fi` dim 1 (= ", dims(eps_h_fi)[1], ") does not match `T` (= ", T, "), inferred from `eps_birth` dim 1. `T` sizes: `eps_birth` dim 1 (= ", dims(eps_birth)[1], "), `eps_sex` dim 1 (= ", dims(eps_sex)[1], "), `eps_h_sw` dim 1 (= ", dims(eps_h_sw)[1], "), `eps_h_fi` dim 1 (= ", dims(eps_h_fi)[1], "), `eps_ca` dim 1 (= ", dims(eps_ca)[1], "), `eps_placental` dim 1 (= ", dims(eps_placental)[1], "), `hq_sw` dim 1 (= ", dims(hq_sw)[1], "), `hq_fi` dim 1 (= ", dims(hq_fi)[1], ").");
if (dims(eps_ca)[1] != T) reject("seal_state: dim mismatch — `eps_ca` dim 1 (= ", dims(eps_ca)[1], ") does not match `T` (= ", T, "), inferred from `eps_birth` dim 1. `T` sizes: `eps_birth` dim 1 (= ", dims(eps_birth)[1], "), `eps_sex` dim 1 (= ", dims(eps_sex)[1], "), `eps_h_sw` dim 1 (= ", dims(eps_h_sw)[1], "), `eps_h_fi` dim 1 (= ", dims(eps_h_fi)[1], "), `eps_ca` dim 1 (= ", dims(eps_ca)[1], "), `eps_placental` dim 1 (= ", dims(eps_placental)[1], "), `hq_sw` dim 1 (= ", dims(hq_sw)[1], "), `hq_fi` dim 1 (= ", dims(hq_fi)[1], ").");
if (dims(eps_placental)[1] != T) reject("seal_state: dim mismatch — `eps_placental` dim 1 (= ", dims(eps_placental)[1], ") does not match `T` (= ", T, "), inferred from `eps_birth` dim 1. `T` sizes: `eps_birth` dim 1 (= ", dims(eps_birth)[1], "), `eps_sex` dim 1 (= ", dims(eps_sex)[1], "), `eps_h_sw` dim 1 (= ", dims(eps_h_sw)[1], "), `eps_h_fi` dim 1 (= ", dims(eps_h_fi)[1], "), `eps_ca` dim 1 (= ", dims(eps_ca)[1], "), `eps_placental` dim 1 (= ", dims(eps_placental)[1], "), `hq_sw` dim 1 (= ", dims(hq_sw)[1], "), `hq_fi` dim 1 (= ", dims(hq_fi)[1], ").");
if (dims(hq_sw)[1] != T) reject("seal_state: dim mismatch — `hq_sw` dim 1 (= ", dims(hq_sw)[1], ") does not match `T` (= ", T, "), inferred from `eps_birth` dim 1. `T` sizes: `eps_birth` dim 1 (= ", dims(eps_birth)[1], "), `eps_sex` dim 1 (= ", dims(eps_sex)[1], "), `eps_h_sw` dim 1 (= ", dims(eps_h_sw)[1], "), `eps_h_fi` dim 1 (= ", dims(eps_h_fi)[1], "), `eps_ca` dim 1 (= ", dims(eps_ca)[1], "), `eps_placental` dim 1 (= ", dims(eps_placental)[1], "), `hq_sw` dim 1 (= ", dims(hq_sw)[1], "), `hq_fi` dim 1 (= ", dims(hq_fi)[1], ").");
if (dims(hq_fi)[1] != T) reject("seal_state: dim mismatch — `hq_fi` dim 1 (= ", dims(hq_fi)[1], ") does not match `T` (= ", T, "), inferred from `eps_birth` dim 1. `T` sizes: `eps_birth` dim 1 (= ", dims(eps_birth)[1], "), `eps_sex` dim 1 (= ", dims(eps_sex)[1], "), `eps_h_sw` dim 1 (= ", dims(eps_h_sw)[1], "), `eps_h_fi` dim 1 (= ", dims(eps_h_fi)[1], "), `eps_ca` dim 1 (= ", dims(eps_ca)[1], "), `eps_placental` dim 1 (= ", dims(eps_placental)[1], "), `hq_sw` dim 1 (= ", dims(hq_sw)[1], "), `hq_fi` dim 1 (= ", dims(hq_fi)[1], ").");
if (dims(hs_fi)[1] != nd) reject("seal_state: dim mismatch — `hs_fi` dim 1 (= ", dims(hs_fi)[1], ") does not match `nd` (= ", nd, "), inferred from `hs_sw` dim 1. `nd` sizes: `hs_sw` dim 1 (= ", dims(hs_sw)[1], "), `hs_fi` dim 1 (= ", dims(hs_fi)[1], "), `population_init` dim 1 (= ", dims(population_init)[1], ").");
if (dims(population_init)[1] != nd) reject("seal_state: dim mismatch — `population_init` dim 1 (= ", dims(population_init)[1], ") does not match `nd` (= ", nd, "), inferred from `hs_sw` dim 1. `nd` sizes: `hs_sw` dim 1 (= ", dims(hs_sw)[1], "), `hs_fi` dim 1 (= ", dims(hs_fi)[1], "), `population_init` dim 1 (= ", dims(population_init)[1], ").");
int n_demo = nd;
real 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_offset, male_adult_offset);
vector[(2 * n_age)] S_diag = exp((-mu_m));
matrix[n_demo, n_demo] aging = create_aging_matrix(n_demo, n_age);
real min_bbr = (min_bbr_sc * max_bbr);
vector[dims(bbr_logit)[1]] baseline_bbr = (min_bbr + ((max_bbr - min_bbr) * inv_logit(bbr_logit)));
real dd_scaled = birth_rate_at_carrying_capacity(phi_a, mu_m, n_age);
real dd_intercept = compute_density_dependence_intercept(max_bbr, dd_scaled);
real dd_slope = (-log(carrying_capacity));
vector[dims(population_init)[1]] pop_first = initialize_population(
population_init,
pop_init_size,
burn_in,
baseline_bbr[1],
aging,
S_diag,
n_age
);
vector[dims(eps_placental)[1]] pi_s = (report_placental_mean * exp(((-eps_placental) * report_placental_sd)));
vector[dims(eps_ca)[1]] pi_c = (report_ca_mean * exp(((-eps_ca) * report_ca_sd)));
matrix[(3 * n_demo), T] transition_noise_raw = to_matrix(tnoise, (3 * n_demo), T);
return run_state_process(
T,
n_age,
pop_first,
baseline_bbr[1],
sum(pop_first),
baseline_bbr,
dd_intercept,
dd_slope,
aging,
S_diag,
mu_m,
hs_sw,
hs_fi,
hq_sw,
hq_fi,
he_sd_sw,
he_sd_fi,
eps_h_sw,
eps_h_fi,
t_mate,
t_hunt,
eps_birth,
eps_sex,
transition_noise_raw,
pi_s,
pi_c,
prob_of_ca,
ode_init,
ode_times
);
}
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;
}
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);
}
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 neg_binomial_2_int_rng(
int anontok__1,
vector a,
real b
) {
int n = anontok__1;
if((n == 0)) {
array[n] int rv;
return rv;
} else {
return neg_binomial_2_rng(a, b);
}
}
vector aerial_mean(
real mu,
vector poptot,
array[] int years
) {
int na = dims(years)[1];
return (mu * poptot[years]);
}
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 normal_vector_rng(
int anontok__1,
vector a,
vector b
) {
int n = anontok__1;
if((n == 0)) {
vector[n] rv;
return rv;
} else {
return to_vector(normal_rng(a, b));
}
}
vector harvest_mean(
vector total,
array[] int years
) {
int nb = dims(years)[1];
return total[years];
}
vector harvest_sd(
vector total,
array[] int years,
real cv
) {
int nb = dims(years)[1];
return (cv * total[years]);
}
real brm_multinomial_lpmf(
array[, ] int obs,
matrix probs,
array[] int N
) {
int nrow = dims(obs)[1];
int K = dims(obs)[2];
if (dims(probs)[1] != nrow) reject("brm_multinomial_lpmf: dim mismatch — `probs` dim 1 (= ", dims(probs)[1], ") does not match `nrow` (= ", nrow, "), inferred from `obs` dim 1. `nrow` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `probs` dim 1 (= ", dims(probs)[1], "), `N` dim 1 (= ", dims(N)[1], ").");
if (dims(N)[1] != nrow) reject("brm_multinomial_lpmf: dim mismatch — `N` dim 1 (= ", dims(N)[1], ") does not match `nrow` (= ", nrow, "), inferred from `obs` dim 1. `nrow` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `probs` dim 1 (= ", dims(probs)[1], "), `N` dim 1 (= ", dims(N)[1], ").");
if (dims(probs)[2] != K) reject("brm_multinomial_lpmf: dim mismatch — `probs` dim 2 (= ", dims(probs)[2], ") does not match `K` (= ", K, "), inferred from `obs` dim 2. `K` sizes: `obs` dim 2 (= ", dims(obs)[2], "), `probs` dim 2 (= ", dims(probs)[2], ").");
real ll = 0.0;
for(i in 1:nrow) {
ll = (ll + multinomial_lpmf(obs[i] | to_vector(probs[i, :]), N[i]));
}
return ll;
}
real multinomial_lpmf(
array[] int obs,
vector probs,
int N
) {
int K = dims(obs)[1];
if (dims(probs)[1] != K) reject("multinomial_lpmf: dim mismatch — `probs` dim 1 (= ", dims(probs)[1], ") does not match `K` (= ", K, "), inferred from `obs` dim 1. `K` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `probs` dim 1 (= ", dims(probs)[1], ").");
if((sum(obs) != N)) {
reject("multinomial: explicit N must equal sum(obs)");
}
return multinomial_lpmf(obs | probs);
}
vector brm_multinomial_lpmfs(
array[, ] int obs,
matrix probs,
array[] int N
) {
int nrow = dims(obs)[1];
int K = dims(obs)[2];
if (dims(probs)[1] != nrow) reject("brm_multinomial_lpmfs: dim mismatch — `probs` dim 1 (= ", dims(probs)[1], ") does not match `nrow` (= ", nrow, "), inferred from `obs` dim 1. `nrow` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `probs` dim 1 (= ", dims(probs)[1], "), `N` dim 1 (= ", dims(N)[1], ").");
if (dims(N)[1] != nrow) reject("brm_multinomial_lpmfs: dim mismatch — `N` dim 1 (= ", dims(N)[1], ") does not match `nrow` (= ", nrow, "), inferred from `obs` dim 1. `nrow` sizes: `obs` dim 1 (= ", dims(obs)[1], "), `probs` dim 1 (= ", dims(probs)[1], "), `N` dim 1 (= ", dims(N)[1], ").");
if (dims(probs)[2] != K) reject("brm_multinomial_lpmfs: dim mismatch — `probs` dim 2 (= ", dims(probs)[2], ") does not match `K` (= ", K, "), inferred from `obs` dim 2. `K` sizes: `obs` dim 2 (= ", dims(obs)[2], "), `probs` dim 2 (= ", dims(probs)[2], ").");
vector[nrow] lls;
for(i in 1:nrow) {
lls[i] = multinomial_lpmf(obs[i] | to_vector(probs[i, :]), N[i]);
}
return lls;
}
array[, ] int brm_multinomial_int_rng(
tuple(int, int) anontok__1,
matrix probs,
array[] int N
) {
int nrow = anontok__1.1;
int K = anontok__1.2;
if (dims(probs)[1] != nrow) reject("brm_multinomial_rng: dim mismatch — `probs` dim 1 (= ", dims(probs)[1], ") does not match `nrow` (= ", nrow, "), inferred from `anontok__1` dim 1. `nrow` sizes: `anontok__1` dim 1 (= ", anontok__1.1, "), `probs` dim 1 (= ", dims(probs)[1], "), `N` dim 1 (= ", dims(N)[1], ").");
if (dims(N)[1] != nrow) reject("brm_multinomial_rng: dim mismatch — `N` dim 1 (= ", dims(N)[1], ") does not match `nrow` (= ", nrow, "), inferred from `anontok__1` dim 1. `nrow` sizes: `anontok__1` dim 1 (= ", anontok__1.1, "), `probs` dim 1 (= ", dims(probs)[1], "), `N` dim 1 (= ", dims(N)[1], ").");
if (dims(probs)[2] != K) reject("brm_multinomial_rng: dim mismatch — `probs` dim 2 (= ", dims(probs)[2], ") does not match `K` (= ", K, "), inferred from `anontok__1` dim 2. `K` sizes: `anontok__1` dim 2 (= ", anontok__1.2, "), `probs` dim 2 (= ", dims(probs)[2], ").");
array[nrow, K] int rv;
for(i in 1:nrow) {
rv[i, :] = multinomial_rng(to_vector(probs[i, :]), N[i]);
}
return rv;
}
matrix comp_hunted(
matrix hunted,
vector total,
array[] int years
) {
int K = dims(hunted)[1];
int T = dims(hunted)[2];
int nc = dims(years)[1];
if (dims(total)[1] != T) reject("comp_hunted: dim mismatch — `total` dim 1 (= ", dims(total)[1], ") does not match `T` (= ", T, "), inferred from `hunted` dim 2. `T` sizes: `hunted` dim 2 (= ", dims(hunted)[2], "), `total` dim 1 (= ", dims(total)[1], ").");
matrix[nc, K] p;
for(i in 1:nc) {
p[i, :] = ((hunted[:, years[i]] ./ total[years[i]])');
}
return p;
}
matrix comp_bycatch(
matrix bycatch,
vector bias,
array[] int years
) {
int K = dims(bycatch)[1];
int nc = dims(years)[1];
if (dims(bias)[1] != K) reject("comp_bycatch: dim mismatch — `bias` dim 1 (= ", dims(bias)[1], ") does not match `K` (= ", K, "), inferred from `bycatch` dim 1. `K` sizes: `bycatch` dim 1 (= ", dims(bycatch)[1], "), `bias` dim 1 (= ", dims(bias)[1], ").");
vector[dims(bias)[1]] w = exp(bias);
matrix[nc, K] p;
for(i in 1:nc) {
p[i, :] = (((w .* bycatch[:, years[i]]) ./ dot_product(w, bycatch[:, years[i]]))');
}
return p;
}
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 binomial_int_rng(
int anontok__1,
array[] int N,
vector p
) {
int n = anontok__1;
if (dims(N)[1] != n) reject("binomial_rng: dim mismatch — `N` dim 1 (= ", dims(N)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `N` dim 1 (= ", dims(N)[1], ").");
if((n == 0)) {
array[n] int rv;
return rv;
} else {
return binomial_rng(N, p);
}
}
matrix comp_repro(
matrix rp,
array[] int years
) {
int nc = dims(years)[1];
matrix[nc, 4] p;
for(i in 1:nc) {
p[i, :] = (rp[:, years[i]]');
}
return p;
}
}
data {
int herring_index_2_n;
int herring_index_1_n;
vector[herring_index_1_n] herring_index_1;
vector[herring_index_2_n] herring_index_2;
int n_year;
int year_idx_n;
array[year_idx_n] int year_idx;
int n_demo;
int demo_idx_n;
array[demo_idx_n] int demo_idx;
int n_noise_cell;
int noise_cell_idx_n;
array[noise_cell_idx_n] int noise_cell_idx;
int hunting_quota_finland_n;
int population_init_n;
vector[population_init_n] population_init;
int hunting_quota_sweden_n;
array[hunting_quota_sweden_n] int hunting_quota_sweden;
array[hunting_quota_finland_n] int hunting_quota_finland;
real t_mate_to_preg;
real t_birth_to_end_hunt;
int population_burn_in;
int n_age;
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_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_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_pregnancy_count_n;
array[obs_pregnancy_count_n] int obs_pregnancy_count;
int pregnancy_sample_size_n;
array[pregnancy_sample_size_n] int pregnancy_sample_size;
int pregnancy_count_year_n;
array[pregnancy_count_year_n] int pregnancy_count_year;
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[herring_index_2_n, (2 + 1)] X_bbr_logit = hcat(rep_vector(1.0, num_elements(herring_index_1)), herring_index_1, herring_index_2);
int pop_bbr_logit_n_covariates = (2 + 1);
}
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<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, 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;
real<lower=0.0, upper=1.0> aerial_mu;
real<lower=0.0> phi_aerial;
real<lower=0.0> harvest_bag_cv;
vector[pop_bbr_logit_n_covariates] pop_bbr_logit_beta_pop;
real r_eps_birth_year_log_scale;
vector[n_year] r_eps_birth_year_xi;
real r_eps_sex_year_log_scale;
vector[n_year] r_eps_sex_year_xi;
real r_eps_h_sw_year_log_scale;
vector[n_year] r_eps_h_sw_year_xi;
real r_eps_h_fi_year_log_scale;
vector[n_year] r_eps_h_fi_year_xi;
real r_eps_ca_year_log_scale;
vector[n_year] r_eps_ca_year_xi;
real r_eps_placental_year_log_scale;
vector[n_year] r_eps_placental_year_xi;
real r_hs_sw_demo_log_scale;
vector[n_demo] r_hs_sw_demo_xi;
real r_hs_fi_demo_log_scale;
vector[n_demo] r_hs_fi_demo_xi;
real r_bycatch_bias_demo_log_scale;
vector[n_demo] r_bycatch_bias_demo_xi;
real r_tnoise_noise_cell_log_scale;
vector[n_noise_cell] r_tnoise_noise_cell_xi;
}
transformed parameters {
vector[herring_index_2_n] pop_bbr_logit = (X_bbr_logit * pop_bbr_logit_beta_pop);
vector[herring_index_2_n] bbr_logit = pop_bbr_logit;
vector[year_idx_n] r_eps_birth_year = (exp(r_eps_birth_year_log_scale) * r_eps_birth_year_xi[year_idx]);
vector[year_idx_n] eps_birth = r_eps_birth_year;
vector[year_idx_n] r_eps_sex_year = (exp(r_eps_sex_year_log_scale) * r_eps_sex_year_xi[year_idx]);
vector[year_idx_n] eps_sex = r_eps_sex_year;
vector[year_idx_n] r_eps_h_sw_year = (exp(r_eps_h_sw_year_log_scale) * r_eps_h_sw_year_xi[year_idx]);
vector[year_idx_n] eps_h_sw = r_eps_h_sw_year;
vector[year_idx_n] r_eps_h_fi_year = (exp(r_eps_h_fi_year_log_scale) * r_eps_h_fi_year_xi[year_idx]);
vector[year_idx_n] eps_h_fi = r_eps_h_fi_year;
vector[year_idx_n] r_eps_ca_year = (exp(r_eps_ca_year_log_scale) * r_eps_ca_year_xi[year_idx]);
vector[year_idx_n] eps_ca = r_eps_ca_year;
vector[year_idx_n] r_eps_placental_year = (exp(r_eps_placental_year_log_scale) * r_eps_placental_year_xi[year_idx]);
vector[year_idx_n] eps_placental = r_eps_placental_year;
vector[demo_idx_n] r_hs_sw_demo = (exp(r_hs_sw_demo_log_scale) * r_hs_sw_demo_xi[demo_idx]);
vector[demo_idx_n] hs_sw = r_hs_sw_demo;
vector[demo_idx_n] r_hs_fi_demo = (exp(r_hs_fi_demo_log_scale) * r_hs_fi_demo_xi[demo_idx]);
vector[demo_idx_n] hs_fi = r_hs_fi_demo;
vector[demo_idx_n] r_bycatch_bias_demo = (exp(r_bycatch_bias_demo_log_scale) * r_bycatch_bias_demo_xi[demo_idx]);
vector[demo_idx_n] bycatch_bias = r_bycatch_bias_demo;
vector[noise_cell_idx_n] r_tnoise_noise_cell = (exp(r_tnoise_noise_cell_log_scale) * r_tnoise_noise_cell_xi[noise_cell_idx]);
vector[noise_cell_idx_n] tnoise = r_tnoise_noise_cell;
tuple(
vector[hunting_quota_finland_n],
vector[hunting_quota_finland_n],
vector[hunting_quota_finland_n],
matrix[demo_idx_n, hunting_quota_finland_n],
matrix[demo_idx_n, hunting_quota_finland_n],
matrix[demo_idx_n, hunting_quota_finland_n],
vector[hunting_quota_finland_n],
vector[hunting_quota_finland_n],
matrix[4, hunting_quota_finland_n]
) state = seal_state(
bbr_logit,
eps_birth,
eps_sex,
eps_h_sw,
eps_h_fi,
eps_ca,
eps_placental,
hs_sw,
hs_fi,
tnoise,
phi_a_sc,
phi_sc,
survival_shape,
male_pup_survival_offset,
male_adult_survival_offset,
carrying_capacity,
max_baseline_birth_rate,
min_baseline_birth_rate_sc,
hunting_effort_sd_sweden,
hunting_effort_sd_finland,
population_init_size,
report_ca_mean,
report_placental_mean,
report_ca_sd,
report_placental_sd,
prob_of_ca,
population_init,
hunting_quota_sweden,
hunting_quota_finland,
t_mate_to_preg,
t_birth_to_end_hunt,
population_burn_in,
n_age,
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);
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);
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);
aerial_mu ~ uniform(0.0, 1.0);
phi_aerial ~ lognormal(0.0, 1.0);
harvest_bag_cv ~ lognormal(0.0, 1.0);
pop_bbr_logit_beta_pop ~ std_normal();
r_eps_birth_year_log_scale ~ std_normal();
r_eps_birth_year_xi ~ std_normal();
r_eps_sex_year_log_scale ~ std_normal();
r_eps_sex_year_xi ~ std_normal();
r_eps_h_sw_year_log_scale ~ std_normal();
r_eps_h_sw_year_xi ~ std_normal();
r_eps_h_fi_year_log_scale ~ std_normal();
r_eps_h_fi_year_xi ~ std_normal();
r_eps_ca_year_log_scale ~ std_normal();
r_eps_ca_year_xi ~ std_normal();
r_eps_placental_year_log_scale ~ std_normal();
r_eps_placental_year_xi ~ std_normal();
r_hs_sw_demo_log_scale ~ std_normal();
r_hs_sw_demo_xi ~ std_normal();
r_hs_fi_demo_log_scale ~ std_normal();
r_hs_fi_demo_xi ~ std_normal();
r_bycatch_bias_demo_log_scale ~ std_normal();
r_bycatch_bias_demo_xi ~ std_normal();
r_tnoise_noise_cell_log_scale ~ std_normal();
r_tnoise_noise_cell_xi ~ std_normal();
obs_aerial_count ~ neg_binomial_2(aerial_mean(aerial_mu, state.3, aerial_year), phi_aerial);
obs_hunting_bag_sweden ~ normal(
harvest_mean(state.7, hunting_bag_year_sweden),
harvest_sd(state.7, hunting_bag_year_sweden, harvest_bag_cv)
);
obs_hunting_bag_finland ~ normal(
harvest_mean(state.8, hunting_bag_year_finland),
harvest_sd(state.8, hunting_bag_year_finland, harvest_bag_cv)
);
obs_hunting_comp_sweden ~ brm_multinomial(
comp_hunted(state.4, state.7, hunting_comp_year_sweden),
hunting_comp_sample_size_sweden
);
obs_hunting_comp_finland ~ brm_multinomial(
comp_hunted(state.5, state.8, hunting_comp_year_finland),
hunting_comp_sample_size_finland
);
obs_bycatch_comp ~ brm_multinomial(comp_bycatch(state.6, bycatch_bias, bycatch_comp_year), bycatch_comp_sample_size);
obs_pregnancy_count ~ binomial(pregnancy_sample_size, state.2[pregnancy_count_year]);
obs_reproductive_signs_finland ~ brm_multinomial(comp_repro(state.9, reproductive_signs_year), reproductive_signs_sample_size);
}
generated quantities {
vector[obs_aerial_count_n] obs_aerial_count_likelihood = neg_binomial_2_lpmfs(obs_aerial_count, aerial_mean(aerial_mu, state.3, aerial_year), phi_aerial);
array[obs_aerial_count_n] int obs_aerial_count_gen = neg_binomial_2_int_rng(obs_aerial_count_n, aerial_mean(aerial_mu, state.3, aerial_year), phi_aerial);
vector[obs_hunting_bag_sweden_n] obs_hunting_bag_sweden_likelihood = normal_lpdfs(
obs_hunting_bag_sweden,
harvest_mean(state.7, hunting_bag_year_sweden),
harvest_sd(state.7, hunting_bag_year_sweden, harvest_bag_cv)
);
vector[obs_hunting_bag_sweden_n] obs_hunting_bag_sweden_gen = normal_vector_rng(
obs_hunting_bag_sweden_n,
harvest_mean(state.7, hunting_bag_year_sweden),
harvest_sd(state.7, hunting_bag_year_sweden, harvest_bag_cv)
);
vector[obs_hunting_bag_finland_n] obs_hunting_bag_finland_likelihood = normal_lpdfs(
obs_hunting_bag_finland,
harvest_mean(state.8, hunting_bag_year_finland),
harvest_sd(state.8, hunting_bag_year_finland, harvest_bag_cv)
);
vector[obs_hunting_bag_finland_n] obs_hunting_bag_finland_gen = normal_vector_rng(
obs_hunting_bag_finland_n,
harvest_mean(state.8, hunting_bag_year_finland),
harvest_sd(state.8, hunting_bag_year_finland, harvest_bag_cv)
);
vector[hunting_comp_sample_size_sweden_n] obs_hunting_comp_sweden_likelihood = brm_multinomial_lpmfs(
obs_hunting_comp_sweden,
comp_hunted(state.4, state.7, hunting_comp_year_sweden),
hunting_comp_sample_size_sweden
);
array[hunting_comp_sample_size_sweden_n, demo_idx_n] int obs_hunting_comp_sweden_gen = brm_multinomial_int_rng(
(obs_hunting_comp_sweden_m, obs_hunting_comp_sweden_n),
comp_hunted(state.4, state.7, hunting_comp_year_sweden),
hunting_comp_sample_size_sweden
);
vector[hunting_comp_sample_size_finland_n] obs_hunting_comp_finland_likelihood = brm_multinomial_lpmfs(
obs_hunting_comp_finland,
comp_hunted(state.5, state.8, hunting_comp_year_finland),
hunting_comp_sample_size_finland
);
array[hunting_comp_sample_size_finland_n, demo_idx_n] int obs_hunting_comp_finland_gen = brm_multinomial_int_rng(
(obs_hunting_comp_finland_m, obs_hunting_comp_finland_n),
comp_hunted(state.5, state.8, hunting_comp_year_finland),
hunting_comp_sample_size_finland
);
vector[bycatch_comp_sample_size_n] obs_bycatch_comp_likelihood = brm_multinomial_lpmfs(
obs_bycatch_comp,
comp_bycatch(state.6, bycatch_bias, bycatch_comp_year),
bycatch_comp_sample_size
);
array[bycatch_comp_sample_size_n, demo_idx_n] int obs_bycatch_comp_gen = brm_multinomial_int_rng(
(obs_bycatch_comp_m, obs_bycatch_comp_n),
comp_bycatch(state.6, bycatch_bias, bycatch_comp_year),
bycatch_comp_sample_size
);
vector[obs_pregnancy_count_n] obs_pregnancy_count_likelihood = binomial_lpmfs(obs_pregnancy_count, pregnancy_sample_size, state.2[pregnancy_count_year]);
array[pregnancy_sample_size_n] int obs_pregnancy_count_gen = binomial_int_rng(obs_pregnancy_count_n, pregnancy_sample_size, state.2[pregnancy_count_year]);
vector[reproductive_signs_sample_size_n] obs_reproductive_signs_finland_likelihood = brm_multinomial_lpmfs(
obs_reproductive_signs_finland,
comp_repro(state.9, reproductive_signs_year),
reproductive_signs_sample_size
);
array[reproductive_signs_sample_size_n, 4] int obs_reproductive_signs_finland_gen = brm_multinomial_int_rng(
(obs_reproductive_signs_finland_m, obs_reproductive_signs_finland_n),
comp_repro(state.9, reproductive_signs_year),
reproductive_signs_sample_size
);
}julia
Turing unsupported for this BRM example
Turing backend: response and predictor row counts differThe Turing pane is intentionally retained even though this structural model is outside the current Turing executor; its build-time construction error is part of the comparison rather than being hidden.
Provenance and scope
Reproduction: research/seal/grey_seal_brm.jl, gated by test/seal_brm.jl. The 13 state-process UDFs and run_state_process are verbatim from the SlicTranspiler research model; the fixture is its example dataset (3 age classes × 2 sexes × 3 state years). Verified at the same level as the @slic reference — transpile + stanc + compiles. Like the reference it is not sampled: the mechanistic scan's simplex constraints make a naive init degenerate for both forms (the @slic reference is non-finite at the origin too), so the deliverable is the model. The StanBlocks-native case study carries the same model in its own idiom.
Run the reproduction
After bootstrapping the repository's test environment:
sh
julia --startup-file=no --project=test test/seal_brm.jl