Comparing Posterior Samples - Conceptual Question

Questions
2

Hi,
if you have two models for the same dataset, you’d typically use a model comparison metric. There are some examples in Model comparison — PyMC3 3.11.4 documentation

If it’s actually the same model and you’re interested in the probability that a parameter (e.g. sales) for company A is higher than for company B, it depends:

  • If you did one pm.sample(), fitting the PyMC model to all data at the same time, the posterior draws are correlated and you can simply do p_B_greater_A = numpy.mean(posterior.sales_B > posterior.sales_A).
  • If you did multiple pm.sample(), fitting the same model to two different datasets, you’ll not have any correlation between sales_A and sales_B. Now you want to get the p_B_greater_A given let’s say 10_000 posterior draws of sales from a model fitted to company A data and say 5_000 posterior draws of sales from a model fitted to company B data. Here you can’t do an elementwise comparison, but need to consider all 10_000 * 5_000 combinations. So “What’s the probability that a random posterior sample from B is greater than a random posterior sample from A?”. Doing this exactly for N>2 candidates to compare is surprisingly difficult to implement, but you can just use my pyrff package for that: p_best = pyrff.sampling_probabilities([posteriorA.sales, posteriorB.sales], correlated=False). See the test case for some simple examples.

cheers