Unobserved Random Variables#

Every unobserved RV has the following calling signature: name (str), parameter keyword arguments. Thus, a normal prior can be defined in a model context like this:

withpm.Model():x=pm.Normal("x",mu=0,sigma=1)

As with the model, we can evaluate its logp:

pm.logp(x,0).eval()

array(-0.91893853)

Observed Random Variables#

Observed RVs are defined just like unobserved RVs but require data to be passed into theobserved keyword argument:

withpm.Model():obs=pm.Normal("x",mu=0,sigma=1,observed=rng.standard_normal(100))

observed supports lists,numpy.ndarray andpytensor data structures.

Deterministic transforms#

PyMC allows you to freely do algebra with RVs in all kinds of ways:

withpm.Model():x=pm.Normal("x",mu=0,sigma=1)y=pm.Gamma("y",alpha=1,beta=1)plus_2=x+2summed=x+ysquared=x**2sined=pm.math.sin(x)

Though these transformations work seamlessly, their results are not stored automatically. Thus, if you want to keep track of a transformed variable, you have to usepm.Deterministic:

withpm.Model():x=pm.Normal("x",mu=0,sigma=1)plus_2=pm.Deterministic("x plus 2",x+2)

Note thatplus_2 can be used in the identical way to above, we only tell PyMC to keep track of this RV for us.

Lists of RVs / higher-dimensional RVs#

Above we have seen how to create scalar RVs. In many models, we want multiple RVs. Users will sometimes try to create lists of RVs, like this:

withpm.Model():# bad:x=[pm.Normal(f"x_{i}",mu=0,sigma=1)foriinrange(10)]

This works, but it is slow and not recommended. Instead, we can usecoordinates:

coords={"cities":["Santiago","Mumbai","Tokyo"]}withpm.Model(coords=coords)asmodel:# good:x=pm.Normal("x",mu=0,sigma=1,dims="cities")

x is now a array of length 3 and each of the 3 variables within this array is associated with a label. This will make it very easy to distinguish the 3 different variables when we go to look at results. We can index into this array or do linear algebra operations on it:

withmodel:y=x[0]*x[1]# indexing is supportedx.dot(x.T)# linear algebra is supported

Initialize Random Variables#

Though PyMC automatically initializes models, it is sometimes helpful to define initial values for RVs. This can be done via theinitval kwarg:

withpm.Model(coords={"idx":np.arange(5)})asmodel:x=pm.Normal("x",mu=0,sigma=1,dims="idx")model.initial_point()

{'x': array([0., 0., 0., 0., 0.])}

withpm.Model(coords={"idx":np.arange(5)})asmodel:x=pm.Normal("x",mu=0,sigma=1,dims="idx",initval=rng.standard_normal(5))model.initial_point()

{'x': array([-0.36012097, -0.16168135,  1.07485641, -0.08855632, -0.03857412])}

This technique is sometimes useful when trying to identify problems with model specification or initialization.

3.1 Sampling#

The main entry point to MCMC sampling algorithms is via thepm.sample() function. By default, this function tries to auto-assign the right sampler(s).pm.sample() returns anarviz.InferenceData object.InferenceData objects can easily be saved/loaded from a file and can carry additional (meta)data such as date/version and posterior predictive samples. Take a look at theArviZ Quickstart to learn more.

withpm.Model()asmodel:mu=pm.Normal("mu",mu=0,sigma=1)obs=pm.Normal("obs",mu=mu,sigma=1,observed=rng.standard_normal(100))idata=pm.sample(2000)

Auto-assigning NUTS sampler...Initializing NUTS using jitter+adapt_diag...Multiprocess sampling (2 chains in 2 jobs)NUTS: [mu]

Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 4 seconds.

As you can see, with model that exclusively contains continuous variables, PyMC assigns the NUTS sampler, which is very efficient even for complex models. PyMC also runs initial tuning to find good starting parameters for the sampler. Here we draw 2000 samples from the posterior in each chain and allow the sampler to adjust its parameters in an additional 1500 iterations.

If not set via thechains kwarg, the number of chains is determined from the number of available CPU cores.

idata.posterior.dims

Frozen({'chain': 2, 'draw': 2000})

The tuning samples are discarded by default. Withdiscard_tuned_samples=False they can be kept and end up in a separate group within theInferenceData object (i.e.,idata.warmup_posterior).

You can control how the chains are run in parallel using thechains andcores kwargs:

withpm.Model()asmodel:mu=pm.Normal("mu",mu=0,sigma=1)obs=pm.Normal("obs",mu=mu,sigma=1,observed=rng.standard_normal(100))idata=pm.sample(cores=4,chains=6)

Auto-assigning NUTS sampler...Initializing NUTS using jitter+adapt_diag...Multiprocess sampling (6 chains in 4 jobs)NUTS: [mu]

