functions{ real integrand(real log_lambda, real notused, array[] real theta, array[] real log_lambda, array[] int yi) { real alpha = theta[1] ; real beta = theta [2] ; real sigma= theta [3] ; real log_lambda_i= theta [4] ; return exp(normal_lpdf(log_lambda|alpha+beta*log_lambda_i, sigma) + poisson_log_lpmf(yi|log_lambda)); } } data { int N; array[N] int y; } parameters { real alpha; real beta; real sigma; vector[N] log_lambda; } model { alpha ~ normal(0, 1); beta ~ normal(0, 1); sigma ~ exponential(1); //Latent true expected abundance log_lambda[1] ~ normal(0,5); for(n in 2:N){ log_lambda[n] ~ normal(alpha + beta * log_lambda[n-1], sigma); } for(n in 1:N){ target += poisson_log_lpmf(y[n] | log_lambda[n]); } } generated quantities{ vector[N] log_lik; vector[N] log_lambda_pred; array[N] int y_pred; log_lambda_pred[1] = log_lambda[1]; for(t in 2:N){ log_lambda_pred[t] = normal_rng(alpha + beta*log_lambda[t-1], sigma); } for(n in 1:N){ //log_lik[n] = poisson_log_lpmf(y[n] | log_lambda[n]); log_lik[n] = log(integrate_1d(integrand, negative_infinity(), positive_infinity(), {alpha, beta,sigma, log_lambda[n-1]}, log_lambda[n], {y[n]},0.00001)); if(log_lambda_pred[n]>20){ y_pred[n] = poisson_log_rng(20); }else{ y_pred[n] = poisson_log_rng(log_lambda_pred[n]); } } }