# pyMC produces values outside range of uniform distribution while sampling from Bayesian hierarchical model [closed]

I have a hierarchical Bayesian model consisting of a Uniform prior distribution, between a minimum and maximum value (hyperparameters) at the top level of the hierarchy. I sample a "mean" from the top level and generate a new uniform distribution for the next level, centred around that mean. I also constrain lower level distributions to lie between the minimum and maximum specified at the top level using python's min and max functions. I would like to do this by truncating the lower level Uniform distributions, but truncating Uniforms is not allowed in pyMC.

While looking at samples from the posterior, I noticed that occasionally, samples from nodes below the top level of the hierarchy were coming from outside the range specified by the hyperparameters.

The minimial code example below demonstrates the problem (at least running pyMC 2.3.6, on my machine). What have I done wrong here? How can I ensure that K_a and K_b always lie between K_min and K_max?

#Simplified hierarchical model

import pymc as mc

# Hyperparameters
# No sample at any level of hierarchy should be outside this range
K_min,K_max=0.0,0.5

priors={}

K=mc.Uniform('K',K_min,K_max)
K_delta=mc.Uniform('K_delta',0,(K_max-K_min)/2)

priors["K"]=K
priors["K_delta"]=K_delta

grplabs=["a","b"]
for grplab in grplabs:
priors["K_{}".format(grplab)]=mc.Uniform("K_{}".format(grplab),max(K_min,K-K_delta),min(K_max,K+K_delta))

M=mc.MCMC(priors)

K=[]
K_a=[]
K_b=[]

for i in xrange(1000):
M.draw_from_prior()
K.append(float(M.K.get_value()))
K_a.append(float(M.K_a.get_value()))
K_b.append(float(M.K_b.get_value()))

assert min(K)>=0 # This passes
assert min(K_a)>=0 # This fails
assert min(K_b)>=0 # This fails


Ok, it seems that max and min are not among the "elementary operations" which have syntactic sugar support when building deterministic variables in pyMC. A method to generate a prior with the expected samples by using python decorators to define a new deterministic operation is written below.

It might be safer to have some kind of warning when a user attempts to carry out an unsupported operation on pyMC objects, rather than just ignoring the operation?

#Simplified hierarchical model

import pymc as mc
import numpy as np

# Hyperparameters
# No sample at any level of hierarchy should be outside this range
K_min,K_max=0.0,0.5

priors={}

K=mc.Uniform('K',K_min,K_max)
K_delta=mc.Uniform('K_delta',0,(K_max-K_min)/2)

priors["K"]=K
priors["K_delta"]=K_delta

@mc.deterministic(plot=False)
def cliplowK(K=K,K_min=K_min,K_delta=K_delta):
return(np.maximum(K_min,K-K_delta))

@mc.deterministic(plot=False)
def cliphighK(K=K,K_min=K_min,K_delta=K_delta):
return(np.minimum(K_max,K+K_delta))

grplabs=["a","b"]
for grplab in grplabs:
priors["K_{}".format(grplab)]=mc.Uniform("K_{}".format(grplab),cliplowK,cliphighK)

M=mc.MCMC(priors)

K_vals=[]
K_a_vals=[]
K_b_vals=[]

for i in xrange(100000):
M.draw_from_prior()
K_vals.append(float(M.K.get_value()))
K_a_vals.append(float(M.K_a.get_value()))
K_b_vals.append(float(M.K_b.get_value()))

assert min(K_vals)>=0 # This passes
assert min(K_a_vals)>=0 # This passes
assert min(K_b_vals)>=0 # This passes