Here is a toy version of a population model

```
### STAN CODE - simplePopMod.stan ######################
data {
int popInit;
real propYInit;
real preg;
real Ma;
real Mj;
real sigma;
int years;
// real ObsA[years];
}
parameters {
real<lower=0, upper=1> Mjuv;
// real pObs;
}
transformed parameters {
vector<lower=0>[2] pop[years] ;
real propY[years];
pop[1,1] = popInit * propYInit;
pop[1,2] = popInit * (1 - propYInit);
for (year in 1:(years-1)){
pop[year+1,1] = pop[year,2] * preg * (1 - Mjuv) ;
pop[year+1,2] = (pop[year,1] + pop[year,2]) * (1 - Ma);
}
}
model {
Mjuv ~ normal(Mj, sigma);
// pObs ~ uniform(0,1);
// for (year in 1:years){
// ObsA[year] ~ binomial(pop[year,2], pObs);
// }
}
```

```
### R CODE #################################################
library(rstan)
fitmodel = "simplePopMod.stan"
stan.data <- list(years = 100, popInit = 100, propYInit = 0.3, preg = 0.4, Ma = 0.05, Mj = 0.6, sigma = 0.05)#, ObsA = ObsA)
init_fun <- function() {list(
Mj=runif(1, 0.5, 0.9)
)}
params <- c("Mjuv", "pop")
nsamples <- 1000
nburnin <- 1000
ncore <- 1
out <- stan(
file = fitmodel,
data = stan.data,
pars = params,
init= init_fun,
chains = ncore,
warmup = nburnin,
iter = nburnin+nsamples,
cores = ncore,
refresh = 100,
)
```

Now imagine you have for each year observation data of the number of adults (ObsA). You know that you have only a certain probability (pObs) to observe the truth (pop[year, 2]) and you know the relationship between the real number and the observed number is the result of a binomial process.

So I think I can write the following:

```
pObs ~ uniform(0,1);
for (year in 1:years){
ObsA[year] ~ binomial(pop[year,2], pObs);
}
```

However, in Stan, pop[year,2] is a â€śrealâ€ť and the binomial function only accept integers (and you cannot convert a real into integer using the round function).

So if someone can provide a way to fit such relationship in Stan, that would be great !