Your marginalized model is working for me. However, I don’t think your data generation matches your models.
This produces a 10 x 100 array of true occupation data (pres_vec) and a 10 x 100 array of detection data (dets_vec). Based on the occupation data, it would seem that the animal is truly moving in and out of each site, as if you have one sensor at each site and you’re collecting a time series.
However, your models (and, I think, the two resources you linked) suggest that they’re either in the whole time or out, as if you’re collecting data from 100 sensors per site at one time point. For example
specifies only one true occupancy datum per site.
Assuming that model is correct (and I apologize, as I’m not familiar with these models so I might be making some naïve assumptions), then you probably want something more like this, where you first simulate a single occupancy observation per site and then simulate your detection data per visit based on the occupancy datum.
pres_vec = rand(Bernoulli(ψ), N_sites)
dets_vec = [rand.(Bernoulli.(pres * p), N_visits) for pres in pres_vec]