# Estimating Bernoulli success rate (and confidence interval) from success count in sample?

I have a sample of 0s and 1s from a Bernoulli distribution with unknown success rate $p$. I know the function FindDistributionParameters can be used to estimate $p$:

FindDistributionParameters[{1, 0, 0, 1, 0, 0, 0}, BernoulliDistribution[p]]

But how can I get confidence intervals?

• Talke a look at DistributionFitTest Commented Dec 16, 2015 at 16:09
• DistributionFitTest computes $p$-values, but I don't think it computes confidence intervals.
– a06e
Commented Dec 17, 2015 at 0:03
• Yep. You're right. Sorry Commented Dec 17, 2015 at 0:28
• @belisariushassettled Nevermind. I did not know about DistributionFitTest, and learning about it was useful, so thanks for pointing it out.
– a06e
Commented Dec 17, 2015 at 13:22

Just in case it is of interest. In this toy example "coin with probability of heads 0.4"...using a uniform prior and conjugacy of Beta distribution. In the following the red grid line is correct value of p, the blue line is uniform distribution, the orange line is the posterior distribution and the other gridlines the quantiles 0.025 to 0.0975

Manipulate[
sample = RandomVariate[BernoulliDistribution[0.4], n];
b = BetaDistribution[1 + Count[sample, 1], Count[sample, 0] + 1];
m = NExpectation[x, x \[Distributed] b];
ci = x /. First@Quiet[NSolve[CDF[b, x] == #, x]] & /@ {0.025, 0.975};
Plot[{PDF[BetaDistribution[1, 1], x], PDF[b, x]}, {x, 0, 1},
PlotRange -> All, Frame -> True,
GridLines -> {{ci[[1]], m, ci[[2]], {0.4, {Red, Thick}}}, None}],
{n, {10, 100, 1000, 10000}},
Button["new sample",
sample = RandomVariate[BernoulliDistribution[0.4], n]]]

Using the example of the common "95%" confidence interval there are a variety of methods to construct approximate 95% confidence intervals and one "exact" method (although the exact method is conservative and guarantees "at least 95% confidence" as opposed to "approximately 95% confidence").

See https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval for more details. Here is some code to implement the exact (Clopper-Person) interval and the most basic approximation for equal-tail area confidence intervals:

data = {1, 0, 0, 1, 0, 0, 0};
x = Total[data]; (* Number of successes *)
n = Length[data]; (* Number of trials *)
alpha = 0.05;

(* Clopper-Pearson interval - At least 100(1-alpha)% confidence *)
lower = InverseCDF[BetaDistribution[x, n - x + 1], alpha/2]
(* 0.036692566176085545 *)
upper = InverseCDF[BetaDistribution[x + 1, n - x], 1 - alpha/2]
(* 0.7095791362626571 *)

(* Simple normal approximation - Approximately 100(1-alpha)% confidence *)
lower = Max[0, x/n - 1.96 (x (n - x)/n^3)^0.5]
(* 0 *)
upper = Min[1, x/n + 1.96 (x (n - x)/n^3)^0.5]
(* 0.6203782963279159 *)

Update

As pointed out by @becko, the above Clopper-Pearson confidence interval formula only works when $x \neq 0$ and $x \ne n$. Below is an adjustment to include those extreme values. Essentially when $x = 0$ or $x = n$, then I turn things into a one-tailed confidence interval.

clopperPearson[n_, x_, alpha_] :=
If[x == 0, {0, InverseCDF[BetaDistribution[1, n], 1 - alpha]},
If[x == n, {lower = InverseCDF[BetaDistribution[n, 1], alpha], 1},
{InverseCDF[BetaDistribution[x, n - x + 1], alpha/2],
InverseCDF[BetaDistribution[x + 1, n - x], 1 - alpha/2]}]]

clopperPearson[7, 0, 0.05]
(* {0,0.34816365513116077} *)
clopperPearson[7, 1, 0.05]
(* {0.0036102968619006193,0.5787231970431952} *)
clopperPearson[7, 2, 0.05]
(* {0.036692566176085545,0.7095791362626571} *)
clopperPearson[7, 3, 0.05]
(* {0.09898827844250789,0.8159484323599169} *)
clopperPearson[7, 4, 0.05]
(* {0.184051567640083,0.9010117215574921} *)
clopperPearson[7, 5, 0.05]
(* {0.29042086373734266,0.9633074338239145} *)
clopperPearson[7, 6, 0.05]
(* {0.4212768029568048,0.9963897031380994} *)
clopperPearson[7, 7, 0.05]
(* {0.6518363448688391,1} *)

I say "I turn things into a one-tailed confidence interval" with an emphasis on "I" because it seems most (if not all) of the online Clopper-Pearson confidence limit calculators leave just alpha/2 in one tail and zero in the other when $x=0$ or $x=n$. That seems too conservative to me which is why I changed one of the tail areas to alpha. But I might be in the minority doing so (and maybe wrong in some sense).

• The Clopper-Pearson blows up in the extreme cases x == n or x==0. Can you fix this?
– a06e
Commented Dec 16, 2015 at 18:59