def bs_tt_bblock (times, tstarts, tstops, p0=0.05, nbootstrap=512):
    """Bayesian Blocks for time-tagged events with bootstrapping uncertainty
    assessment. THE UNCERTAINTIES ARE NOT VERY GOOD! Arguments:

    tstarts        - Array of input bin start times.
    tstops         - Array of input bin stop times.
    times          - Array of event arrival times.
    p0=0.05        - Probability of preferring solutions with additional bins.
    nbootstrap=512 - Number of bootstrap runs to perform.

    Returns a Holder with:

    blockstarts - Start times of output blocks.
    bsrates     - Mean event rate in each bin from bootstrap analysis.
    bsrstds     - ~Uncertainty: stddev of event rate in each bin from bootstrap analysis.
    counts      - Number of events in each output block.
    finalp0     - Final value of p0, after iteration to minimize `nblocks`.
    ledges      - Times of left edges of output blocks.
    midpoints   - Times of midpoints of output blocks.
    nblocks     - Number of output blocks.
    ncells      - Number of input cells/bins.
    origp0      - Original value of p0.
    rates       - Event rate associated with each block.
    redges      - Times of right edges of output blocks.
    widths      - Width of each output block.

    """
    times = np.asarray (times)
    tstarts = np.asarray (tstarts)
    tstops = np.asarray (tstops)

    nevents = times.size
    if nevents < 1:
        raise ValueError ('must be given at least 1 event')

    info = tt_bblock (tstarts, tstops, times, p0)

    # Now bootstrap resample to assess uncertainties on the bin heights. This
    # is the approach recommended by Scargle+.

    bsrsums = np.zeros (info.nblocks)
    bsrsumsqs = np.zeros (info.nblocks)

    for _ in range (nbootstrap):
        bstimes = times[np.random.randint (0, times.size, times.size)]
        bstimes.sort ()
        bsinfo = tt_bblock (tstarts, tstops, bstimes, p0)
        blocknums = np.minimum (np.searchsorted (bsinfo.redges, info.midpoints),
                                bsinfo.nblocks - 1)
        samprates = bsinfo.rates[blocknums]
        bsrsums += samprates
        bsrsumsqs += samprates**2

    bsrmeans = bsrsums / nbootstrap
    mask = bsrsumsqs / nbootstrap <= bsrmeans**2
    bsrstds = np.sqrt (np.where (mask, 0, bsrsumsqs / nbootstrap - bsrmeans**2))
    info.bsrates = bsrmeans
    info.bsrstds = bsrstds
    return info