import numpy as np

def coverage(a,cl=68):
    ''' For an array a return width for which cl percent of the distribution is covered'''
    b = np.percentile(abs(a - np.median(a)), cl)
    return b



#
# Example application of the estimator, using 10000 random Gaussian events
#

Nevt = 10000
g = np.random.normal(0.,1.0,Nevt)

#
#  Add estimation of statistical uncertainties using bootstrap method.
#  Namely, resample the distribution of the generatted events
#  using Possion distribution  with mu=1.  Each event is then taken
#  0, 1, 2, ... N times according to Poisson(1) probability.
#  For more details see
#
#    B. Efron, Bootstrap Methods: Another Look at the Jackknife, Ann. Statist. 7 (1979) 1.
#
#  In case of histogram based input, use Poisson(Nentries) where Nentries
#  is the number of events in each bin for re-sampling. 
#

Nrs = 100

# prepare bootstrap weights:
bs = np.random.poisson(1.,(Nevt,Nrs))

# Loop over bootstrap replica, estimate 68% and 95% coverage for them:
p68 = []
p95 = []
for i in range(Nrs):
    resampled = np.repeat(g,bs[:,i])    # take event 0,1, ... N times following Poisson(1)
    p68.append(coverage(resampled))
    p95.append(coverage(resampled,95.))

# Direct:
print ("68% Coverage of g:", coverage(g))
print ("95% Coverage of g:", coverage(g,95.))

# Mean over bootstrap +- rms:
print ("68% coverage, mean +- rms:", np.mean(p68),"+-",np.std(p68))
print ("95% coverage, mean +- rms:", np.mean(p95),"+-",np.std(p95))