Rhat Nan Metropolis

I’m trying to see whether there was a switch point in Syrian conflict casualties (2011-2014, grouped by casualties per week).

with pm.Model() as model:
    alpha = 1.0/a.casualties.mean() 
    lambda_1 = pm.Exponential("lambda_1", alpha)    
    lambda_2 = pm.Exponential("lambda_2", alpha)    
    tau = pm.DiscreteUniform("tau", lower=0, upper=count_conflict_time-1)#switchpoint
   
    idx = np.arange(count_conflict_time) # Index
    lambda_ = pm.math.switch(tau > idx, lambda_1, lambda_2)
    observation = pm.Poisson("obs", lambda_, observed=a.casualties)

with model:
    trace = pm.sample(10000, tune=5000)

pymc3 automatically chooses the right sampler:

CompoundStep
NUTS: [lambda_2_log__, lambda_1_log__]
Metropolis: [tau]

But when look at summary (pm.summary(trace)) I see Rhat is NaN:

     mean    sd 	mc_error 	hpd_2.5 	hpd_97.5 	n_eff 	Rhat
 tau  78.0, 0.0      0.0,      78.0,         78.0,      20000.0, NaN

And by plotting trace we see:

No matter what I do, I always get scale reduction factors (rhat) = NaN. Does this mean the tau didn’t converge? How can this be corrected?

As you can see from the trace, tau is completely stuck. This is likely due to the inefficient of random walk Metropolis in high dimension: all proposals are rejected.

While Discrete switch point model could be sampled in PyMC3 after some careful tuning, our recommendation is to rewrite them as a continuous model. You can have a look of this notebook to see how to do it:

1 Like

thank you, junpenglao. When I try to run your code above, I get:

Can’t get attribute ‘Composed’ on <module ‘main’ (built-in)>

I’ve been trying for hours to fix that to no avail. Do you have any idea what is wrong?

This suggestion says to write the function into a different file, but it didn’t seem to work.

1 Like