Defining constraint on prior with potential class

I have written an MCMC code in order to estimate parameters Xpos, Ypos, MASS and concentration with a set of input data gal_pos, observed_g (the parameter which the model gets compared) and the uncertainty in this parameter, shearerr. The problem is that the MCMC chain will break down if MASS would be less than $10^{13}$ due to the properties of the likelihood function nfwmodeltools.shearprofile_like and it follows the exponential and LogNormal prior shape. In order to put the constraint on MASS and concentration, I have defined a potential class but I don't know how it should be used and where, in order to enforce MASS to follow lower bound.

import pymc as pm
import nfwmodeltools
import numpy as np
import math
import random
import matplotlib.pyplot as plt
from scipy.stats import expon
import astropy.cosmology
from astropy.cosmology import FlatLambdaCDM
cos = FlatLambdaCDM(H0=70, Om0=0.3)
from astropy import units as u
cosmo=nfwmodeltools.global_cosmology()
print "Reading input catalogues as observed values......."
for i in range(len(data[:,0])):
if (i<1):
gal_pos=np.array([data[i,0],data[i,1]])
observed_g=np.array([shear[i,2],shear[i,3]])
shearerr=np.array([data[i,5]])
else:
gal_pos=np.vstack((gal_pos, np.array([data[i,0],data[i,1]])))
observed_g=np.vstack((observed_g, np.array([shear[i,2],shear[i,3]])))
shearerr=np.vstack((shearerr,np.array([data[i,5]])))

g_err=shearerr.reshape(shearerr.shape[0],)
print "Reading redshift distributions of galaxies, where we will marginalize over this parameter when we will compute beta and rho_c_over_sigma_c....."
#Normalize the probability
zpdf=np.array([redshift_pdf[j,:]/sum(redshift_pdf[j,:]) for j in xrange(gal_pos.shape[0])] )
#The whole range of redshift
z=np.arange(0,1.5,0.001)
z_halo=0.15
x=1*u.Mpc
angular_separation=(x*180)/np.pi/ cos.angular_diameter_distance(z_halo)
print "convert radian to degree for angular separation in the redshift of cluster and compute 1 Mpc seperation in degree:"
print angular_separation
#one of the input parameter for loglikelihood
rho_crit=cosmo.cosmology.rho_crit(z_halo)
#prepare the requirement for the mcmc model
print "We multiply the probability with the computed function for the given array z to marginalize over z"
rho_c_over_sigma_c=np.dot(zpdf,cosmo.cosmology.RhoCrit_over_SigmaC(z, z_halo))
print rho_c_over_sigma_c.shape
beta=np.dot(zpdf, cosmo.cosmology.beta_s(z, z_halo))
print "Defining the model as nfw model for MCMC..."
def nfw(gal_pos,observed_g,g_err,rho_crit,rho_c_over_sigma_c,beta):
"""
gal_pos:background galaxy positions
gt: tangential reduced shear
g_err: shear error
rho_crit: critical density
rho_c_over_sigma_c: critical density over critical surface mass density
beta:Dds/Ds
"""
@pm.stochastic(dtype=np.float, observed=False, trace=True)
def Xpos(value=25.4,x_l=25.25,x_h=25.55):
"""The probable region of the position of halo centre"""
if ((value>x_h) or (value<x_l)):
return -np.inf
else:
return -np.log(x_h-x_l+1)

@pm.stochastic(dtype=np.float, observed=False, trace=True)
def Ypos(value=-10.01,y_l=-10.3,y_h=-9.9):
"""The probable region of the position of halo centre"""
if ((value>y_h) or (value<y_l)):
return -np.inf
else:
return -np.log(y_h-y_l+1)
#------------------------------------------------------------
# Based on a simulated or observed mass function we adopt a simple exponential prior
@pm.stochastic
def MASS(value=math.pow(10,14), rate = math.pow(10,15)):
"""mass is a stochastic parameter with exponential distribution.p(M)~exp(-M/10^15)"""
return pm.exponential_like(value, rate)
@pm.potential
def MASS_bound(MASS=MASS):
if ((MASS >= math.pow(10,13)) and (MASS < math.pow(10,16))):
return 0.0
else:
return -np.inf

#MASS=math.pow(10,30)*pm.Exponential('mass', beta=math.pow(10,15))
@pm.deterministic
def sigma( name='sigma_concentration' ,M=MASS, trace=True, plot= True):
if M < 10**15:
return .09
else:
return .06

