So, I’m trying to wrap my head around how to to the following:
I’ve generated two prediction distributions, P1 & P2. Now I want to get a new distribution that’s basically P1/P2.
In the frequentist way I would just have divided the means. But in a Bayesian fashion, I’d like to have the distributions instead.
One way I can think of is just to go row by row in each prediction array and do division. But this doesn’t feel like the right way to do it.
Let’s say X \sim P_1 and Y \sim P_2 with the range of possible values identical between the two. Do you want to generate a distribution for X/Y or are you actually interested in the ratio of densities P_1(x) / P_2(x)? The second case is not going to be a distribution in general. It seems like you are more interested in the former. In that case, if you have posterior samples of the predictive distributions for X and Y (e.g. from the trace), you can indeed just take the ratio of the samples. This is one of the nice things about Monte Carlo - we can just treat functions of the sampled values as sampled values of the functions.
You are correct, it was the former. Thanks!
Of course, for completeness, should you want an exact expression and sample directly from the resulting distribution, that is also possible, albeit with substantially more effort. Typically the expression for the quotient of two random variates (and/or its asymptotic expansion) is obtained through inverse Mellin transforms (see https://projecteuclid.org/download/pdf_1/euclid.aoms/1177730201 for a brief overview). You would then implement a custom distribution family for X/Y, and sample directly from that, instead of sampling from X and Y separately.