Solved – How do we define log-normal prior and a multivariate posterior log-likelihood in PyMC

bayesiancomputational-statisticshierarchical-bayesianmarkov-chain-montecarlopymc

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 are observational data and correspondence error respectively and 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

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)

Best Answer

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)

Related Solutions

How to Compute a Posterior Density Estimate from a Prior and Likelihood?

You have several things mixed up. The theory talks about multiplying the prior distribution and the likelihood, not samples from the prior distribution. Also it is not clear what you have the prior of, is this a prior on the mean of something? or something else?

Then you have things reversed in the likelihood, your observations should be x with either prior draws or known fixed constants as the mean and standard deviation. And even then it would really be the product of 4 calls to dnorm with each of your observations as x and the same mean and standard deviation.

What is really no clear is what you are trying to do. What is your question? which parameters are you interested in? what prior(s) do you have on those parameters? are there other parameters? do you have priors or fixed values for those?

Trying to go about things the way you currently are will only confuse you more until you work out exactly what your question is and work from there.

Below is added due after the editing of the original question.

You are still missing some pieces, and probably not understanding everything, but we can start from where you are at.

I think you are confusing a few concepts. There is the likelihood that shows the relationship between the data and the parameters, you are using the normal which has 2 parameters, the mean and the standard deviation (or variance, or precision). Then there is the prior distributions on the parameters, you have specified a normal prior with mean 0 and sd 1, but that mean and standard deviation are completely different from the mean and standard deviation of the likelihood. To be complete you need to either know the likelihood SD or place a prior on the likelihood SD, for simplicity (but less real) I will assume we know the likelihood SD is $\frac12$ (no good reason other than it works and is different from 1).

So we can start similar to what you did and generate from the prior:

> obs <- c(0.4, 0.5, 0.8, 0.1)
> pri <- rnorm(10000, 0, 1)

Now we need to compute the likelihoods, this is based on the prior draws of the mean, the likelihood with the data, and the known value of the SD. The dnorm function will give us the likelihood of a single point, but we need to multiply together the values for each of the observations, here is a function to do that:

> likfun <- function(theta) {
+ sapply( theta, function(t) prod( dnorm(obs, t, 0.5) ) )
+ }

Now we can compute the likelihood for each draw from the prior for the mean

> tmp <- likfun(pri)

Now to get the posterior we need to do a new type of draw, one approach that is similar to rejection sampling is to sample from the prior mean draws proportional to the likelihood for each prior draw (this is the closest to the multiplication step you were asking about):

> post <- sample( pri, 100000, replace=TRUE, prob=tmp )

Now we can look at the results of the posterior draws:

> mean(post)
[1] 0.4205842
> sd(post)
[1] 0.2421079
> 
> hist(post)
> abline(v=mean(post), col='green')

and compare the above results to the closed form values from the theory

> (1/1^2*mean(pri) + length(obs)/0.5^2 * mean(obs))/( 1/1^2 + length(obs)/0.5^2 )
[1] 0.4233263
> sqrt(1/(1+4*4))
[1] 0.2425356

Not a bad approximation, but it would probably work better to use a built-in McMC tool to draw from the posterior. Most of these tools sample one point at a time not in batches like above.

More realistically we would not know the SD of the likelihood and would need a prior for that as well (often the prior on the variance is a $\chi^2$ or gamma), but then it is more complicated to compute (McMC comes in handy) and there is no closed form to compare with.

The general solution is to use existing tools for doing the McMC calculations such as WinBugs or OpenBugs (BRugs in R gives an interface between R and Bugs) or packages such LearnBayes in R.

Solved – Posterior autocorrelation in Pymc. How to interpret it

Autocorrelation dictates the amount of time you have to wait for convergence. If autocorrelation is high, you will have to use a longer burn-in, and you will have to draw more samples after the burn-in to get a good estimate of the posterior distribution.
Low autocorrelation means good exploration. Exploration and convergence are essentially the same thing. If a sampler explores well, then it will converge fast, and vice versa. Part of the confusion here is that the book does not explain the concept of convergence very well. You can find a much better definition of convergence on Wikipedia.

Best Answer

Related Solutions

How to Compute a Posterior Density Estimate from a Prior and Likelihood?

Solved – Posterior autocorrelation in Pymc. How to interpret it

Related Question