It’s curious that the summary statistics are (extremely) locally sufficient, though. When p is set to precisely the true value, everything works totally fine. But I also tried setting p \sim N(6, 1), and even this caused the model to fail. I would have expected a tight prior around the true value to at least give something, but it didn’t.
I struggle with the notion of identification in Bayesian models. I am often stuck in a frequentist “you can’t have more parameters than observations” mindset, but this is not always the case, nor is it even the case here: the model estimates 502 parameters perfectly from 500 data points, given the true value of p. Things are evidently much more subtle in MCMClandia than in the People’s Republic of Matrix Inversion .
The functional form of the dependence between p and \gamma in the pair plot is also tantalizing. I wonder if one could do some algebra on the model to discover precisely what is going on here. A mystery for another day, I suppose.