I don’t see any immediate problem with the code. Do the divergences go away with adapt_delta = 0.999? If so, that seems like the solution :-)
You might want to reparameterize
normal(0, 1) and define
mu = 100 + 5 * m0;.
You can try putting stronger priors on things—looks like the problem you’re having is that you have this exponential
weibull function calculating the mean of your normal and vaues of
mu around 100. Then you plug that into an exponential. How big are your
Also, the function
weibull should be vectorized so that you can vectorize the sampling for
y—it’ll be a lot faster.