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]}]

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

The result is: boostrap

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

  • 1
    $\begingroup$ 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. $\endgroup$
    – JimB
    Nov 12 '21 at 18:00
  • $\begingroup$ 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) $\endgroup$
    – dmh
    Nov 12 '21 at 18:18
  • $\begingroup$ 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). $\endgroup$
    – JimB
    Nov 12 '21 at 21:11
  • $\begingroup$ 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 $\endgroup$
    – dmh
    Nov 12 '21 at 21:40