How to quantify the confidence on a Bayesian posterior probability?

Question

Consider a physical system that depends on a parameter $0\leq \phi <\infty$. I want to (i) find the probability that this parameter is smaller than a critical value: $\phi\leq \phi_c$, and (ii) quantify my confidence on the probability that I produce in (i)---see below. To this aim, I use Bayesian approach. A priori, the probability distribution of this parameter is $p(\phi)$. I measure some observable, say $X$. By repeating this measurement $m$ times, I collect the dataset of outcomes of $X$, which I show with the vector ${\bf x} = \{x_1,x_2,\dots,x_m\}$. From the Bayes rule, I can update the probability distribution $p(\phi) \to p(\phi|{\bf x})$ as follows: $$ p(\phi|{\bf x}) = \frac{p({\bf x}|\phi)p(\phi)}{\int d\phi ~ p({\bf x}|\phi)p(\phi)}.$$ Thus, the answer to (i) is simply $p(\phi\leq\phi_c) = \int_0^{\phi_c} d\phi~ p(\phi|{\bf x})$. It remains to answer (ii), i.e., to quantify how reliable is our answer to (i). More precisely, I want to quantify the error made in estimating $p(\phi\leq\phi_c)$. Such error should depend on the number of measurements $m$, and vanish as $m\to \infty$ (most likely it scales with $1/\sqrt{m}$).

I would appreciate any hint or references on quantifying the described error. I hope that I was clear enough about the problem...if not, please let me know.

Answer 1

As you stated, using the data that you have, you can obtain a posterior distribution and calculate the integral $$I = \int_0^{\phi_c} p(\phi \mid \boldsymbol{x}) d\phi$$ explicitly. This doesn't come with an uncertainty estimate, but what we would like to know is: thinking about the data sets that could have been drawn, what would our estimates for $I$ look like? We want to use bootstrapping here. Bootstrapping involves repeatedly resampling the data, calculating an estimate for $I$ under each resampling, and then aggregating our estimates for $I$ at the end. The distribution of these estimates of $I$ allows us to estimate our uncertainty in the estimate.

You can use the Bayesian bootstrap specifically (there is a good description of it here) to approximate $I$ in a Bayesian fashion.

I'll make a specific example here. Suppose that the prior for $\phi$ is a gamma distribution with rate shape $\alpha$ and rate $\beta$ (I'll use $\alpha = \beta = 1$ in my analysis), so $$p(\phi) = \frac{\beta^{\alpha - 1}}{\Gamma(\alpha)}\phi^{\alpha - 1}e^{-\beta\phi}, \phi > 0$$ The likelihood will be Poission with parameter $\phi$, so the likelihood of $\boldsymbol{x}$ given $\phi$ is $$\frac{e^{-\phi m} \phi^{\sum_i x_i}}{\prod_i x_i !},$$ where each $x_i \in \{0, 1, 2, \dots\}.$

Hence the posterior distribution is a gamma distribution: $$p(\phi \mid \boldsymbol{x}) = \frac{(\beta + m)^{\alpha + \sum_i x_i - 1}}{\Gamma(\alpha + \sum_i x_i)} \phi^{\alpha + \sum_i x_i - 1}e^{-(\beta + m)\phi}, \phi > 0,$$ and we can calculate $I$ using lookup tables.

I then simulated this setup with $\phi = 0.95$ and $\phi_c = 1$ (so $I$ should get close to $1$ as we get more information about $\phi$). Running this in R with $m = 10$ I got an estimate of $I = 0.54$. Running $1{,}000$ resamples using the bayesboot library, I got the following output: We can see that the distribution is heavier towards $1$, but it's reasonably uniform. Now, using a sample of size $m = 100$ I got an estimate of $I = 0.32$, and the following BB output: The distribution has not changed much. Now I try running this with $m = 10{,}000$ data points. I got an estimate of $I = 0.96$ and the following BB output: Now we are really certain that $I$ is close to $1$!

The estimates have the properties that you suggested they would: they depend on the sample size $m$ and the uncertainty decreases as $m$ gets larger. I don't think you will be able to find a general formula for how quickly the uncertainty will decrease. For example, suppose that the data were generated from a uniform $U(0, \phi)$ distribution. Then if we ever observe some $\phi_c < x_i$ then we know that $P(\phi < \phi_c) = 0$ immediately, since it's impossible to observe $\phi < \phi_c < x_i$.

Here is the R code that I used to generate the plots:

install.packages("bayesboot")
library(bayesboot)
phi = 0.95
priorShape = 1
priorRate = 1
postProbLessThanOne = function(x) {
  postShape = priorShape + sum(x)
  postRate = priorRate + length(x)
  postProb = pgamma(1, shape = postShape, rate = postRate)
  return(postProb)
}

set.seed(1)
x10 = rpois(10, phi)
postProbLessThanOne(x10) # 0.54
bootstraps10 = bayesboot(data = x10, statistic = postProbLessThanOne, R = 1000,
                         R2 = 10, use.weights = FALSE, .progress = "text")
hist(bootstraps10$V1, main = "BB Distribution of I: m = 10", xlab = "Bootstrap Values of I")

x100 = rpois(100, phi)
postProbLessThanOne(x100) # 0.32
bootstraps100 = bayesboot(data = x100, statistic = postProbLessThanOne, R = 1000,
                          R2 = 100, use.weights = FALSE, .progress = "text")
hist(bootstraps100$V1, main = "BB Distribution of I: m = 100", xlab = "Bootstrap Values of I")

x10000 = rpois(10000, phi)
postProbLessThanOne(x10000) # 0.96
bootstraps10000 = bayesboot(data = x10000, statistic = postProbLessThanOne, R = 1000,
                            R2 = 10000, use.weights = FALSE, .progress = "text")
hist(bootstraps10000$V1, main = "BB Distribution of I: m = 10,000", xlab = "Bootstrap Values of I")