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.
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.
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.
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.