Difference between 'jitter+adapt_diag' and 'adapt_diag'?

Questions
1

hello everyone,
I am a beginner at Pymc3. I built a nonlinear regression model y=a(x-b)^c to fit the observed data (x,y). This model assumes that the observation of x is equal to the true value while the observation of y is not. I have constructed prior distribution of 5 parameters and the likelihood using my expertise. I found a strange phenomenon in my simulation and I don’t know why. Hoping to get an explanation. My code is as follows:

data = data.sort_values('Hgauged')
Hgauged = np.array(data['Hgauged'])
Qgauged = np.array(data['Qgauged'])
sigma_Qmeasurement = data['Qgauged']*(data['uQ[%]']/100)/1.96
with pm.Model() as model:
    hydraulic_param_a = pm.Normal('hydraulic_param_a', mu=82.16, sd=22.46)
    hydraulic_param_b = pm.Normal('hydraulic_param_b', mu=-0.5, sd=0.25)
    hydraulic_param_c = pm.Normal('hydraulic_param_c', mu=1.67, sd=0.025)
    error_param_gamma = pm.Uniform('error_param_gamma', 0, 10000, shape=2)
    Qratingcurve = pm.Deterministic('Qratingcurve', hydraulic_param_a*(Hgauged - hydraulic_param_b)**hydraulic_param_c)
    sigma_Fratingcurve = error_param_gamma[0] + error_param_gamma[1]*Qratingcurve
    sigma_total = np.sqrt(sigma_Fratingcurve**2 + sigma_Qmeasurement**2)
    Q_pred2 = pm.Normal('Q_pred2', mu=Qratingcurve, sd=sigma_total, observed=Qgauged)
    step = pm.NUTS()
    trace=pm.sample(20000, step=step, cores=1)

At first, I use the default init (init=‘jitter+adapt_diag’) like above，the iteration is very slow at first chain( more than 15 mins) in 1 jobs, while the second chain is obviously faster (about 3 mins). But the traceplot is like this:

traceplot2
traceplot2.png856×424 58.7 KB

Then, I choose the init=‘adapt_diag’, the calculating speed is up to about 3 mins for one chain. The traceplot is like this:
traceplot
traceplot.png7200×3600 1.01 MB

But when I shut down my Spyder(Python 3.7), and turned it on again to run my model with init=‘adapt_diag’, the calculation speed again became very slow, and the result was the same as the first picture.
I repeat the above process several times, and get the same traceplot ficture.
I want to know the difference between ‘jitter+adapt_diag’ and ‘adapt_diag’, and how to set the parameters of pm.sample() to get the satisfied results like the second figure and the iteration speed can also be fast. Thank you very much!

2

jitter+adapt_diag makes a jitter on the starting value, so that different chain started at a different point in the parameter space. It is the recommended option, as it helps you discover bad parameterization. For example, in your case there are indication that some parameter space is difficult to sample from - narrow geometry etc - which makes some chain very slow (as it takes smaller and larger quantities of leapfrog steps).
So unless you have some edge cases where jitter+adapt_diag makes the staring value invalid (nan or inf logp), you should always use jitter+adapt_diag.

In your code, the likely problem is error_param_gamma = pm.Uniform('error_param_gamma', 0, 10000, shape=2) - HalfNormal is usually a much better choice of prior for sigma.

1 Like
3

thanks junpeng very much! I solved the problem using a HalfNormal distribution to describe the prior of my error parameters. the iteration speed is also good, thanks very much!