cExpected = 5.26/(1+z_halo)*(MASS/math.pow(10,14))**(-.1) # based on Neto et al. 2007
concentration = pm.Lognormal("concentration", cExpected, sigma)
@pm.potential
def concentration_bound(concentration=concentration):
if (concentration <= 20):
return 0.0
else:
return -np.inf

model_pars=[Xpos, Ypos, MASS, concentration]
@pm.stochastic( name='reduced_shear', dtype=float,observed=True, trace = True )
def reduced_shear(value=observed_g, model_pars=model_pars):
Xpos=model_pars[0]
Ypos=model_pars[1]
mass=model_pars[2]
conc=model_pars[3]
#seperation in the scale of Mpc
dist=np.sqrt((gal_pos[:,0]-Xpos)**2+(gal_pos[:,1]-Ypos)**2)*angular_separation
phi=np.arctan2((gal_pos[:,1]-Ypos), (gal_pos[:,0]-Xpos))
gtan=-(value[:,0]*np.cos(2*phi)+value[:,1]*np.sin(2*phi))
gcros=-value[:,0]*np.sin(2*phi)+value[:,1]*np.cos(2*phi)
d=np.vstack((dist,gtan,g_err,beta,rho_c_over_sigma_c))
#concatenate input data and sort data based on their distance from cluster center
ndata=np.array(sorted(d.T, key=lambda  l:l[0]))
#choose objects closer than 2.5 Mpc to the cluster center
ndata=ndata[ndata[:,0]<=2.5*angular_separation]
#binning data again based on their distance
bins = np.linspace(ndata[0,0], ndata[-1,0], 10)
#compute the mean in each bin for different input parameters
digitized = np.digitize(ndata[:,0], bins)
bin_r_mpc= np.array([ndata[digitized == i,0].mean() for i in range(1, len(bins))])
bin_shear= np.array([ndata[digitized == i,1].mean() for i in range(1, len(bins))])
bin_shearerr= np.array([ndata[digitized == i,2].mean() for i in range(1, len(bins))])
avebeta= ndata[:,3].mean()
avebeta2=(ndata[:,3]**2).mean()
ave_rho_c_over_sigma_c=ndata[:,4].mean()
loglikelihood = nfwmodeltools.shearprofile_like(mass,
conc,
bin_r_mpc,
bin_shear,
bin_shearerr,
avebeta,
avebeta2,
rho_crit,
ave_rho_c_over_sigma_c)
return loglikelihood
return locals()

if __name__ == '__main__':

M = pm.MCMC(nfw(gal_pos,observed_g,g_err,rho_crit,rho_c_over_sigma_c,beta),db='pickle',dbname='NFWTracer.pickle')
mass = M.step_method_dict[M.MASS][0]
M.isample(100000,10000,100)
[pm.Matplot.plot(s) for s in [ M.MASS, M.Xpos, M.Ypos, M.concentration]]


I will appreciate if somebody give me advise that where it should be used and which kind of change I need to make in my code..

running the code with above configuration I got this error:

ERROR: ZeroProbability: Potential concentration_bound forbids its parents' current values [pymc.PyMCObjects]
Traceback (most recent call last):
File "MASS_NFW_MCMC.py", line 148, in <module>
M = pm.MCMC(nfw(gal_pos,observed_g,g_err,rho_crit,rho_c_over_sigma_c,beta),db='pickle',dbname='NFWTracer.pickle')
File "MASS_NFW_MCMC.py", line 102, in nfw
def concentration_bound(concentration=concentration):
File "/vol/anaconda/lib/python2.7/site-packages/pymc/InstantiationDecorators.py", line 215, in potential
return instantiate_pot(__func__)
File "/vol/anaconda/lib/python2.7/site-packages/pymc/InstantiationDecorators.py", line 208, in instantiate_pot
return Potential(parents=parents, **kwds)
File "/vol/anaconda/lib/python2.7/site-packages/pymc/PyMCObjects.py", line 275, in __init__
if not isinstance(self.logp, float):
File "/vol/anaconda/lib/python2.7/site-packages/pymc/PyMCObjects.py", line 327, in get_logp
raise ZeroProbability(self.errmsg)
pymc.Node.ZeroProbability: Potential concentration_bound forbids its parents' current values