# Confidence intervals when using Bayes' theorem

I'm computing some conditional probabilities, and associated 95% confidence intervals. For many of my cases, I have straightforward counts of x successes out of n trials (from a contingency table), so I can use a Binomial confidence interval, such as is provided by binom.confint(x, n, method='exact') in R.

In other cases though, I don't have such data, so I use Bayes' theorem to compute from information I do have. For example, given events $a$ and $b$:

$$P(a|b) = \frac{P(b|a) \cdot P(a)}{P(b)}$$

I can compute a 95% confidence interval around $P(b|a)$ using $\textrm{binom.confint}(\#\left(b\cap{}a),\#(a)\right)$, and I compute the ratio $P(a)/P(b)$ as their frequency ratio $\#(a)/\#(b)$. Is it possible to derive a confidence interval around $P(a|b)$ using this information?

Thanks.

-
$a$ and $b$ are events. In my case, $a$ is a system failure (which is quite rare, so relatively hard to find "in the wild"), and $b$ is a pre-failure alarm, so I'm measuring the probability of failure given an alarm. –  Ken Williams Aug 21 '12 at 18:56
The above comment was in response to someone who asked for more background on what $a$ and $b$ were, but seems to have deleted that comment. –  Ken Williams Aug 21 '12 at 19:35
Well you can't just take the confidence interval for p(b|a) and scale it by p(a)/p(b) because of the uncertainty in the estimate of that ratio. If you can construct a 100(1-α)% confidence interval for p(a)/p(b) call it [A, B] then take the lower bound for a 100(1-α)% confidence interval for p(b|a) and multiply it by A and take the upper bound for p(b|a) and multiple it by B. That should give at an interval that has at least a 100(1-α)$^2$% confidence level for p(a|b). –  Michael Chernick Aug 21 '12 at 19:37
Could work... getting a confidence interval for $P(a)/P(b)$ isn't obvious to me though - do you feel like moving this into the "Answer" area? I promise at least one upvote. =) –  Ken Williams Aug 21 '12 at 21:46
Don't you want a Bayesian credible interval instead? That is directly computable from the posterior distribution of $a$. –  whuber Aug 21 '12 at 22:44
show 1 more comment

Well, you can't just take the confidence interval for $p(b|a)$ and scale it by $p(a)/p(b)$ because of the uncertainty in the estimate of that ratio. If you can construct a $100(1-\alpha)\%$ confidence interval $[A, B]$ for $p(a)/p(b)$, then take the lower bound for a $100(1-\alpha)\%$ confidence interval for $p(b|a)$ and multiply it by $A$ and take the upper bound for $p(b|a)$ and multiply it by $B$. That should give at an interval that has at least a $100(1-\alpha)^2\%$ confidence level for $p(a|b)$.
That seems workable at least as a first stab. But I'm not aware of a method for deriving confidence intervals for the ratio of two Bernoulli probabilities $P(a)/P(b)$ given observed counts of $a$ and $b$ in a sample population. –  Ken Williams Aug 21 '12 at 22:33
There's one purported method for computing an interval on $P(a)/P(b)$ by Katz et al., 1978. I can't find the original paper without paywalls, but this citer seems to show the method: jstor.org/discover/10.2307/2531405 –  Ken Williams Aug 22 '12 at 15:59
If I haven't erred, here's a function to do that ratio-of-Bernoullis interval estimation: binrat.confint <- function(x, y, n=Inf, m=n, p=0.95) { s2 <- 1/x - 1/n + 1/y - 1/m; x/y * exp(c(-1:1)*pnorm((1+p)/2)*sqrt(s2)) }  –  Ken Williams Aug 24 '12 at 19:05