Difference in posterior between Metropolis and NUTS Sampler

Questions v5
, ,
1

Hello,
I try to sample different polymer-chains with length N in spherical-coordinates, and than calculate the center of mass in Cartesian coordinates. Theta and phi (angles between monomers) have a uniform prior distribution.
When running this with NUTS, I get a Gausian distribution for center of mass, with mean at 0. When i run this with Metropolis i get a Gausian, but for the x axis the mean is shifted. Seen in the plots below. I use 10 chains and 8000 draws. For me it looks like it converged(I am quite new to this). How can this happen, do I have a mistake in my parametrization?
(Using pymc 5.0.2)
My Model:

N=40
processors=10
with pm.Model() as chain:
    
    theta = pm.Uniform("theta", lower=-np.pi/2, upper=np.pi/2, size=N)
    phi = pm.Uniform("phi", lower=0, upper=2*np.pi, size=N)
    
    ct = theta.cos()
    st = theta.sin()
    
    cp = phi.cos()
    sp = phi.sin()
    
    _x, _ = pt.scan(lambda _ct, _st, _cp, _sp, prior: prior + [_cp*_ct, _sp*_ct, _st],
               sequences = [ct, st, cp, sp],
               outputs_info = pt.tensor.zeros(3))

    x = _x[:,0]
    y = _x[:,1]
    z = _x[:,2]
      
    
    pm.Deterministic("z_mean", z.mean())
    pm.Deterministic("x_mean", x.mean())
    pm.Deterministic("y_mean", y.mean())

    #This lets the energy-less distribution be uniform on the sphere
    pm.Potential("sphere", ct.log().sum())

    
    #step=[pm.step_methods.Metropolis([theta, phi])]
    step=[pm.NUTS([theta, phi], target_accept=0.9)]

    #datat = pm.smc.sample_smc(draws=500,cores=processors)
    datat = pm.sample(draws = 2000, cores = processors, step= step)

NUts_forum_01
NUts_forum_011192Ă—592 245 KB

metropolis_without_potential
metropolis_without_potential1192Ă—592 208 KB

1 Like
2

Welcome!

It is not entirely clear what is going on here. If you have a runnable example (e.g., something simpler than what you are currently using) someone might be able to investigate further.

A couple of suggestions. First, I would install the most recent version of PyMC if you can mange to do so. Second, I would suggest avoiding explicit specification of the step methods and allow pm.sample() to automatically select step methods (including compound methods). There is not particular reason to use Metropolis except in relatively specific scenarios. If NUTS works for you model, you should probably use it. Third, 80,000 samples is likely excessive unless your model is sampling very very poorly (at which point you should just fix your model).

1 Like
3

Thanks for the quick response.
Updating to 5.6.1 solved the problem.
I guess the problem was fixed some time ago.

Regarding why i use Metropolis instead of NUTS;
I want to use SMC(for tempering), where i had the same problem(updating also solved this problem).

1 Like
4

Maybe i spoke to soon.

Updated to 5.11.0 it solved the problem for the Metropolis sampler.
Now i tried to use the SMC sample, resulting in the same problem. As seen in the pictures below, x_mean goes to 0 the more draws I use.

I don’t know if I should open a new thread for this?

Questions I have:
Why is this happening for x and not for y?
The different chains should be uncorrelated, so why are they converging to the same posterior which then converges to mean=0?(maybe I misunderstand how the SMC alg. works)

Also to note updating to 5.11.0 solved the problem for the “normal” Metropolis sampler but not for SMC.

This time I also added a exported jupyter notebook with my model. I couldn’t think of something more simpler. I excluded the z-axes.

smc_2000
smc_20001192Ă—392 147 KB

smc_5000draws
smc_5000draws1191Ă—392 113 KB

For_forum_2d.py (1.7 KB)