Frequently Asked Questions

[junpenglao]

• Is find_MAP still encouraged in practice? i notice a lot of code samples stopped using find_MAP (asked by @kpei)

find_MAP is pretty much discouraged now in most situations. The MAP is not a representative point in high dimensions, because the MAP is usually not in the typical set (Betancourt’s paper is the best reference https://arxiv.org/abs/1701.02434). Moreover, it can make life difficult for NUTS, as the gradients are small at the MAP.

• Is NUTS sampler with ADVI always encouraged? I know that it only supports continuous (asked by @kpei)

Use NUTS when you can, it scales reasonably in higher dimensions. How to properly initialize it is a bit of an open question right now, but advi is probably the best we have at the moment: it gives a very good scaling matrix for NUTS. And in general, the estimation from ADVI is quite close to the posterior from MCMC (you can even get better result using FullRankADVI)
[EDIT]: We now have a new default for initialising NUTS, which is close to what STAN does by adapting the mass matrix during tuning.

[junpenglao]

• I see a flat posterior distribution using ADVI, what is happening?

This is usually an indication that the gradient is zero for that random variable. It happens if some of the manipulation of the random variables breaks the smooth geometry property. Even if all your random variables are continuous, some operation will break the gradient, for example:

. casting a variable to int: c = a.astype('int32')
. using switch a = pm.math.switch(tau >= X, a1, a2)

Sometimes it even cause unexpected consequences, for example, poor initialization for NUTS (e.g., see here).

One of the workaround is reparametrization, for example, instead of using switch:

tau = pm.Uniform('tau', lower=0, upper=100)
a = pm.math.switch(tau >= X, a1, a2)

Appoximated the switch point with a Sigmoid function

tau = pm.Uniform('tau', lower=0, upper=100)
weight = tt.nnet.sigmoid(2 * (X - tau))
a = weight * a1 + (1 - weight) * a2

[junpenglao]

• Help my model is slow! It takes hours to sample!

anwsered by @aseyboldt (see Possible Performance Issue Using Sparse Matrices in Bipartite Noisy-Or Model)
The first thing I do in cases like that is usually to try and find out why it is slow. There are two basic reasons for the slowness. Either the gradient evaluations take a lot of time or we need a lot of gradient evaluations for each sample. If it is the first, the problem is the shape of the posterior, so you need to look into reparametrizations that reduce funnels or correlations. If it is the second, you can use theano’s profiling tools to figure out where the problem is.

You can access the profile using something along those lines:

func = model.logp_dlogp_function(profile=True)
func.set_extra_values({})
x = np.random.randn(func.size)
%timeit func(x)

func.profile.summary()

You can access the number of grad evals using trace.tree_size. There are also a couple of other diagnostics that you could look at, there is a list in the doc for pymc3.NUTS.

Also check out this reply, which includes some more information: PyMC3 slows rapidly with increasing numbers of parameters

[junpenglao]

• Getting NaN or Inf during NUTS or VI (e.g., 'Bad initial energy: inf. The model might be misspecified., NaN occurred in optimization.)

Get nan or inf from model.logp (model.test point) is an attestation of incorrectly configured model?

A few thoughts: 1, In your model, the observed part you are doing: with model: ... mz = pm.Uniform("mz", lower=0, upper=1, observed=mz_freq) M = pm.Binomial("M", n=Yp, p=mz, observed=loss_shifte However, the node mz is not doing anything, since there is no prior on mz, in effect, this is equivalent to mz = mz_freq 2, Since we got -inf from model.logp (model.logp(model.test_point)), we want to identify which RV node is causing the problem. @aseyboldt suggest doing model.bijection.…

NaN occurred in optimization with NUTS

This is an issue I sometimes has as well, a few example on Github: FloatingPointError: NaN occurred in optimization with NUTS #2272 Issues getting variable in hierarchical model to update #2271 My first instinct is usually: model misspecification. For example, commenting on issue #2271: I would check the logp of each random variable and the model mu.logp(model.test_point) model.logp(model.test_point) Make sure that no logp is inf. Another quick try, is to reduce the learning rate (1e-3 i…

[junpenglao]

• NaN occurred in optimization at the first iteration in VI (i.e., seeing something like 0%| | 0/5000 [00:00<?, ?it/s])

Check the initial approximation parameters: approx.approx.params[0].eval(). If it returns nan, that is what causing the error.

Setting the value manually to zero could solve this. Below is an example if you are doing ADVI:

mu_approx = approx.approx.params[0]
rho_approx = approx.approx.params[1]

mu_approx.set_value(np.zeros(mu_approx.eval().shape))
rho_approx.set_value(np.zeros(rho_approx.eval().shape))

[junpenglao]

• How can I get the Model logp?

Evaluate logposterior at sample points

Is there a convenient / idiomatic way to evaluate logposterior (bonus: logprior/likelihood) at sample points? (I also asked the same here)

How to save logp values during sampling?

I asked a similar question here, but this time I’d like know if there’s a way to add the logp values of accepted points to the trace as we sample, so that we don’t have to evaluate it twice?

[junpenglao]

This comes up a couple time:

How to do sample_ppc after fitting the model that contains pm.Minibatch?

Answer: you cannot do that directly. But there is a workaround.

Evaluating minibatch trained model on test set

see the comment by @ferrine in this issue. Copied below: To sample ppc after minibatch inference is the following: Initialize model with shared data, not minibatch. Create minibatches for inference using pm.Minibatch Use more_replacements argument in fit like approx = pm.fit(more_replacements={full_x: minibatch_x, full_y:minibatch_y}, ...) After that you can get trace from approximation and pass it to sample_ppc. Note that it will perform ppc for examples you defined in step 1, not miniba…

[adam]

My model was also very slow. I got the following error: “g++ not detected ! Theano will be unable to
execute optimized C-implementations (for both CPU and GPU) and will default to
Python implementations.”
What solved the problem for me was to install libpython (conda install mingw libpython).

[junpenglao]

• The input of DensityDist is not just the observed, but the dictionary including potentially the parameter and the observed:

NOTE: This is no longer correct for PyMC > 4.0

How to set up a custom likelihood function for two variables

Your set up is not quite right in both cases. Put it this way, your logp should be evaluted based on the input of both the observed data and the parameters. So if you think of your logp as a function, it should has (f_1, f_2, theta_1, theta_2) as inputs. The differences between Potential and Densitydist is that, in Potential you dont distinguish between data and parameters (i.e., everything is just input to the logp), whereas in Densitydist you set up the parameterization, and evaluate on the da…

[junpenglao]

How to track a 'nan energy'?

The energy problem, if it is not from invalid start value (i.e., model.test_point) causing non-finite logp, it is usually due to gradient being non-finite. It could be difficult to diagnose, so here would be all the possible step to identify the problem: with pm.Model() as model: # your model definition # make sure all test_value are finite print(model.test_point) # make sure all logp are finite model.check_test_point() with model: step = pm.HamiltonianMC() q0 = step._logp_dlogp_fun…

[junpenglao]

• How can I marginalized a model with latent discrete variables?

Problem with hierachical occupancy model

I would do it this way: First reshape everything so within the pm.Model we are not dealing with dict: index = [] W_r = [] y_r = [] for i in range(n): index.append(np.ones((W[i].shape[0], 1))*i) W_r.append(W[i]) y_r.append(y[i][:, np.newaxis]) index = np.concatenate(index, axis=0).astype(int64) W_r = np.concatenate(W_r, axis=0) y_r = np.concatenate(y_r, axis=0) Then I will use matrix multiplication instead of nested for loops: with pm.Model() as m: beta = pm.Normal('beta', mu=…

Naive Bayes model with PyMC3

So here is how I would approach it: First, I would rewrite the explicit model, so that the observed is a flatten array - it is easier to handle shape wise: # flatten and index data data_flatten = np.reshape(data, data.shape[0]*data.shape[1]) data_index = np.repeat(np.arange(data.shape[0]), data.shape[1]) with pm.Model() as model: # Global topic distribution theta = pm.Dirichlet("theta", a=alpha) # Word distributions for K topics phi = pm.Dirichlet("phi", a=beta, shape=(K,…

[ricardoV94]

• I am getting shape errors when doing posterior predictive sampling after calling set_data

By default observed variables take their shape from the shape of the observed data. This is a convenience, so users don’t have to specify shape manually (via shape argument or dims: see Distribution Dimensionality — PyMC 5.5.0 documentation)

import pymc as pm

with pm.Model() as m:
  x = pm.MutableData("x", [0.0, 1.0, 2.0])
  data = pm.ConstantData("data", [0.0, 2.0, 4.0])

  beta = pm.Normal("beta")
  mu = beta * x
  y = pm.Normal("y", mu, observed=data)
  # This is equivalent to
  # y = pm.Normal("y", mu, observed=data, shape=data.shape)

  idata = pm.sample()

with m:
  pm.set_data({"x": [0.0, 1.0, 2.0, 3.0]})
  pm.sample_posterior_predictive(idata)

ValueError: shape mismatch: objects cannot be broadcast to a single shape.  Mismatch is between arg 0 with shape (3,) and arg 1 with shape (4,).

This happens because y has a fixed shape of (3,).
The error can be obtained directly from numpy like this:

import numpy as np

np.random.normal([0.0, 1.0, 2.0, 3.0], size=(3,))

The recommended fix is to specify how the shape of y depends on the parameter we intend to change

import pymc as pm

with pm.Model() as m:
  x = pm.MutableData("x", [0.0, 1.0, 2.0])
  data = pm.ConstantData("data", [0.0, 2.0, 4.0])

  beta = pm.Normal("beta")
  mu = beta * x
  y = pm.Normal("y", mu, observed=data, shape=mu.shape)

  idata = pm.sample()

with m:
  pm.set_data({"x": [0.0, 1.0, 2.0, 3.0]})
  pm.sample_posterior_predictive(idata)

Examples available in the docs: pymc.set_data — PyMC 5.5.0 documentation

[ricardoV94]

• Multi-core sampling is very slow / using more cores than specified.

PyMC graphs can invoke computations from other libraries (such as scipy, numpy, PyTensor) that are natively multi-threaded by default. This can cause the number of cores used to be larger than the number of cores specified in sample and may slowdown sampling altogether. Calculating Deterministics on the main thread can also use cores beyond those needed by the individual chains.

See more details and tips about how to control this behavior in the following threads:

Regarding the use of multiple cores

Maybe it helps to explain a bit how the parallelization works: When you specify cores>1 in pm.sample, pymc will start one new process for each chain. The main process then tells cores of those processes to start sampling, the others will just wait and do nothing. When one of the processes is finished, one of the waiting processes is told to start sampling. There will never be more than cores processes working at the same time. If you have n cores in your computer, it makes sense for most models…

Parallelizing chains with custom likelihood on multiple cores

I think the original model is so small, that the parallelization doesn’t do anything useful. If we try something a bit more costly it does: with pm.Model() as model: x = pm.Normal("x", shape=(50, 50)) y = pt.linalg.solve(x + x.T, np.ones(50)).sum() pm.Normal("y", mu=y, observed=100.) times = [] for cores in [1, 4, 8]: start = time.time() with model: pm.sample(cores=cores, chains=cores, tune=100, draws=100, compute_convergence_checks=False) end = time.time() …