Warning: Please double check everything I say because on the one hand I feel like I am on to something. On the other hand, I have no special skills in any of this type of modeling and I am pretty sure I messed up somewhere.
I am not really sure how this code can generate N < max_y if N = k and k is between max_y and K. Are you sure that you get N < max_y.
I think the code that you have is accounting for the fact that the poisson is trunctated at K and you cannot have any probability above K (- poisson_lcdf(K[i])). You have to make the same adjustment for the fact that the poisson is also truncated below at max_iter.
However - and this is the most speculative part -, the loglikelihood that is in the model is not only truncated. It’s also perturbed by the binomial. So to take that into account I think this is how you generate N equivalent to your model. The code is in R because it’s easier for me to proto type in R.
For each observation and mcmc draw, start with generating p according to the model. max_y, y, and K are data. lambda is an estimated parameter. So you have a number for each. The first likelihood is the binomial, the second is the poisson in your model (I hope I got the relation between y and max_y right here). Multiplying the likelihoods (or summing the log_likelihoods as in the stan model) gives you the perturbed poisson. Than standardize the likelihood to get the probabilities for each k (This is the step I am most unsure about). Sample N from the probabilities of all ks. If you don’t really need exact simulated N but you are more interested in the distribution of the N for each observation I think you could also calculate the expected value for N (meanN). Remember that is the expected value for each observation and each mcmc draw.
p <- .75
max_y <- 50
lambda <- 40
y <- 40
K <- 10
ks <- max_y + (0:K)
lik_binom <- dbinom(x = y, size = ks, p = p)
lik_poisson <- dpois(x = ks, lambda = lambda)
lik = lik_binom * lik_poisson
probs <- lik/sum(lik)
N <- sample(x = ks, size = 1, prob = probs)
meanN = sum(ks * probs)