I need some help troubleshooting divergencies (534 out of 30,000 after warmup iterations, 4 chains). It seems that all parametrization guidelines were followed. pairs() plot shows this but I have no clue what to do to remove those divergencies except making adapt_delta 0.9999.

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 m0 as normal(0, 1) and define mu[1] = 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[1] around 100. Then you plug that into an exponential. How big are your y values?

Also, the function weibull should be vectorized so that you can vectorize the sampling for y—it’ll be a lot faster.