Ok, I translated the DGP to Stan and follow the instructions here:

```
functions {
array [,] int simulate_introduction_rng(int N, int B_total, array [] int sampling_vector, real b0_sim, real b1_sim){
array[N,7] int out_mat;
vector[N] rate_sim;
vector[N] prob_invasive;
array[N] int samples;
array[N] int yearly_introduced;
array[N] int detected_i_species;
array[N] int undetected_i_before;
array[N] int undetected_i_after;
array[N] int detected_i_so_far;
array[N] int detected_n_species;
array[N] int undetected_n_before;
array[N] int undetected_n_after;
array[N] int detected_n_so_far;
array[N] int t;
for (i in 1:N) {
t[i] = i;
rate_sim[i] = exp(b0_sim + b1_sim * i);
}
yearly_introduced = poisson_rng(rate_sim);
detected_i_species = rep_array(0, N);
undetected_i_before = rep_array(0, N);
undetected_i_after = rep_array(0, N);
detected_i_so_far = rep_array(0, N);
detected_n_species = rep_array(0, N);
undetected_n_before = rep_array(0, N);
undetected_n_after = rep_array(0, N);
detected_n_so_far = rep_array(0, N);
undetected_i_before[1] = yearly_introduced[1];
undetected_n_before[1] = B_total;
prob_invasive[1] = 1.0 * undetected_i_before[1] / (undetected_n_before[1] + undetected_i_before[1]);
samples[1] = binomial_rng(sampling_vector[1], prob_invasive[1]);
detected_i_species[1] = samples[1];
detected_n_species[1] = sampling_vector[1] - samples[1];
undetected_i_after[1] = undetected_i_before[1] - detected_i_species[1];
undetected_n_after[1] = undetected_n_before[1] - detected_n_species[1];
detected_i_so_far[1] = cumulative_sum(detected_i_species)[1];
detected_n_so_far[1] = cumulative_sum(detected_n_species)[1];
for (i in 2:N){
undetected_i_before[i] = yearly_introduced[i] + undetected_i_after[i-1];
undetected_n_before[i] = undetected_n_after[i-1];
prob_invasive[i] = 1.0 * undetected_i_before[i] / (undetected_n_before[i] + undetected_i_before[i]);
samples[i] = binomial_rng(sampling_vector[i], prob_invasive[i]);
detected_i_species[i] = min(samples[i], undetected_i_before[i]);
detected_n_species[i] = min((sampling_vector[i] - samples[i]), undetected_n_before[i]);
undetected_i_after[i] = undetected_i_before[i] - detected_i_species[i];
undetected_n_after[i] = undetected_n_before[i] - detected_n_species[i];
detected_i_so_far[i] = cumulative_sum(detected_i_species)[i];
detected_n_so_far[i] = cumulative_sum(detected_n_species)[i];
}
out_mat[,1] = t;
out_mat[,2] = yearly_introduced;
out_mat[,3] = detected_i_species;
out_mat[,4] = detected_i_so_far;
out_mat[,5] = undetected_n_before;
out_mat[,6] = detected_n_species;
out_mat[,7] = detected_n_so_far;
return out_mat;
}
}
data {
int <lower = 1> N; // number of rows in the data
int <lower = 1> B_total; // assumed number of ind in group B
array[N] int <lower = 1> sampling_vector; // number of
real b0_mu; // prior for betas
real b1_mu; // prior for betas
}
transformed data {
real b0_sim = normal_rng(b0_mu, 0.1); // params for data creation
real b1_sim = normal_rng(b1_mu, 0.001); // params for data creation
array[N,7] int <lower = 0> sampling_mat;
sampling_mat = simulate_introduction_rng(N, B_total, sampling_vector, b0_sim, b1_sim);
}
parameters {
real b0;
real b1;
}
transformed parameters {
vector<lower = 0>[N] unrecorded_invasives;
vector[N] rate;
for (i in 1:N){
rate[i] = exp(b0 + b1 * i);
unrecorded_invasives[i] = cumulative_sum(rate)[i] - sampling_mat[i,4];
}
}
model{
for (i in 1:N){
sampling_mat[i,3] ~ beta_binomial(sampling_vector[i], unrecorded_invasives[i], sampling_mat[i,5]);
}
//priors
b0 ~ normal(b0_mu, 0.1);
b1 ~ normal(b1_mu , 0.001);
}
generated quantities {
array[2] int <lower = 0, upper = 1> lt_sim
= { b0 < b0_sim , b1 < b1_sim };
}
```

I’m left with a rank for this specific run. So, the next step would be something like this:

```
data_scb <- list(
N = 150,
sampling_vector = rpois(150, 10),
B_total = 1000,
b0_mu = 1,
b1_mu = 0.01
)
n_trials = 30
ranks <- matrix(NA, n_trials , 2)
for (n in 1:n_trials ){
samples <- scb_model$sample(data = data_scb,
chains = 3,
parallel_chains = 3,
iter_warmup = 200,
iter_sampling = 1000, thin = 2,
show_messages = F, refresh = 0)
specific_rank <- as_draws_rvars(samples$draws(variables = "lt_sim"))
specific_rank <- draws_of(specific_rank $lt_sim)
ranks[n,] <- colSums(specific_rank )
}
```

If I understand, then I am hoping to get a uniform distribution for the ranks, between 0 and number of iterations. That indicates that the model “works”.

But still, it seems like it doesn’t answer my question, does it?

But the more I think about it, these simulations should assume my prior is correct…