How to index pm.Data with a Rv's output -->Data[rv]?

1

I believe I might be approaching this the wrong way. I wanted to index pm.Data using the output from a random variable (RV) within the model, but the random variable appears to be stuck at a single value, indicating that it’s not sampling properly. Here’s an example of the code I used (there are also some warning generated in the logs):

import numpy as np
import pymc as pm
import pytensor.tensor as pt
import pytensor
import arviz as az
import seaborn as sns


with pm.Model() as m:
  mean_sales = pm.MutableData("mean_sales",[100, 0., 50000])
  categories = pm.Categorical("categories", p=[0.33, 0.33, 0.34])

  sales = pm.Normal("sales",mu=mean_sales[categories],sigma=1) #i <<-indexing here

  idata=pm.sample(10000)

Sampling 2 chains, 0 divergences ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00 / 0:00:20
/usr/local/lib/python3.11/dist-packages/arviz/stats/diagnostics.py:596: 
RuntimeWarning: invalid value encountered in scalar divide
(between_chain_variance / within_chain_variance + num_samples - 1) / (num_samples)

When I plot the posterior samples from the random variable categories, it appears to be stuck at a single value which has similar effects on its descendant nodes.

image
image708×506 17.3 KB

2

Sampling from a prior distribution using MCMC can often be harder than you expect. In this case, the sampler is faced with the following likelihood function:

Untitled
Untitled556×413 8.01 KB

Since the standard deviation of the Normal is so small, it essentially has 3 spikes at the three possible means, and a sea of zero everywhere else. Depending on the initial value, the sampler will find the nearest mean and wander around it. It will not be able to make a proposal that crosses the large gulf between the means, and will thus never draw samples from them. This blog post by @colcarroll has some nice illustrations that will be relevant for building intuition about this issue. This animation is also one I really like.

If you crank up the sigma to something like 10,000, you get a likelihood function that is more tractable:

Untitled
Untitled547×428 12.9 KB

And posterior samples from categories reflect this:

Untitled
Untitled555×443 8.75 KB

Another way to remedy the situation is to condition the likelihood on some generated data.

3

Yeah, Thanks for sharing Jesse. Meanwhile using pm.sample_prior_predictive instead of pm.sample