setting up MCMC with log-likelihood and log-normal prior with PyMC

Question

I am a newbie with pyMC and I am not still able to construct the structure of my MCMC with pyMC. I would like to establish a chain and I am confused how to define my parameters and log-likelihood function together. My chi-squared function is given by:

where and are observational data and correspondence error respectively and is the model with four free parameter and the parameters are non-linear.

The prior for X and Y are uniform like:

import pymc as pm
import numpy as np
import math
import random

@pm.stochastic(dtype=np.float, observed=False, trace=True)
def Xpos(value=1900,x_l=1851,x_h=1962):
    """The probable region of the position of halo centre"""
    def logp(value,x_l,x_h):
        if ((value>x_h) or (value<x_l)):
       return -np.inf
    else:
       return -np.log(x_h-x_l+1)
    def random(x_l,x_h):
        return np.round((x_h-x_l)*random.random())+x_l

@pm.stochastic(dtype=np.float, observed=False, trace=True)
def Ypos(value=1900,y_l=1851,y_h=1962):
    """The probable region of the position of halo centre"""
    def logp(value,y_l,y_h):
        if ((value>y_h) or (value<y_l)):
       return -np.inf
    else:
       return -np.log(y_h-y_l+1)
    def random(y_l,y_h):
        return np.round((y_h-y_l)*random.random())+y_l

but for M and C are given as following:

where the mean of C is computed via

For M and C, the priors should look like this:

    M=math.pow(10,15)*pm.Exponential('mass', beta=math.pow(10,15))

    @pm.stochastic(dtype=np.float, observed=False, trace=True)
    def concentration(value=4, zh, M200):
        """logp for concentration parameter"""
        def logp(value=4.,zh, M200):
            if (value>0):
           x = np.linspace(math.pow(10,13),math.pow(10,16),200 )
           prob=expon.pdf(x,loc=0,scale=math.pow(10,15))
           conc = [5.26/(1.+zh)*math.pow(x[i]/math.pow(10,14),-0.1) for i in range(len(x))]
           mu_c=0
           for i in range(len(x)):
               mu_c+=prob[i]*conc[i]/sum(prob)
           if (M200 < pow(10,15)):
              tau=1./(0.09*0.09)
           else:
              tau=1./(0.06*0.06)
               return  pm.lognormal_like(value, mu_c, tau)
            else
               return -np.inf
        def random(mu_c,tau):
            return np.random.lognormal(mu_c, tau, 1)

The parameter z is also a constant in C prior. I am wondering how I could define my likelihood for , and should it be referred as @Deterministic variable? Have I defined M and C as priori information in a correct way or not?

I will be grateful if somebody gives me some tips that how I can combine these parameters with given priors.

Answer 1

#priors
@pm.stochastic(dtype=np.float, observed=False, trace=True)
def Xpos(value=1900,x_l=1800,x_h=1950):
    """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=1750,y_l=1200,y_h=2000):
    """The probable region of the position of halo centre"""
    def logp(value,y_l,y_h):
        if ((value>y_h) or (value<y_l)):
       return -np.inf
    else:
       return -np.log(y_h-y_l+1)


M=math.pow(10,15)*pm.Exponential('mass', beta=math.pow(10,15))


@deterministic
def sigma(value = 1, M=M): 
   if M < 10**15:
       return .09
   else:
       return .06

cExpected = 5.26/(1+z)*(M/math.pow(10,14))**(-.1) # based on Neto et al. 2007
concentration = Lognormal("concentration", cExpected, sigma)


#model
@pm.deterministic( name='reduced_shear', dtype=np.float, observed=False, trace = True )
def reduced_shear(x=Xpos,y=Ypos,mass=M,conc=concentration):
    nfw = NFWHalo(mass,conc,zh=0.128,[x,y])
    g1tot=0;g2tot=0
    for i in range(len(z)):
        g1,g2,magnification=nfw.getLensing( gal_pos, z[i])
        g1tot+=g1*redshift_pdf[i]/sum(redshift_pdf)
        g2tot+=g2*redshift_pdf[i]/sum(redshift_pdf)
    theta=arctan2(gal_ypos - Ypos, gal_xpos - Xpos)
    value=-g1tot*cos(2*theta)-g2tot*sin(2*theta) #tangential shear
    return value

@pm.deterministic( name='reduced_shear', dtype=np.float, observed=False, trace = True )
def tau_shear(Xpos,Ypos,M,concentration):
    nfw = NFWHalo(M,concentration,zh=0.128,[Xpos,Ypos])
    g1tot=0;g2tot=0
    for i in range(len(z)):
        g1,g2,magnification=nfw.getLensing( gal_pos, z[i])
        g1tot+=g1*redshift_pdf[i]/sum(redshift_pdf)
        g2tot+=g2*redshift_pdf[i]/sum(redshift_pdf)
        theta=arctan2(gal_ypos - Ypos, gal_xpos - Xpos)
    gt=-g1tot*cos(2*theta)-g2tot*sin(2*theta)
    g_squared=g1tot**2+g2tot**2
    delta_abse=sqrt(delta_e1**2+delta_e1**22)
    value=(1-g_squared)*delta_abse
    return value


tau = pm.Normal('tau', tau_shear, 0.2)
#likelihood
obs = pm.Normal("obs", mu=reduced_shear, tau, value=data, observed=True)