0
$\begingroup$

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
$\endgroup$

closed as off-topic by Xi'an, Peter Flom Apr 25 at 12:21

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – Xi'an, Peter Flom
If this question can be reworded to fit the rules in the help center, please edit the question.

0
$\begingroup$

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
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.