Sampling 6 chains for 1_000 tune and 1_000 draw iterations (6_000 + 6_000 draws total) took 7 seconds.

idata.posterior["mu"].shape

(6, 1000)

# get values of a single chainidata.posterior["mu"].sel(chain=2).shape

(1000,)

3.2 Analyze sampling results#

The most common used plot to analyze sampling results is the so-called trace-plot:

withpm.Model()asmodel:mu=pm.Normal("mu",mu=0,sigma=1)sd=pm.HalfNormal("sd",sigma=1)obs=pm.Normal("obs",mu=mu,sigma=sd,observed=rng.standard_normal(100))idata=pm.sample()

Auto-assigning NUTS sampler...Initializing NUTS using jitter+adapt_diag...Multiprocess sampling (2 chains in 2 jobs)NUTS: [mu, sd]

Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 4 seconds.

az.plot_trace(idata);

Another common metric to look at is the Gelman-Rubin statistic, or R-hat:

az.summary(idata)

R-hat is also presented as part of theaz.plot_forest:

az.plot_forest(idata,r_hat=True);

Finally, for a plot of the posterior that is inspired by[Kruschke, 2014], you can use the:

az.plot_posterior(idata);

For high-dimensional models it becomes cumbersome to look at the traces for all parameters. When usingNUTS we can look at the energy plot to assess problems of convergence:

withpm.Model(coords={"idx":np.arange(100)})asmodel:x=pm.Normal("x",mu=0,sigma=1,dims="idx")idata=pm.sample()az.plot_energy(idata);

Auto-assigning NUTS sampler...Initializing NUTS using jitter+adapt_diag...Multiprocess sampling (2 chains in 2 jobs)NUTS: [x]

100.00% [4000/4000 00:04<00:00 Sampling 2 chains, 0 divergences]

Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 5 seconds.

For more information on sampler stats and the energy plot, seeSampler Statistics. For more information on identifying sampling problems and what to do about them, seeDiagnosing Biased Inference with Divergences.

3.3 Variational inference#

PyMC supports various Variational Inference techniques. While these methods are much faster, they are often also less accurate and can lead to biased inference. The main entry point ispymc.fit().

withpm.Model()asmodel:mu=pm.Normal("mu",mu=0,sigma=1)sd=pm.HalfNormal("sd",sigma=1)obs=pm.Normal("obs",mu=mu,sigma=sd,observed=rng.standard_normal(100))approx=pm.fit()

Finished [100%]: Average Loss = 142

The returnedApproximation object has various capabilities, like drawing samples from the approximated posterior, which we can analyse like a regular sampling run:

idata=approx.sample(1000)az.summary(idata)

arviz - WARNING - Shape validation failed: input_shape: (1, 1000), minimum_shape: (chains=2, draws=4)

Thevariational submodule offers a lot of flexibility in which VI to use and follows an object oriented design. For example, full-rank ADVI estimates a full covariance matrix:

mu=pm.floatX([0.0,0.0])cov=pm.floatX([[1,0.5],[0.5,1.0]])withpm.Model(coords={"idx":np.arange(2)})asmodel:pm.MvNormal("x",mu=mu,cov=cov,dims="idx")approx=pm.fit(method="fullrank_advi")

Finished [100%]: Average Loss = 0.012772

An equivalent expression using the object-oriented interface is:

withpm.Model(coords={"idx":np.arange(2)})asmodel:pm.MvNormal("x",mu=mu,cov=cov,dims="idx")approx=pm.FullRankADVI().fit()

Finished [100%]: Average Loss = 0.020531

withpm.Model(coords={"idx":np.arange(2)})asmodel:pm.MvNormal("x",mu=mu,cov=cov,dims="idx")approx=pm.FullRankADVI().fit()

Finished [100%]: Average Loss = 0.014143

plt.figure()idata=approx.sample(10000)az.plot_pair(idata,var_names="x",coords={"idx":[0,1]});

<Figure size 432x288 with 0 Axes>

Stein Variational Gradient Descent (SVGD) uses particles to estimate the posterior:

w=pm.floatX([0.2,0.8])mu=pm.floatX([-0.3,0.5])sd=pm.floatX([0.1,0.1])withpm.Model()asmodel:pm.NormalMixture("x",w=w,mu=mu,sigma=sd)approx=pm.fit(method=pm.SVGD(n_particles=200,jitter=1.0))

withpm.Model()asmodel:pm.NormalMixture("x",w=[0.2,0.8],mu=[-0.3,0.5],sigma=[0.1,0.1])

plt.figure()idata=approx.sample(10000)az.plot_dist(idata.posterior["x"]);

For more information on variational inference, seeVariational Inference.