I worked on this awhile but didn’t get too far. I one issue might be with the if(time > tlag) conditional.
I don’t think the derivatives of the Weibull function are always zero at zero (going from pictures here: https://en.wikipedia.org/wiki/Weibull_distribution). Functions we put in HMC have to be continuous and have first derivatives, otherwise it can mess up the sampler, so it might not be totally okay to put a zero on the left side of the Weibull.
If we’re shifting the timescale, that’s a problem. It’s not a problem for the other parameters cause evaluations of that function either stay on one side or the other.
The most interesting thing I got was moving the tlag away from the shape/scale parameters:
y = m0 * (1 - exp(-(time / scale) ^ shape) - tlag_not_quite);
If it sits on the outside like this, the parameters seem to sample a lot better than if it sits inside with the time (you can verify this with a smaller model that only works on one lot if you want).
If you parameterize this way, you could do algebra to compute effectively what the t0 would be to reproduce the effective argument to the exp, but I dunno how far that gets you. The priors on tlag_not_quite wouldn’t be that easy to interpret, but if you’re just wanting to fit the curves, might be worth a try.