Unexpected behavior in Bootstrap confidence interval in Mathematica [closed]

EDIT: I think this is expected. There was a bug in my python code

I'm computing a bootstrap confidence interval using Mathematica.

f[p_, nIter_, m_] := Module[{bootstrap, samples, avg, deltas},
samples  = RandomVariate[BernoulliDistribution[p], m];
bootstrap  = RandomVariate[BernoulliDistribution[p], {nIter, m}];
avg = Mean[samples];
deltas = ArrayReduce[Mean, bootstrap, 1] - avg;
{avg - ArrayReduce[Quantile[#, 0.975] &, deltas, 1],
avg - ArrayReduce[Quantile[#, 0.025] &, deltas, 1]}]

ListPlot[
Table[Count[Table[Between[p, N[f[p, 100, 5]]], 10], True] / 10 , {p,
0, 1, 0.01}]]

This is showing that the bootstrap CI for a bernoulli mean estimate with 100 iterations from a sample of size 5 has a coverage of close to 0 for most $$p$$. But for a few $$p$$ regularly spaced (notably 0.2, 0.4, 0.6, 0.8) we have much higher coverage. This is unexpected since I would expect the bootstrap to have a coverage near 0 for all $$p$$.

I have implemented the same code in python and it does not exhibit this unusual behavior:

def f(p, n_iter, N):
samples = np.random.rand(N) < p
bootstrap = np.random.rand(n_iter, N) < p
avg = samples.mean()
deltas = bootstrap.mean(0) - avg
return (avg - np.quantile(deltas, 0.975), avg - np.quantile(deltas, 0.025))

[np.mean([p in f(p, 100, 5) for i in range(100)]) for p in np.arange(0, 1, 0.01)]

gives all values near 0.

Could it be related to some numerical precision issue?

EDIT: the correct python code is:

[np.mean([np.prod(np.sign(np.array(f(p, 100, 5)) - p)) ==-1 for i in range(100)]) for p in np.arange(0, 1, 0.01)]

which gives a result very similar to the mathematica one

• This is not an answer to your question but it appears that neither your Mathematica nor your python code are performing bootstrap sampling. A bootstrap sample would be a random sample taken with replacement from samples and not a new sample from a Bernoulli distribution.
– JimB
Nov 12 '21 at 18:00
• Good point :) But even if I change the bootstrap line to bootstrap = RandomChoice[samples, {nIter, m}] I still have the same issue (although the python code exhibits the same behavior)
– dmh
Nov 12 '21 at 18:18
• I think you've made some logic mistakes when attempting to apply a bootstrap procedure: (1) As already mentioned you need to use RandomChoice, (2) You're generating a "summary" figure that's based on a single (and unidentified) sample for each value of p, and (3) (most importantly) nIter=100 bootstrap samples should end up with 100 values but you end up with only 5 in deltas because you need ArrayReduce[Mean, bootstrap, 2] rather than ArrayReduce[Mean, bootstrap, 1]. With all of that I'm voting to close the question (but I'll retract that vote if you convince me otherwise).
– JimB
Nov 12 '21 at 21:11
• That's right, the implementation is not correct. If I change the reduce axis to 2 then I get a result that makes more sense. Thank you! and yes I agree the question should be closed since I had the same problem in python
– dmh
Nov 12 '21 at 21:40