Vrftsjtryk

Consider a physical system that depends on a parameter $0leq phi <infty$. I want to (i) find the probability that this parameter is smaller than a critical value: $phileq 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 dphi ~ p({bf x}|phi)p(phi)}.$$
Thus, the answer to (i) is simply $p(phileqphi_c) = int_0^{phi_c} dphi~ 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(phileqphi_c)$. Such error should depend on the number of measurements $m$, and vanish as $mto 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.

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}) dphi$$
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^{-betaphi}, 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")