Variance and Parameter Ranking in Bayesian Inference

[jessegrabowski]

@gil_kanfer

When I tried to execute the “dummy” model, it was flagging that other parameters like sigma and Mu were missing. The model started to function when I assigned pm.FLAT to these parameters, but I’m not completely sure why this step is necessary. My understanding is that the values of these parameters are already inferred from the fitted model, so why do we need to allocate flat priors?

I should have linked you this blog post that goes deeper into how this dummy stuff works. I also tried to explain it here, maybe that will help?

The dummy model needs to have all information on how to compute stuff up to the nodes you specify in var_names. In my example, the var_names was only beta_dist, which needs to know beta, which (looking at the DAG) needs to know beta_bar, sigma_bar, and z_beta. If you instead asked for a variable deeper into the dag, like count, you’d have to specify everything upstream of that, which would include sigma and mu.

@gil_kanfer

Additionally, I have been struggling to decipher the RV structure in the PPA. Usually, I run the model and then utilize beta.eval().shape to comprehend the dimensions, but this approach seems to fail when it comes to PPA analysis. As a result, I don’t understand why we use beta[:, None]; from the fitted model, beta has 20 sgRNAs, 1000 draws, and 4 chains.

While the variables are just pytenor symbols on a computation graph – that is, before you compile and run the model – you don’t have to reason about the batch dimension (chains and draws). You just have to think about the core dimensions. So on the graph, beta is just a (20,) array, so beta[:, None] is a (20, 1) column vector. Pytensor’s whole raison d’etre is to then handle the vectorization. Just focus on the shapes as you would understand them if you were writing out your model on a piece of paper. Or better yet, use dims absolutely everywhere, and reason about the named dimensions instead of the shapes.

Actually there’s nothing special about sample_posterior_predictive. If all you wanted were pairwise differences between the betas, you can just do it with the (4, 1000, 20) posterior with something like idata.posterior.beta.values[..., None] - idata.posterior.beta.values[None, ...]. Here’s an example where I first tried to do this, then realized it was dumb and overly complicated and did it elegantly with sample_posterior_predictive in a dummy model.

@gil_kanfer

In an attempt to calculate beta_dist on a gene level, how do I integrate the gene indexing vector into the post-predictive analysis?

Exactly the same way you put it into your model. It should be a 1-to-1 map. Copy+paste everything you need up until you compute beta in your primary model, then use that beta to compute some statistic of interest (beta_dist in my example) , then pm.sample_posterior_predictive(idata, var_names=[statistic_of_interest]). If you get stuck I’m happy to look at specific code.

As for size, you’re just doing pairwise subtraction so it should scale pretty well. If it gets really bad you could try to think about how to compute only the unique pairs with some clever indexing, since the distance matrix is symmetrical.