1
$\begingroup$

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:

enter image description here

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

The priors for X and Y are uniform but for M and C are given as following:

enter image description here ;

where the probability of c follows log-normal distribution while the expectation value of c is computed with the above formula and is the function of M and $\sigma$ is 0.09 if $M < 10^{15}$ otherwise $\sigma=0.06$:

enter image description here

enter image description here

for each C the parameter z is constant. I am wondering how I could define my likelihood for enter image description here , and should it be referred as @Deterministic variable? Did I define 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.

import pymc as pm
import numpy as np
import math
import random
from scipy.stats import expon

@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"""
    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=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)
    def random(y_l,y_h):
        return np.round((y_h-y_l)*random.random())+y_l


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): #c parameter
    """logp for concentration parameter"""
    def logp(value=4.,zh, M):
        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 (M < 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)
$\endgroup$
  • $\begingroup$ Is this a self-study question? $\endgroup$ – user44764 May 21 '14 at 23:48
  • $\begingroup$ Well @Matthew I am learning by myself how to use pyMC to code my problem. I am trying to define my variables correctly before starting any MCMC chain. $\endgroup$ – Dalek May 22 '14 at 7:41
  • $\begingroup$ $X, Y, M, C$ all have probability distributions, so they would be @Stochastic class. All you said about $\hat g$ is that it is a 'model', can you write out what that means exactly? Does $\hat g$ have random components? It seems that you have defined the prior of $C, M$ already? $\endgroup$ – user44764 May 22 '14 at 16:05
  • $\begingroup$ @Matthew $\hat{g}$ is a non-linear function of $(X,Y,M,C)$. I assume I can define a uniform prior for $X,Y$, since I have prior knowledge where they are roughly located and it might prevent of degeneracy between the results if I run MCMC with different initial conditions. I also think Jeffreys prior can describe $M$ very well but I do not know how to define $M,C$ priors and I appreciate if you can help. $\endgroup$ – Dalek May 22 '14 at 16:17
  • $\begingroup$ @Matthew I modify my question by adding what I have written so far and add some more details. I would like that somebody with the experience coding with pymc take a look and confirm I have coded my problem in the right way. $\endgroup$ – Dalek May 28 '14 at 11:35
0
$\begingroup$

Your @stochastic uses are not correct. Notice that your functions don't return anything. When using the decorator, you're supposed to return value of the logp. See here for example usage.

If you're going to use @stochastic I think you probably want something like this for each of your @stochastic uses.

@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)

(directly returning the logp)

If you need to provide a random() function (I suspect you don't), I think you can pass it to stochastic.

However, since you just want a uniform prior for Xpos and Ypos, you can just use Uniform instead.

Xpos = Uniform("Xpos", 1800, 1950)

Your concentration stochastic seems very complicated given that C is pretty straightforward. I would expect it to directly follow your definition of C.

g_hat should definitely be a @deterministic, since you say its a function. If so, it shouldn't have a likelihood of its own.

For concentration, you want something like

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

cExpected = const/(1+z)*M**-.1
concentration = Lognormal("concentration", cExpected, sigma)
$\endgroup$
  • $\begingroup$ I add more description to the problem. For x and y, I was following the pyMC manual and I thought if I want to instantiate Stochastic directly I need to also generate the random data for MCMC sampling later. On the other hand my main concern is how to define prior for c. $\endgroup$ – Dalek May 28 '14 at 20:14
  • $\begingroup$ Where were you looking in the manual? Perhaps that section can be improved. I've updated the answer. $\endgroup$ – John Salvatier May 28 '14 at 22:52
  • $\begingroup$ The example for Decorator Stochastic class section 4.2 for parameter s in 'PyMC:Bayesian Stochastic Modelling in Python'. I was wondering why the function random is defined, e.g. because in sampling parameter space MCMC needs to generate random steps and then test for the acceptance rate or not? $\endgroup$ – Dalek May 29 '14 at 7:27
  • $\begingroup$ if we go back to my problem, M parameter has an exponential distribution, then when I define cExpected which is the Expected value of c with a given functionality of M and its distribution, shouldn't I compute the expectation value with considering the exponential distribution of M and insert it in the concentration formula? What is the role of random() function in decorator shape of Stochastic class? $\endgroup$ – Dalek May 29 '14 at 15:21
  • $\begingroup$ No. M does not neccessarily have an exponential distribution, only an exponential prior distribution. I think I see what you're asking, whether you have to correct for the first parameter being the location parameter and not the expectation, and yes you do (pymc-devs.github.io/pymc/…). I don't think you need random, don't worry about it. Or look it up. $\endgroup$ – John Salvatier May 29 '14 at 16:01

Your Answer

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

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