How are you using the bernoulli_loglh in your model? Via pm.DensityDist? One way to get only integer values is to create a class Distribution that inherits from pm.Discrete, but the only thing this will do is tell the Metropolis sampler to round it’s proposed values before evaluating the logp. You can achieve the same if you round your value in your _logp and then save the rounded accepted values in a pm.Deterministic Something like
with pm.Model() as m:
my_dist = pm.DensitiyDist('my_dist', bernoulli_loglh...)
my_values = pm.Deterministic('my_values', np.round(my_dist))
Same thing if you are using a pm.Potential…