NUTS sampler converges to wrong value

Good question! I am actually suspicious of it being multimodal…

Here is the model:

raw_py is my data:

prior_l = np.zeros(len(raw_py))
prior_u = raw_py/(eff*pFE)

trace_out = []
resmat = resp
truthdim = len(resmat)
model = mc.Model()    
    
with model:

truthprior = []

priormethod = getattr(mc,'Uniform')

for bin,(l,u) in enumerate(zip(prior_l,prior_u)):
    name = 'truth%d'%bin
    default_args = dict(name=name,lower=l,upper=u)
    args = dict(list(default_args.items()))

    prior = priormethod(**args)

    truthprior.append(prior)

truth = mc.math.stack(truthprior)
            
tresmat = np.array(resmat)
    
reco = theano.dot(truth, tresmat)
out = reco
  


unfolded = mc.Poisson('unfolded', mu=out,
                              observed=np.array(raw_py))  



step = mc.NUTS() 
trace = mc.sample(2000, step=step, init='advi_map', tune=1000, cores=2, 
                          chains=3, nuts_kwargs=None,
                          discard_tuned_samples=True,
                          progressbar=True)

If you need more information just ask!