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 + error_param_gamma*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:
Then, I choose the init=‘adapt_diag’, the calculating speed is up to about 3 mins for one chain. The traceplot is like this:
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!