Consider the traces on the image. These are 8 NUTS chains with 90000 samples each, first third being the tuning samples.
6 out of those 8 have converged (to the correct values), two have not. But what is of interest is the difference between the converged chains: the last two chains look really nice, but chains 0,2,3,4 are random walking(!).
Now what I would like to do is to take the scale matrix and integration time parameters used in chains 6 or 7 and use those in ordinary HMC without tuning.
How do I do that?
I think this is a reason to concern about your model when you see chains like this. It is likely there are some hidden problems in your model parameterization that makes the sampler having so many difficulties after 30000(!) of tunings.
Also using the scale matrix from NUTS doesnt means that you can have good performance in HMC. The advantage of NUTS is the tree building that allows the exploration of the parameter space where fix integration time of HMC cannot archive good performance at.
All in all this is a bad idea, but if you insist, try the following:
chains = 4
cov = np.atleast_1d(pm.trace_cov(trace))
start = list(np.random.choice(trace, chains))
potential = quadpotential.QuadPotentialFull(cov)
step = pm.HamiltonianMC(potential=potential)
trace_hmc = pm.sample(1000, tune=1000, step=step, start=start)
Thank you @junpenglao . I will start a new topic to discuss my model.
@junpenglao As I was writing up the description of what I was doing in the model I actually found that I coded not what was intended! So thanks again, you were right all this time (you know what I mean) that the problem is in the model. Sorry for being so reluctant to heed of your advice.
LOL, no problem. Glad it helps!