3
$\begingroup$

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......."
shear=np.loadtxt('/vol/cluster.shear.cat') #my shear catalogue
data=np.loadtxt('/vol/Extended.MockFidSurvey.dat')
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....."
redshift_pdf=np.loadtxt('/vol/data/pdf.redshift.MockFidSurvey.dat')
#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')
    M.use_step_method(pm.AdaptiveMetropolis, M.model_pars ,verbose=1)
    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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.