How does PYMC sample a univariate distribution with no observed values?

Hi Pymc3 community,

I was wondering if anyone could point me to the documentation or paper that explains in the simplest terms how PYMC sample from a simple univariate distribution.
For example

import pymc3 as pm
import numpy as np

def log_likelihood(x, alpha):
    """
    Power law mass function
    """
    return np.log(x**-alpha)

with pm.Model() as model:
   x= pm.Uniform(0.01, 1)
   alpha= pm.Uniform(1, 5)
   like = pm.Potential('likelihood', log_likelihood(x, alpha))
   trace= pm.sample(draws=1000,  cores=4,   init='advi')

I think I roughly understand how this would work for a multidimensional belief network with observed values, but I’m having a hard time understanding how this works for this simple example.

Thank you,

1 Like

There’s a couple of steps.

  1. After the variable is added to model, PyMC3 calculates the log-posterior as a function of x. In this case, since the prior density is uniform and there is no likelihood, the log posterior is constant.

  2. Since the uniform random variable x is bounded, PyMC3 creates a derived variable which is unbounded by applying a transformation f such that f(x) \in (-\infty, \infty) even though x\in [0,1]. This makes it easier to do sampling later.

  3. Sampling is performed with Markov chain Monte Carlo. Specifically, a variant of Hamiltonian Monte Carlo is used which calculates gradients of the log posterior (again, equal to zero here since the posterior is constant).

Note that this approach won’t work for all priors - some of the priors are improper meaning they don’t integrate to 1. For these cases, we can’t sample. If you change Uniform to Flat in the above example, you will get an error.

3 Likes

Thank you! I didn’t mean to make the posterior constant, I forgot to add a potential term. (See edited version)

1 Like

A comment here would be that I wasn’t able to figure it out from the documentation or the paper https://arxiv.org/abs/1507.08050