Bayesian stats need simulation
The simplest Bayesian test would be a Beta-Binomial.
N <- 10000
a_samples <- rbeta(N,2+1,13+1) # 2 success, 13 failure
b_samples <- rbeta(N,3+1,11+1) # 3 success, 11 failure
sum(b_samples > a_samples)/N # 0.7028
To get the probability we are essentially doing a Monte-Carlo simulation by drawing from two Beta distributions.
Finding an analytical solution would entail solving the following:
\[P(B > A) = P(B - A > 0) = \int_{x=0}^{\text{max}}\text{pdf}_{X}(x)\]Where the pdf can’t be Beta because the domain is [-1, 1] and not [0, 1]. Without knowing the exact form of this distribution we are unable to solve the integral analytically. So we resort to MC simulations. This is not a hack but a numerical method of solving the problem. Maybe not as accurate and efficient as solving the integral, but easy to implement and understand. And with increasing sophistication of Bayesian models the integration only gets more complicated which is why we need MCMC samplers like Stan and PyMC3.
References: https://www.countbayesie.com/blog/2020/8/16/why-bayesian-stats-need-monte-carlo-methods
Notes mentioning this note
There are no notes linking to this note.