I’m struggling to specify my model. I’m new to hierarchical models, so may be missing some conceptual clarity too.

I’m modeling the output of an automated classifier for detecting a species in an audio clip as a 2-component gaussian mixture. A component with higher mean for the positives (clips that contain the species of interest), and a component with a lower mean for the negatives (clips that don’t contain the species of interest). I want to estimate the mean number of animals per site. I assume the number of animals ( N_i ) at a site is poisson distributed, each animal has a fixed probability (p) of vocalizing. So the probability of an animal being present in a clip is given by 1-(1-p)^{N_i} . It’s basically a Royle-Nichols model, but with observed data (y) that are ‘scores’ drawn from the gaussian mixture rather than true 1/0 detections.

The model looks something like:

N_i \sim Poisson(\lambda)

n\_positives_i \sim Binomial(n\_clips_i, (1-(1-p)^{N_i})

\theta_i = \frac{n\_positives_i}{n\_clips_i}

y_{ij} \sim \theta_i N(\mu^{1}_i, \sigma^{1}_i) + (1-\theta_i) N(\mu^{2}_i, \sigma^{2}_i)

Here’s my code so far, which is 1. Compiling (but wrong - looking for help in specifying the model) and 2. Very slow

```
data {
int<lower=1> M; // number of sites
int<lower=1> J; // number of scores
array[M,J] real y; // response variable
int<lower=0> N_max; // Max possible number of animals at a site to marginalize up to
}
parameters {
real mu1; // gaussian params. Could use ordered[2] instead.
real mu2;
array[2] real<lower=0> sigma;
real<lower=0,upper=1> p; // per-capita p of detection
real lambda; // average abundance
}
transformed parameters {
array[N_max] real<lower=0, upper=1> p_det; // probability of detection given n animals
for (n in 1:N_max){
p_det[n] = 1 - (1 - p)^n;
}
model {
// priors
mu1 ~ normal(-1, 2);
mu2 ~ normal(5, 3);
sigma ~ normal(0, 3);
p ~ beta(1, 1);
lambda ~ uniform(0, 15);
for (site in 1:M){
for (j in 1:J){ //scores
for (i in 0:N_max){ // possible number of animals
if (i == 0){. // All scores are negatives
target += poisson_lpmf(i | lambda) + normal_lpdf(y[site, j] | mu1, sigma[1]);
} else {
target += poisson_lpmf(i | lambda);
target += log_mix(p_det[i],
normal_lpdf(y[site, j] | mu1, sigma[1]),
normal_lpdf(y[site, j] | mu2, sigma[2]));
}
}
}
}
}
```

Any help is much appreciated.