Sorry, I pasted a reply meant for another thread.
I was going to go with bootstrap sampling from each set of posterior means - does that make sense?
You should always work directly with the posterior samples.
Sounds like you’re doing the 2nd bullet point from the scenarios I mentioned.
This should do the trick:
pyrff.sampling_probabilities([
samplesA, # 👈 NumPy array with shape (chains*draws,)
samplesB, # also a 1D array; doesn't need to be the same shape though
], correlated=False)