Choosing Histogram Bins for Time-Indexed Events

Consider a hypothetical service that writes one log line per request, and an engineer who wants to know when during the day traffic arrives in bursts. She loads the day’s timestamps into a histogram and has to pick a bin width. With two-hour bins the day is a gentle daytime hump. With 15-minute bins a sharp spike appears late in the afternoon, and with 3-minute bins the plot is a thicket of peaks. The same 1,286 events support all three pictures, and the plotting call gives no hint about which one describes the traffic.

A defensible answer treats the histogram as an estimate of the event rate and asks which bin width makes that estimate closest to the true rate. The true rate is unknown, but for events that arrive independently the distance can be estimated from the bin counts alone, and the same estimate can report that no width resolves anything. The derivation rests on one property of the Poisson process, that the variance of a count equals its mean, and it fails in a predictable way when events cluster. The statistical background is in probability for machine learning; the rest is built from expectation and variance.

A histogram of timestamps estimates a rate

Let events occur at times in a window \([0, T)\) and suppose they follow a Poisson process with rate \(\theta(t)\), the expected number of events per unit time near \(t\). Counts in disjoint intervals are independent Poisson variables whose mean is the integral of \(\theta\) over the interval; the glossary entry on the Poisson process states the definition and the equal-mean-and-variance property used throughout. A histogram with \(N\) equal bins of width \(\Delta = T/N\) counts \(k_i\) events in bin \(i\) and estimates the rate on that bin as

\[ \hat{\theta}(t) = \frac{k_i}{n\Delta} \quad \text{for } t \text{ in bin } i, \]

where \(n\) is the number of windows pooled into the counts: \(n = 1\) for a single day, and \(n = 7\) if the same day is observed over a week and the counts are added. Dividing by \(n\Delta\) keeps the estimate in events per minute whatever the width, so histograms with different widths estimate the same curve and can be compared with each other and with the truth. Given the number of events, the event times of a Poisson process are independent draws from the density \(\theta/\int\theta\), so tools for estimating a probability density from a sample apply. The Poisson model adds one thing: the counts carry their own noise level, which the derivation below exploits.

The examples use a synthetic day. The true rate is a constant 0.5 events per minute, plus a broad daytime hump centered at minute 780 with a standard deviation of 200 minutes, plus a short burst at minute 1013 with a standard deviation of 6 minutes. The burst adds about 75 events to roughly 1,250 in total. The day is treated as a cycle, so bin grids can be rotated past midnight; the rate at both ends of the day is nearly the same, which makes this harmless. Events are drawn by thinning: propose events at a rate \(\theta_{\max}\) that is never exceeded, then keep a proposal at time \(t\) with probability \(\theta(t)/\theta_{\max}\). The programs use Python 3.12 and NumPy 2.5.3 with fixed random seeds, and the data are generated, not collected. The first block is saved as ssbins.py; later programs start with from ssbins import *.

import numpy as np

T = 1440.0  # observation window in minutes: one day, treated as a cycle


def rate(t):
    """Hypothetical true rate in events per minute: baseline, a daily hump, a short burst."""
    t = np.asarray(t, dtype=float)
    base = 0.5
    hump = 0.9 * np.exp(-0.5 * ((t - 780.0) / 200.0) ** 2)
    burst = 5.0 * np.exp(-0.5 * ((t - 1013.0) / 6.0) ** 2)
    return base + hump + burst


def simulate(rng, rate_fn=rate, rate_max=7.0):
    """Inhomogeneous Poisson process on [0, T) by thinning."""
    n_cand = rng.poisson(rate_max * T)
    cand = rng.uniform(0.0, T, n_cand)
    keep = rng.uniform(0.0, rate_max, n_cand) < rate_fn(cand)
    return np.sort(cand[keep])

The counts function in the same module does the binning. A single day drawn with a fixed seed contains 1,286 events, and the next program prints it at three widths as rows of block characters, one character per bin, with height proportional to the rate (events per minute, full height at 6), followed by the number of local maxima at eight widths.

def counts(times, n_bins, shift=0.0):
    """Counts in n_bins equal bins of the cycle [0, T), with the bin grid shifted by `shift`."""
    t = (np.asarray(times) - shift) % T
    return np.bincount((t / (T / n_bins)).astype(int), minlength=n_bins)[:n_bins]

import numpy as np
from ssbins import *

rng = np.random.default_rng(8)
x = simulate(rng)
print("events:", len(x), "  true peak rate: %.2f per minute" % rate(1013.0))

bars = "▁▂▃▄▅▆▇█"
def spark(k, width):                         # one character per bin; height = rate / 6 per minute
    r = k / width
    return "".join(bars[min(7, int(v / 6.0 * 8))] for v in r)

for w in (120, 30, 15):
    k = counts(x, int(T / w))
    print("%4d min  %s" % (w, spark(k, w)))

def local_max(k):                            # bins strictly higher than both neighbours (cyclic)
    a = np.r_[k[-1], k, k[0]]
    return int(((a[1:-1] > a[:-2]) & (a[1:-1] > a[2:])).sum())

print("\nwidth(min)  bins  local maxima  busiest bin (events/min)")
for w in (3, 6, 12, 15, 30, 60, 120, 240):
    k = counts(x, int(T / w))
    print("%9d %6d %12d %14.2f" % (w, len(k), local_max(k), k.max() / w))

events: 1286   true peak rate: 5.96 per minute
 120 min  ▁▁▁▂▂▂▂▂▃▂▁▁
  30 min  ▁▁▁▁▁▁▁▁▁▁▁▁▂▁▂▂▁▁▂▂▂▃▂▂▂▂▃▂▂▂▂▂▂▅▂▁▂▂▂▁▁▁▂▁▁▁▁▁
  15 min  ▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂▁▁▂▂▂▂▁▁▁▂▂▂▂▃▂▂▂▃▂▂▃▂▂▂▂▂▃▂▂▂▂▂▂▂▃▂▂▂▂▂▃▆▂▂▂▁▂▁▂▂▁▂▁▁▁▁▁▁▂▁▂▁▁▂▁▁▁▁▁▁

width(min)  bins  local maxima  busiest bin (events/min)
        3    480          118           5.67
        6    240           67           5.50
       12    120           35           4.75
       15     96           30           4.47
       30     48           14           3.13
       60     24            8           2.18
      120     12            2           1.51
      240      6            1           1.30

At 120 minutes the burst is gone: the busiest bin holds 1.5 events per minute while the true peak is 5.96. At 15 minutes the spike is clear, but the bins around it wobble, and the table counts 30 local maxima where the true rate has two humps. At 3 minutes there are 118. The first behavior is oversmoothing, where averaging erases structure; the second is undersmoothing, where noise is drawn as structure. The glossary entry on oversmoothing and undersmoothing describes why both look plausible on a plot. The count of local maxima is not a selection rule, because it falls steadily with width and never reaches a minimum. Choosing among the plots by eye picks the one that matches what the viewer expected to see, so a criterion computed from the data, with a minimum, is needed.

Why the error has a minimum

Define the integrated squared error of a histogram as the average over the window of the squared distance to the true rate, and the mean integrated squared error (MISE) as its expectation over repeated days from the same process:

\[ \mathrm{MISE}(\Delta) = \frac{1}{T}\int_0^T \mathbb{E}\big[(\hat{\theta}(t) – \theta(t))^2\big]\,dt . \]

The glossary entry on mean integrated squared error introduces the quantity. At each time \(t\) the expected squared error splits into variance plus squared bias, the same identity that the bias-variance decomposition derives for predictions. Write \(\bin{\theta}_i\) for the average of \(\theta\) over bin \(i\). The bin’s count is Poisson with mean \(n\Delta\bin{\theta}_i\), so \(\mathbb{E}\hat{\theta} = \bin{\theta}_i\) and \(\operatorname{Var}\hat{\theta} = \bin{\theta}_i/(n\Delta)\). Integrating the variance over the window and dividing by \(T\) gives exactly \(r/(n\Delta)\), where \(r\) is the mean rate over the window. That term falls as bins widen, because each bin collects more events.

The bias comes from the bin averaging the rate over an interval. Inside a bin the rate differs from its average by roughly the local slope \(\theta’\) times the distance from the bin’s center. Averaging the squared distance over a bin of width \(\Delta\) gives \(\theta’^2\Delta^2/12\), so the integrated squared bias is approximately \(R\Delta^2/12\), with \(R = \frac{1}{T}\int \theta'(t)^2\,dt\) the mean squared slope. That term grows as bins widen. Together,

\[ \mathrm{MISE}(\Delta) \approx \frac{r}{n\Delta} + \frac{R\,\Delta^2}{12}, \qquad \Delta^{*} = \left(\frac{6\,r}{n\,R}\right)^{1/3}. \]

The second formula comes from setting the derivative of the first to zero. It is the classical result of Scott for histograms of a density, written in rate units. The approximation of the bias holds when the rate is close to linear inside one bin, so it degrades when the width is comparable to the scale of the sharpest feature. The formula says two useful things: the best width shrinks as \(n^{-1/3}\) when more windows are pooled, and it shrinks as the rate becomes rougher, with the roughness \(R\) dominated by the sharpest features rather than by the overall spread of the events.

Because the data here are simulated, the true rate is available and both terms can be computed exactly by averaging the true rate over each bin. These helpers live in a second module, truth.py, which exists only for simulations because it needs the true rate. It evaluates the rate on a fine grid of 57,600 points and averages the two terms over twelve rotations of the bin grid, so that no result depends on where a bin edge happens to fall.

import numpy as np
from ssbins import T, rate

NGRID = 57600                                  # fine grid for numerical integration over [0, T)
GRID = (np.arange(NGRID) + 0.5) * (T / NGRID)
TH = rate(GRID)                                # true rate on the grid
DG = T / NGRID


def bin_index(n_bins, shift=0.0):
    return (((GRID - shift) % T) / (T / n_bins)).astype(int)

def mise_parts(n_bins, n_windows=1, shifts=1):
    """Exact (integrated variance, integrated squared bias) per unit time of the histogram
    rate estimate, averaged over grid origins. Needs the true rate: simulation only."""
    delta = T / n_bins
    var, bias2 = [], []
    for s in np.arange(shifts) * delta / shifts:
        idx = bin_index(n_bins, s)
        thbar = np.bincount(idx, weights=TH, minlength=n_bins) / np.bincount(idx, minlength=n_bins)
        var.append(thbar.sum() / n_windows / T)
        bias2.append(((thbar[idx] - TH) ** 2).sum() * DG / T)
    return float(np.mean(var)), float(np.mean(bias2))

import numpy as np
from ssbins import *
from truth import *

print("width(min)  integrated variance  integrated bias^2   MISE")
for N in (6, 12, 24, 48, 72, 96, 144, 288, 480):
    v, b = mise_parts(N, shifts=12)
    print("%9.1f %18.4f %18.4f %9.4f" % (T / N, v, b, v + b))

r = TH.mean()                                    # mean rate
R = (np.gradient(TH, GRID) ** 2).mean()          # mean squared slope of the true rate
print("\nmean rate r = %.3f /min,  mean squared slope R = %.5f" % (r, R))
print("small-width approximation  delta* = (6 r / (n R))^(1/3) = %.1f min" % ((6 * r / R) ** (1 / 3)))
Ns = np.arange(2, 361)
m = [mise_exact(int(n), shifts=12) for n in Ns]
print("exact minimiser of MISE over widths T/2 .. T/360: %.1f min" % (T / Ns[int(np.argmin(m))]))
for n in (1, 8, 64):
    print("  n=%2d pooled windows: approximation %.1f min" % (n, (6 * r / (n * R)) ** (1 / 3)))

width(min)  integrated variance  integrated bias^2   MISE
    240.0             0.0036             0.1783    0.1819
    120.0             0.0072             0.1562    0.1634
     60.0             0.0144             0.1272    0.1416
     30.0             0.0288             0.0834    0.1123
     20.0             0.0433             0.0543    0.0976
     15.0             0.0577             0.0364    0.0941
     10.0             0.0865             0.0188    0.1053
      5.0             0.1731             0.0052    0.1782
      3.0             0.2885             0.0019    0.2904

mean rate r = 0.865 /min,  mean squared slope R = 0.00257
small-width approximation  delta* = (6 r / (n R))^(1/3) = 12.6 min
exact minimiser of MISE over widths T/2 .. T/360: 15.3 min
  n= 1 pooled windows: approximation 12.6 min
  n= 8 pooled windows: approximation 6.3 min
  n=64 pooled windows: approximation 3.2 min

At 3 minutes nearly all of the 0.2904 comes from variance, and at 240 minutes nearly all of the 0.1819 comes from bias. The sum is smallest near 15 minutes, at 0.0941. The small-width formula gives 12.6 minutes against the exact 15.3, the same scale and slightly narrow because the burst’s standard deviation of 6 minutes is not small compared with a 15-minute bin. With 8 pooled days the formula moves to 6.3 minutes, and with 64 to 3.2. The minimizer is a property of this rate and this amount of data, which is why no constant width can be correct for every dataset.

Why the textbook rules miss event structure

The usual default rules choose the width from the sample alone. The glossary entry on histogram bin width lists them, and NumPy’s documentation for histogram_bin_edges gives their formulas: Sturges uses \(\log_2 n + 1\) bins; Scott uses a width of \(\sigma(24\sqrt{\pi}/n)^{1/3}\), where \(\sigma\) is the standard deviation of the data; Freedman and Diaconis use \(2\,\mathrm{IQR}/n^{1/3}\), where IQR is the interquartile range; and the auto setting takes whichever of Sturges and Freedman–Diaconis gives more bins. NumPy states that the Sturges rule assumes normal data and is too conservative for larger, non-normal datasets.

Scott’s and the Freedman–Diaconis rules reproduce the \(n^{-1/3}\) scaling of \(\Delta^{*}\). What they replace is the roughness \(R\). Scott substitutes the roughness of a Gaussian with the sample’s standard deviation, and Freedman and Diaconis bound it using the interquartile range. For timestamps over a day, both statistics describe how events are spread across the window, and a narrow burst barely moves them. The next program compares three synthetic days: a daytime hump alone, the hump with the 6-minute burst, and a flash crowd whose burst has a standard deviation of 3 minutes and carries about 700 of its 1,200 events. For each it prints the average standard deviation and IQR of the event times over 50 simulated days, the Scott and Freedman–Diaconis widths computed from them, the true roughness \(R\), and the width \((6r/R)^{1/3}\) from the small-width formula.

import numpy as np
from ssbins import *

rng = np.random.default_rng(4)
hump = lambda t: 0.5 + 0.9 * np.exp(-0.5 * ((np.asarray(t) - 780.0) / 200.0) ** 2)
flash = lambda t: 0.35 + 93.0 * np.exp(-0.5 * ((np.asarray(t) - 1013.0) / 3.0) ** 2)
days = (("hump only", hump, 2.0), ("hump + burst", rate, 7.0), ("flash crowd", flash, 100.0))

g = (np.arange(200000) + 0.5) * T / 200000
print("%-13s %6s %6s %6s %12s %10s %10s %14s" %
      ("day", "events", "sigma", "IQR", "Scott width", "FD width", "R", "(6r/R)^(1/3)"))
for name, fn, rmax in days:
    xs = [simulate(rng, fn, rmax) for _ in range(50)]
    K = np.mean([len(x) for x in xs])
    sd = np.mean([x.std() for x in xs])
    iqr = np.mean([np.subtract(*np.percentile(x, [75, 25])) for x in xs])
    th = fn(g)
    r, R = th.mean(), (np.gradient(th, g) ** 2).mean()
    print("%-13s %6.0f %6.0f %6.0f %9.0f min %7.0f min %10.2g %11.1f min" %
          (name, K, sd, iqr, sd * (24 * np.sqrt(np.pi) / K) ** (1 / 3), 2 * iqr / K ** (1 / 3), R, (6 * r / R) ** (1 / 3)))

day           events  sigma    IQR  Scott width   FD width          R   (6r/R)^(1/3)
hump only       1178    350    484       116 min      92 min    2.5e-06       125.1 min
hump + burst    1241    345    480       112 min      89 min     0.0026        12.6 min
flash crowd     1197    304    156       100 min      29 min        1.8         1.4 min

Adding the burst changes \(\sigma\) from 350 to 345 and the IQR from 484 to 480, so both rules still return widths of about 100 minutes. The roughness, in contrast, rises from \(2.5\times 10^{-6}\) to \(2.6\times 10^{-3}\), a factor of about a thousand, and the width the formula gives falls from 125 to 12.6 minutes. Without the burst the rules and the formula agree (116 and 92 against 125 minutes), which shows that the rules work when the day is smooth, as they were designed to. With a flash crowd the IQR does fall, from 484 to 156 minutes, but the Freedman–Diaconis width of 29 minutes is still far from the 1.4 minutes that the roughness asks for.

The auto rule and the cost function developed below are compared on three days with about the same number of events: a constant rate, the hump-and-burst day, and the flash crowd.

import numpy as np
from ssbins import *

rng = np.random.default_rng(3)
streams = {
    "flat":        (lambda t: np.full(np.shape(t), 0.84), 1.0),
    "hump + burst": (rate, 7.0),
    "flash crowd": (lambda t: 0.35 + 93.0 * np.exp(-0.5 * ((np.asarray(t) - 1013.0) / 3.0) ** 2), 100.0),
}
print("%-13s %6s  %-22s  %s" % ("stream", "events", "numpy 'auto' (bins, width)", "SS width"))
for name, (fn, rmax) in streams.items():
    x = simulate(rng, fn, rmax)
    nb = len(np.histogram_bin_edges(x, bins="auto", range=(0, T))) - 1
    N, _, _ = select_bins(x, grid=np.arange(1, 1441), shifts=8)
    ss = "no resolved width (one bin)" if N == 1 else "%.1f min (%d bins)" % (T / N, N)
    print("%-13s %6d  %3d bins, %5.0f min        %s" % (name, len(x), nb, T / nb, ss))

stream        events  numpy 'auto' (bins, width)  SS width
flat            1171   12 bins,   120 min        no resolved width (one bin)
hump + burst    1244   17 bins,    85 min        13.5 min (107 bins)
flash crowd     1139   58 bins,    25 min        1.6 min (889 bins)

The auto widths follow the event count and the spread: 120 minutes for the constant day, 85 for the hump and burst, 25 for the flash crowd. For the hump-and-burst day the 85-minute bins flatten a burst whose standard deviation is 6 minutes, and the cost-based width is 13.5 minutes. For the flash crowd the cost-based width is 1.6 minutes. For the constant-rate day the cost function returns a single bin, which says there is nothing to resolve; no default rule can say that.

Estimating the error without the true rate

MISE needs the true rate, and the true rate is what is being estimated. Shimazaki and Shinomoto (2007) showed that for a Poisson process the unknown part of MISE can be replaced by a quantity computed from the bin counts, up to a constant that does not depend on the width. The derivation takes expectations under the Poisson model.

Expanding the error

Expanding \((\hat{\theta} – \theta)^2 = \hat{\theta}^2 – 2\hat{\theta}\theta + \theta^2\), the last term does not involve the histogram, so it adds the same constant to MISE at every width. For the other two terms, use the Poisson facts above: \(\mathbb{E}\hat{\theta}^2 = \bin{\theta}_i^2 + \bin{\theta}_i/(n\Delta)\) on bin \(i\), and \(\int_{\text{bin } i}\theta\,dt = \Delta\bin{\theta}_i\). Summing over bins,

\[ \mathrm{MISE} = \frac{1}{T}\Big[\frac{1}{n}\sum_i \bin{\theta}_i – \Delta\sum_i \bin{\theta}_i^{\,2}\Big] + \frac{1}{T}\int_0^T \theta^2\,dt . \]

Now move to count units. Let \(\mu_i = n\Delta\bin{\theta}_i\) be the expected count in bin \(i\), \(m = \frac{1}{N}\sum_i \mu_i = n\Delta r\) the average expected count, and \(s^2 = \frac{1}{N}\sum_i (\mu_i – m)^2\) the spread of the expected counts across bins. Using \(N = T/\Delta\),

\[ \mathrm{MISE} = \frac{m – m^2 – s^2}{(n\Delta)^2} + \frac{1}{T}\int_0^T \theta^2\,dt . \]

The first term of the numerator is the variance part, \(m/(n\Delta)^2 = r/(n\Delta)\). The \(m^2\) term equals \(n^2\Delta^2 r^2\), so after the division it is the constant \(r^2\). The width-dependent structure is \(s^2\): the more the expected counts differ from bin to bin, the more of the true shape the bins resolve, and the lower the error.

What the counts can tell you

Now take the quantities the analyst can compute: the mean \(\bin{k} = \frac{1}{N}\sum_i k_i\) and the variance across bins \(v = \frac{1}{N}\sum_i (k_i – \bin{k})^2\), with divisor \(N\) rather than \(N-1\). If counts are Poisson, \(\mathbb{E}k_i^2 = \mu_i + \mu_i^2\), and \(\mathbb{E}\bin{k}^2 = m^2 + m/N\). Therefore \(\mathbb{E}v = m(1 – 1/N) + s^2\), and

\[ \mathbb{E}\big[2\bin{k} – v\big] = m + \frac{m}{N} – s^2 . \]

Compare the last two displays. The variation \(v\) across bins is the structure \(s^2\) plus Poisson noise of about \(\bin{k}\); subtracting it from \(2\bin{k}\) leaves roughly the noise minus the structure that the bins resolved, which is what MISE measures. Define the cost

\[ C(\Delta) = \frac{2\bin{k} – v}{(n\Delta)^2}. \]

Dividing the expectation of \(2\bin{k} – v\) by \((n\Delta)^2\) and subtracting MISE leaves terms that do not depend on \(\Delta\):

\[ \mathbb{E}\,C(\Delta) = \mathrm{MISE}(\Delta) – \Big[\frac{1}{T}\int_0^T\theta^2\,dt – r^2\Big] + \frac{r}{nT}. \]

The bracket is the variance of the true rate over the window, and \(r/(nT)\) is a constant. So the expected cost has the same minimizer as MISE, though its values are not MISE and can be negative. The only property of the Poisson process used is \(\mathbb{E}k_i^2 = \mu_i + \mu_i^2\), that is, variance equal to mean. The glossary entry on mean integrated squared error makes the same point about the trick: it gives an estimate of the best setting, not a guarantee of an identical choice for every sample.

The identity is an algebraic statement about expectations, so a simulation can test it. The next program draws 20,000 independent days from the synthetic rate and averages \(C\) at six widths, then compares with the right-hand side computed from the exact MISE. It repeats the check with seven pooled days. The same program also verifies a second reading of the cost, described below.

import numpy as np
from ssbins import *
from truth import *

rng = np.random.default_rng(7)
r = TH.mean()
var_t = (TH ** 2).mean() - r ** 2            # variance of the true rate over the window
print("mean rate r = %.4f /min, variance of the rate over time = %.4f" % (r, var_t))

for n_w, reps in ((1, 20000), (7, 2000)):
    sims = [np.concatenate([simulate(rng) for _ in range(n_w)]) for _ in range(reps)]
    print("\n%d pooled window(s), %d simulated datasets" % (n_w, reps))
    print("width(min)   mean cost   (s.e.)    MISE - Var + r/(nT)")
    for N in (6, 24, 60, 96, 144, 288):
        c = np.array([ss_cost(x, N, n_windows=n_w) for x in sims])
        pred = mise_exact(N, n_w) - var_t + r / (n_w * T)
        print("%9.1f %11.4f %9.4f %14.4f" % (T / N, c.mean(), c.std() / np.sqrt(reps), pred))

# leave-one-out identity for a single window: C = L/T + (K/T)^2, L = (2K - sum k^2)/delta
x = simulate(rng); K = len(x)
print("\nleave-one-out check on one dataset of K = %d events" % K)
for N in (24, 96, 288):
    k = counts(x, N); delta = T / N
    L = (2 * K - (k ** 2).sum()) / delta
    print("N=%3d   cost %.6f    LOO score/T + (K/T)^2 = %.6f" % (N, ss_cost(x, N), L / T + (K / T) ** 2))

mean rate r = 0.8654 /min, variance of the rate over time = 0.2982

1 pooled window(s), 20000 simulated datasets
width(min)   mean cost   (s.e.)    MISE - Var + r/(nT)
    240.0     -0.0997    0.0001        -0.0995
     60.0     -0.1540    0.0002        -0.1538
     24.0     -0.1872    0.0003        -0.1868
     15.0     -0.2238    0.0004        -0.2235
     10.0     -0.1937    0.0004        -0.1934
      5.0     -0.1195    0.0004        -0.1194

7 pooled window(s), 2000 simulated datasets
width(min)   mean cost   (s.e.)    MISE - Var + r/(nT)
    240.0     -0.1032    0.0001        -0.1031
     60.0     -0.1667    0.0002        -0.1667
     24.0     -0.2186    0.0003        -0.2183
     15.0     -0.2739    0.0004        -0.2735
     10.0     -0.2686    0.0004        -0.2681
      5.0     -0.2688    0.0004        -0.2682

leave-one-out check on one dataset of K = 1211 events
N= 24   cost -0.165763    LOO score/T + (K/T)^2 = -0.165763
N= 96   cost -0.273090    LOO score/T + (K/T)^2 = -0.273090
N=288   cost -0.190405    LOO score/T + (K/T)^2 = -0.190405

The averaged costs differ from the right-hand side by at most 0.0006, between one and two standard errors, in both tables, and the pooled-days table checks the \(n\Delta\) scaling. All costs are negative; among the six widths shown, the lowest average cost is at 15 minutes with one day and with seven pooled days.

A second reading: leave-one-out cross-validation

For a single window (\(n = 1\)) the cost has a second interpretation. The error \(\int\hat{\theta}^2 – 2\int\hat{\theta}\theta\) (the part of the integrated squared error that depends on the histogram) can be estimated from the events themselves. The first term is \(\sum_i k_i^2/\Delta\). The second, \(\int\hat{\theta}\,\theta\,dt\), is the expected sum of \(\hat{\theta}\) evaluated at the events, but the estimate near an event was built with that event in it, so it must be left out: for the event \(j\) in bin \(i\), the leave-one-out estimate is \((k_i – 1)/\Delta\). Summing over the \(k_i\) events in each bin gives \(\sum_i k_i(k_i – 1)/\Delta\). The leave-one-out score is

\[ L(\Delta) = \frac{1}{\Delta}\sum_i k_i^2 – \frac{2}{\Delta}\sum_i k_i(k_i-1) = \frac{2K – \sum_i k_i^2}{\Delta}, \]

where \(K\) is the number of events. With \(\bin{k} = K/N\) and \(N = T/\Delta\), the cost is \(C(\Delta) = L(\Delta)/T + (K/T)^2\), and \((K/T)^2\) does not depend on \(\Delta\). The two derivations agree: the cost is leave-one-out cross-validation of the squared error, with an event as the unit left out, a form of least-squares cross-validation for histograms. For a Poisson process, removing one event leaves the others distributed as the full process, which is what makes the leave-one-out comparison unbiased. The glossary entry on cross-validation explains the general idea of scoring a model on data it was not fitted to. The last three lines of the program’s output compare \(C\) with \(L/T + (K/T)^2\) on one dataset; they match to six decimals, as the algebra requires.

The selection procedure

The procedure is a scan over bin counts. For each candidate number of bins \(N\) on a grid, count events in \(N\) equal bins, compute \(\bin{k}\) and \(v\), evaluate \(C\), and take the \(N\) with the smallest cost. Three implementation choices matter.

• The grid starts at \(N = 1\). A single bin has \(v = 0\), so for one window its cost is \(2K/T^2\). Including it lets the scan end at “one bin” when nothing beats it, which is what happens on average for a constant rate, as shown below.
• The grid ends where the width reaches the resolution of the timestamps or where the expected count per bin becomes small. A minimum at the last grid point means the scan stopped too early, not that the finest width is best.
• The bin grid has an arbitrary origin, and a burst that straddles a bin edge looks different from one centered in a bin. The cost is therefore averaged over several rotations of the grid. Rotation by wrapping is valid when the window is a cycle, such as a day folded onto itself; for an open window, a partial bin at each end must be dropped instead.

def ss_cost(times, n_bins, n_windows=1, shifts=1, phi=1.0):
    """Cost (2*phi*kbar - v) / (n*delta)^2 averaged over `shifts` grid origins.
    phi = 1 is the Shimazaki-Shinomoto cost for a Poisson process."""
    delta = T / n_bins
    out = []
    for s in np.arange(shifts) * delta / shifts:
        k = counts(times, n_bins, s)
        out.append((2 * phi * k.mean() - k.var()) / (n_windows * delta) ** 2)  # k.var() divides by n_bins
    return float(np.mean(out))


def select_bins(times, n_windows=1, grid=np.arange(1, 361), shifts=1, phi=1.0):
    """Return (best number of bins, grid, cost at each grid point)."""
    cost = np.array([ss_cost(times, int(N), n_windows, shifts, phi) for N in grid])
    return int(grid[np.argmin(cost)]), grid, cost

np.var divides by the number of bins by default, which matches the definition of \(v\). The n_windows argument supplies the \(n\) in \((n\Delta)^2\); the phi argument is 1 for the Poisson cost and is used later. The next program prints the cost curve for the day from the first example at eleven of its widths, once with a single grid origin and once averaged over eight.

import numpy as np
from ssbins import *

rng = np.random.default_rng(8)
x = simulate(rng)                                         # the same day as in the width table
N1, g, c1 = select_bins(x, shifts=1)
N8, _, c8 = select_bins(x, shifts=8)
print("width(min)  bins   cost, one origin   cost, 8 origins")
for N in (1, 2, 6, 12, 24, 48, 72, 96, 144, 240, 360):
    print("%9.1f %5d %16.4f %16.4f" % (T / N, N, c1[N - 1], c8[N - 1]))
print("minimum: one origin %.1f min, 8 origins %.1f min" % (T / N1, T / N8))
print("gain of the minimum over the single bin: %.3f" % (c8[0] - c8.min()))

width(min)  bins   cost, one origin   cost, 8 origins
   1440.0     1           0.0012           0.0012
    720.0     2          -0.0111          -0.0500
    240.0     6          -0.1003          -0.1178
    120.0    12          -0.1215          -0.1337
     60.0    24          -0.1613          -0.1540
     30.0    48          -0.1829          -0.1767
     20.0    72          -0.2052          -0.1878
     15.0    96          -0.1945          -0.1875
     10.0   144          -0.1637          -0.1725
      6.0   240          -0.1182          -0.1220
      4.0   360          -0.0577          -0.0540
minimum: one origin 14.5 min, 8 origins 17.6 min
gain of the minimum over the single bin: 0.193

The one-bin cost is positive and every other value is negative, so only differences between rows mean anything. The cost falls as the width grows from 4 to 15 minutes and as it shrinks from 720 to 20 minutes, and is flat near the bottom: with eight origins the cost at 20 minutes is \(-0.1878\) and at 15 minutes \(-0.1875\). A single origin puts the minimum at 14.5 minutes and eight origins at 17.6. With eight origins the cost stays within about 0.015 of its minimum from 10 to 30 minutes, so the data support a range of widths better than a single value. The gain of the minimum over the single bin is 0.193, the number used in the constant-rate check below.

How well the selected width performs

Because the true rate is known in the simulation, each selected width can be scored by its integrated squared error against the truth. The next program draws 200 independent days and compares the four NumPy rules, the cost-based width with one and with eight origins, and an oracle that picks the best width for each day with hindsight. Every histogram is scored averaged over eight grid origins, so that none benefits from where the burst happens to fall relative to a bin edge. The scoring functions are in the module truth.py; they integrate the squared difference between a histogram and the true rate on a fine grid.

def ise(times, n_bins, shift=0.0):
    from ssbins import counts
    k = counts(times, n_bins, shift)
    est = k[bin_index(n_bins, shift)] / (T / n_bins)
    return float(((est - TH) ** 2).sum() * DG / T)


def ise_avg(times, n_bins, S=8):
    """ISE averaged over S origins of the bin grid, so that no result depends on where a bin edge happens to fall."""
    return float(np.mean([ise(times, n_bins, s * (T / n_bins) / S) for s in range(S)]))

import numpy as np
from ssbins import *
from truth import *

rng = np.random.default_rng(2024)
RULES = ("sturges", "scott", "fd", "auto")
OR_GRID = np.unique(np.geomspace(2, 360, 60).astype(int))      # candidate widths for the oracle

reps = 200
out = {k: [] for k in RULES + ("ss1", "ss8", "oracle")}
w = {k: [] for k in ("ss1", "ss8", "oracle")}
for _ in range(reps):
    x = simulate(rng)
    for rule in RULES:
        N = len(np.histogram_bin_edges(x, bins=rule, range=(0, T))) - 1
        out[rule].append(ise_avg(x, N))
    N1, _, _ = select_bins(x, shifts=1)
    N8, _, _ = select_bins(x, shifts=8)
    out["ss1"].append(ise_avg(x, N1)); out["ss8"].append(ise_avg(x, N8))
    w["ss1"].append(T / N1); w["ss8"].append(T / N8)
    e = np.array([ise_avg(x, int(n)) for n in OR_GRID])
    out["oracle"].append(e.min()); w["oracle"].append(T / OR_GRID[e.argmin()])
print("mean integrated squared error over %d simulated days (events/min)^2" % reps)
for k, v in out.items():
    print("  %-8s %.4f" % (k, np.mean(v)))
print("selected width in minutes, 10th / 50th / 90th percentile")
for k, v in w.items():
    print("  %-8s %s" % (k, np.percentile(v, [10, 50, 90]).round(1)))
print("best single width by exact MISE (origin-averaged):")
Ns = np.arange(2, 361); m = np.array([mise_exact(int(n), shifts=8) for n in Ns])
print("  N=%d, width %.1f min, MISE %.4f" % (Ns[m.argmin()], T / Ns[m.argmin()], m.min()))

mean integrated squared error over 200 simulated days (events/min)^2
  sturges  0.1634
  scott    0.1609
  fd       0.1544
  auto     0.1544
  ss1      0.0982
  ss8      0.0979
  oracle   0.0936
selected width in minutes, 10th / 50th / 90th percentile
  ss1      [13.4 17.6 24.4]
  ss8      [12.5 15.3 21.2]
  oracle   [13.8 15.  18. ]
best single width by exact MISE (origin-averaged):
  N=93, width 15.5 min, MISE 0.0941

The four default rules average 0.154 to 0.163, which is 1.6 to 1.7 times the 0.0941 of the best single width by exact MISE. The cost-based selection averages 0.098 with either number of origins, about 4 percent above that best fixed width, and the day-by-day oracle reaches 0.0936. The selection recovers nearly all of the gap that the defaults leave. The comparison uses one rate function and one sample size; the ratios will differ for other rates, but the mechanism, that the roughness enters the optimal width and the defaults do not measure it, does not.

The selected width varies from day to day: the middle 80 percent of selected widths runs from 12.5 to 21.2 minutes with eight origins, against 13.8 to 18.0 for the oracle. A width chosen from a single day is an estimate with its own sampling error, which the cost curve alone does not show.

How much the chosen width varies

One day gives one curve and one minimizer. The next program looks at the sampling variability of the chosen width on a single day. It splits the events into two random halves, runs the selection on each half, and reports the busiest bin. It then applies the usual bootstrap, resampling events with replacement, and compares with fresh days from the same process. Randomly sending each event to half A or half B is thinning; for a Poisson process the two halves are independent Poisson processes with half the rate.

import numpy as np
from ssbins import *

rng = np.random.default_rng(51)
x = simulate(rng)
N, _, _ = select_bins(x, shifts=8)
print("full day: %d events, selected width %.1f min" % (len(x), T / N))

# Random thinning: each event goes to half A or half B with probability 1/2.
# For a Poisson process the halves are independent Poisson processes with half the rate.
side = rng.random(len(x)) < 0.5
for name, h in (("A", x[side]), ("B", x[~side])):
    n_h, _, _ = select_bins(h, shifts=8)
    k = counts(h, n_h)
    print("half %s: %d events, selected width %.1f min, busiest bin centred at minute %.0f"
          % (name, len(h), T / n_h, (np.argmax(k) + 0.5) * T / n_h))

# Resampling events with replacement, the usual bootstrap, applied to the same day.
widths = []
for _ in range(300):
    xb = rng.choice(x, size=len(x), replace=True)
    n_b, _, _ = select_bins(xb, shifts=8)
    widths.append(T / n_b)
print("bootstrap with replacement: selected width 10/50/90%% = %s min" % np.percentile(widths, [10, 50, 90]).round(1))
print("fresh days from the same process: ", end="")
fresh = [T / select_bins(simulate(rng), shifts=8)[0] for _ in range(300)]
print("selected width 10/50/90%% = %s min" % np.percentile(fresh, [10, 50, 90]).round(1))

full day: 1250 events, selected width 11.3 min
half A: 648 events, selected width 29.4 min, busiest bin centred at minute 1014
half B: 602 events, selected width 14.3 min, busiest bin centred at minute 1019
bootstrap with replacement: selected width 10/50/90% = [4.  4.8 7.1] min
fresh days from the same process: selected width 10/50/90% = [12.5 15.3 20. ] min

On this day the full selection is 11.3 minutes, while fresh days from the same process give a middle 80 percent of 12.5 to 20.0 minutes, so this day falls below the 10th percentile. The halves select 29.4 and 14.3 minutes, a wide disagreement, yet both put the busiest bin at minute 1014 or 1019, within a few minutes of the burst at 1013. The structure replicates across halves even though the width does not, which is the evidence that the burst is more than noise, the check recommended in the glossary entry on oversmoothing and undersmoothing.

The bootstrap result runs the wrong way: resampling with replacement selects 4.0 to 7.1 minutes. A resample of \(K\) events leaves out about a third of the events (a fraction \(1/e\)) and repeats about a quarter of them at least twice. Repeated events sit at identical times, so a resample contains clusters the original does not. The cost assumes variance equal to mean, and clusters violate it in the way the section on non-Poisson events analyzes. The bootstrap is the right tool for many statistics, but this one reads the arrival pattern of individual events, and plain resampling changes that pattern. Thinning, or replicated windows when the data have them, varies the sample without creating ties.

When no width resolves anything

Take a constant rate \(\theta\). Then the bias is zero at every width, the variance of the rate estimate is \(r/(n\Delta)\), and the expected cost is \(r/(n\Delta) + r/(nT)\), which decreases as bins widen. Its minimum is at the widest option, the single bin. A procedure that scans widths will report exactly that on average. In an individual day, noise creates dips in the curve, and the program below measures how often a spurious interior minimum appears and how large it is.

import numpy as np
from ssbins import *

rng = np.random.default_rng(11)
flat = lambda t: np.full(np.shape(t), 1210.0 / T)          # constant rate, about 1210 events a day
grid = np.arange(1, 361)

def best_and_gain(x):
    N, _, c = select_bins(x, grid=grid)
    return N, c[0] - c.min()                                # gain over the single-bin cost C(N=1)

N_null, gain_null = zip(*[best_and_gain(simulate(rng, flat, 1.5)) for _ in range(1000)])
N_null, gain_null = np.array(N_null), np.array(gain_null)
print("1000 constant-rate days:")
print("  minimum at the single bin: %.1f%%   at 2-3 bins: %.1f%%   at 4 or more bins: %.1f%% (largest N = %d)"
      % (100 * (N_null == 1).mean(), 100 * ((N_null >= 2) & (N_null <= 3)).mean(),
         100 * (N_null >= 4).mean(), N_null.max()))
print("  gain over the single bin, 50th / 95th / 99th percentile:", np.percentile(gain_null, [50, 95, 99]).round(4))

print("three days from the hump-and-burst rate:")
for seed in range(3):
    N, g = best_and_gain(simulate(np.random.default_rng(seed)))
    print("  seed %d: chosen width %.1f min, gain %.3f" % (seed, T / N, g))

1000 constant-rate days:
  minimum at the single bin: 65.9%   at 2-3 bins: 17.0%   at 4 or more bins: 17.1% (largest N = 35)
  gain over the single bin, 50th / 95th / 99th percentile: [0.     0.0031 0.0051]
three days from the hump-and-burst rate:
  seed 0: chosen width 14.5 min, gain 0.273
  seed 1: chosen width 20.9 min, gain 0.219
  seed 2: chosen width 21.8 min, gain 0.230

For 1,000 days with a constant rate and about 1,210 events, the minimum is at the single bin in 65.9 percent, at two or three bins in 17.0 percent, and at four or more bins in 17.1 percent, with the largest at 35 bins. An interior minimum alone is therefore weak evidence of structure. What separates the cases is the gain of the minimum over the single-bin cost: for constant-rate days it has a 99th percentile of 0.0051, and for the three hump-and-burst days it is 0.22 to 0.27, about 40 times larger. A usable rule compares the observed gain with a constant-rate simulation: simulate days from a constant rate with the same \(K\) and \(T\), record the gain of each, and report structure only if the observed gain is far beyond that distribution. When it is not, report the data as consistent with a constant rate of \(K/T\) events per minute, with no bin width resolved, and plot a single bin or a flat line instead of a histogram that invites interpretation. The gain test only says that the data do not support a width at this sample size. A weak real feature can be present and too small to detect with that many events.

When events are not Poisson

The cost derivation used \(\mathbb{E}k_i^2 = \mu_i + \mu_i^2\). Real event streams break this when events arrive in groups. In a hypothetical service where every request is retried once a fraction of a minute later, each logical event appears twice, so the counts have variance about twice their mean. Suppose instead that \(\operatorname{Var}k_i = \varphi\,\mu_i\) with a dispersion factor \(\varphi\) common to all bins. Repeating the algebra, \(\mathbb{E}v = \varphi m(1 – 1/N) + s^2\), so the numerator of the cost has expectation \((2 – \varphi)m – s^2\) up to terms without \(\Delta\), while the true MISE has \(\varphi m – s^2\). The noise penalty in the cost is \((2-\varphi)\,r/(n\Delta)\), where the truth is \(\varphi\,r/(n\Delta)\). At \(\varphi = 1\) the two agree. At \(\varphi = 2\) the penalty vanishes, so nothing discourages narrow bins, and above 2 it rewards them.

The next program generates the retry stream with the same expected rate as the original day, measures how far the counts depart from Poisson, and runs the selection on 200 days of each stream. The helper with_retries and the dispersion measure divide each bin’s squared deviation by its true expected count, so a Poisson stream scores about 1.

def with_retries(rng, mean_delay=0.5):
    """Every logical event is delivered twice: the original and a retry a short random time later.
    Originals arrive at half the rate, so the expected event rate is still rate(t)."""
    p = simulate(rng, lambda t: rate(t) / 2, 3.5)
    retry = (p + rng.exponential(mean_delay, len(p))) % T
    return np.sort(np.concatenate([p, retry]))

import numpy as np
from ssbins import *
from truth import *

rng = np.random.default_rng(21)
grid = np.arange(2, 721)

def dispersion(x, N):
    """Mean of (k_i - mu_i)^2 / mu_i over bins, with mu_i the true expected count: 1 for a Poisson process."""
    k = counts(x, N)
    mu = np.bincount(bin_index(N), weights=TH, minlength=N) * DG
    return ((k - mu) ** 2 / mu).mean()

print("%-14s %10s  %-29s %s" % ("stream", "dispersion", "selected width (10/50/90%)", "mean ISE"))
for name, gen in (("Poisson", lambda: simulate(rng)), ("with retries", lambda: with_retries(rng))):
    w, d, e = [], [], []
    for _ in range(200):
        x = gen()
        N, _, _ = select_bins(x, grid=grid)
        w.append(T / N); d.append(dispersion(x, 96)); e.append(ise_avg(x, N))
    print("%-14s %10.2f  %-29s %.3f" % (name, np.mean(d), np.percentile(w, [10, 50, 90]).round(1), np.mean(e)))

stream         dispersion  selected width (10/50/90%)    mean ISE
Poisson              1.01  [13.2 18.1 24.4]              0.098
with retries         1.97  [ 3.5  6.7 13.2]              0.298

The retry stream has a dispersion of 1.97 and the cost selects bins with a median width of 6.7 minutes instead of 18.1; its error is 0.298 against 0.098, three times larger, because the narrow bins chase the retried pairs. The glossary entry on the Poisson process notes that clustered arrivals are overdispersed and that methods built on variance-equals-mean then compute from a premise the data contradict; this is such a case, with a measured cost.

Detecting it in one window

The dispersion above used the true expected counts, which an analyst does not have. A statistic that needs only the counts compares adjacent bins. If two neighboring bins have nearly the same expected count and events are Poisson, then given their total \(a + b\), the count \(a\) is binomial with probability one half, and \((a-b)^2/(a+b)\) has mean 1. A slow change of the rate across the window does not matter, because the two bins are close in time. Clustering raises the average above 1. The function below averages the statistic over pairs of adjacent bins.

def neighbour_dispersion(x, width):
    """Mean of (a - b)^2 / (a + b) over adjacent bin pairs (a, b) with a + b > 0.
    If the rate is nearly equal in the two bins and events are Poisson, each term has mean 1;
    slow changes of the rate between pairs do not matter, clustering pushes the value above 1."""
    k = counts(x, int(T / width))
    a, b = k[0::2], k[1::2]
    keep = (a + b) > 0
    return ((a - b)[keep] ** 2 / (a + b)[keep]).mean()

import numpy as np
from ssbins import *

rng = np.random.default_rng(5)
print("pair width   Poisson   with retries   (mean over 200 days)")
for w in (2.5, 5, 10, 20):
    q_p = np.mean([neighbour_dispersion(simulate(rng), w) for _ in range(200)])
    q_r = np.mean([neighbour_dispersion(with_retries(rng), w) for _ in range(200)])
    print("%9.1f min %8.2f %12.2f" % (w, q_p, q_r))

pair width   Poisson   with retries   (mean over 200 days)
      2.5 min     1.01         1.68
      5.0 min     1.03         1.90
     10.0 min     1.10         2.03
     20.0 min     1.83         2.67

At pair widths of 2.5 to 10 minutes the Poisson stream gives 1.01 to 1.10, and the retry stream gives 1.7 to 2.0. At 20 minutes the Poisson value rises to 1.83, because pairs near the burst are no longer alike, so the statistic needs bins narrow enough for the rate to be locally constant, and a value that grows with width signals that the rate is changing, not that events cluster. Comparing counts across intervals with comparable expected rates is what the glossary entry on the Poisson process recommends, and this is one way to do it.

The estimate can be fed back. With \(\hat{\varphi}\) from the 5-minute pairs (floored at 1), the cost numerator becomes \(2\hat{\varphi}\bin{k} – v\), whose expectation is \(\varphi m – s^2\) up to terms without \(\Delta\), as MISE requires; this is the phi argument of ss_cost. The program compares it with the plain cost.

import numpy as np
from truth import *
from ssbins import *

rng = np.random.default_rng(61)
grid = np.arange(2, 721)
print("%-14s %-22s %-8s %-28s %s" % ("stream", "cost", "mean phi", "selected width (10/50/90%)", "mean ISE"))
for name, gen in (("Poisson", lambda: simulate(rng)), ("with retries", lambda: with_retries(rng))):
    for label in ("phi = 1", "phi from neighbours"):
        w, e, ph = [], [], []
        for _ in range(200):
            x = gen()
            phi = 1.0 if label == "phi = 1" else max(1.0, neighbour_dispersion(x, 5.0))
            N, _, _ = select_bins(x, grid=grid, phi=phi)
            w.append(T / N); e.append(ise_avg(x, N)); ph.append(phi)
        print("%-14s %-22s %-8.2f %-28s %.3f" % (name, label, np.mean(ph), np.percentile(w, [10, 50, 90]).round(1), np.mean(e)))

stream         cost                   mean phi selected width (10/50/90%)   mean ISE
Poisson        phi = 1                1.00     [13.1 18.2 24.4]             0.100
Poisson        phi from neighbours    1.07     [13.9 18.2 24.4]             0.098
with retries   phi = 1                1.00     [ 3.3  6.5 12.6]             0.306
with retries   phi from neighbours    1.89     [14.8 20.  34.3]             0.147

On the retry stream the corrected cost selects 20.0 minutes (median) with an error of 0.147 instead of 0.306, and on the Poisson stream it makes no difference (0.098 against 0.100). The error stays above the Poisson 0.098 because the retried events carry less information than independent ones. The correction assumes one \(\varphi\) for all bins; a stream whose clustering varies in time needs a different treatment. If the retries are not the object of interest, the better remedy is to change what an event is, by deduplicating retries or counting sessions, so that the unit being counted is closer to independent.

Excess variation across pooled windows is a separate phenomenon. The next program pools seven days three ways: independent Poisson days; days that share the daily shape but have a day-level multiplier drawn from a gamma distribution with mean 1; and days with retries. It reports the variance divided by the mean across days for each 60-minute bin, averaged over bins, and the selected width.

import numpy as np
from ssbins import *

rng = np.random.default_rng(31)
n_days = 7
grid = np.arange(2, 361)

def poisson_day(mult=1.0):
    return simulate(rng, lambda t: mult * rate(t), 7.0 * max(mult, 1.0) + 1.0)

def across_day_dispersion(days, N):
    """Variance over days divided by the mean over days, for each of N bins, averaged over bins (1 for Poisson)."""
    K = np.array([counts(d, N) for d in days], dtype=float)         # shape (days, bins)
    m = K.mean(axis=0)
    return (K[:, m > 0].var(axis=0, ddof=1) / m[m > 0]).mean()

scenarios = {
    "independent Poisson days": lambda: [poisson_day() for _ in range(n_days)],
    "day-level multiplier":     lambda: [poisson_day(m) for m in rng.gamma(6.0, 1 / 6.0, n_days)],
    "retry duplicates":         lambda: [with_retries(rng) for _ in range(n_days)],
}
print("%-26s %-18s %s" % ("%d pooled days" % n_days, "across-day Var/mean", "selected width in min (10/50/90%)"))
for name, gen in scenarios.items():
    d, w = [], []
    for _ in range(100):
        days = gen()
        N, _, _ = select_bins(np.concatenate(days), n_windows=n_days, grid=grid)
        d.append(across_day_dispersion(days, 24)); w.append(T / N)
    print("%-26s %-18.2f %s" % (name, np.mean(d), np.percentile(w, [10, 50, 90]).round(1)))

7 pooled days              across-day Var/mean selected width in min (10/50/90%)
independent Poisson days   1.00               [ 6.3  9.9 16.7]
day-level multiplier       9.06               [ 6.4 10.3 16.7]
retry duplicates           1.95               [4.1 4.9 6.6]

The multiplier days have a variance-to-mean ratio of 9.06 across days, far from the Poisson value 1.00, yet they select the same width as independent days: median 10.3 minutes against 9.9. The retry days have a ratio of only 1.95 and select 4.9 minutes. Given the day multipliers, the pooled counts are exactly Poisson with the rate scaled by their sum, so the cost is still working on a Poisson stream whose shape is the daily shape. Retries add noise inside bins, which is what the cost misreads. Overdispersion that comes from a shared level damages something else: pooled counts are not Poisson even though each day is, so error bins that assume Poisson noise are too narrow. The excess variation is a signal to look for clustering, dependence, or a mis-specified rate model; pooling more windows does not by itself establish the model.

Practical choices for event timestamps

Several decisions sit outside the cost function.

• Window and origin. The code wraps the grid around the window, which suits a day folded onto itself. For an open window such as a deployment week, use only whole bins inside the window and do not wrap. Midnight alignment is a convention that makes plots readable; check that a conclusion survives a shifted origin.
• Timestamp resolution. Logs with one-second resolution or coarser contain ties. Many events with an identical timestamp are clusters, as in the retry example, and bins narrower than the resolution count only the rounding. Set the grid’s finest width above the resolution and examine batch-written events separately.
• Which windows to pool. Pooling \(n\) windows shrinks \(\Delta^{*}\) as \(n^{-1/3}\) only if the windows share a rate shape. Mondays can be pooled with Mondays; a holiday, an incident day, or a seasonal change breaks the assumption. The glossary entry on seasonality separates calendar-fixed patterns, which pooling can exploit, from irregular ones, which it cannot.
• Reporting. State the width and the unit. Plot events per minute rather than raw counts, so that plots at different widths remain comparable. Show the chosen width together with a wider and a narrower one, because the cost curve is flat near its minimum.
• One width for everything. A fixed width compromises between flat stretches, where wide bins would be better, and bursts, where narrow bins are needed. When the rate has both, variable-width methods can serve each region separately; Bayesian blocks of Scargle and colleagues segment event data into blocks of constant rate, and Shimazaki and Shinomoto extended the same cost-minimization idea to variable bin widths and to kernel smoothers.

The checks fit together as a routine for a new event stream:

1. Scan the cost over a grid of widths that includes a single bin, averaged over several origins.
2. Compare the gain of the minimum over the single bin with a constant-rate simulation of the same size, and stop with a flat rate if the gain is not far outside it.
3. Test neighboring bins for clustering, and rescale the cost or redefine the event if the statistic is well above 1.
4. Confirm that the features at the chosen width appear in random halves of the data.

References

• Shimazaki, H. and Shinomoto, S. “A method for selecting the bin size of a time histogram.” Neural Computation 19(6), 1503–1527 (2007), with the author’s page that states the cost function and the Poisson assumption.
• Shimazaki, H. and Shinomoto, S. “Kernel bandwidth optimization in spike rate estimation.” Journal of Computational Neuroscience 29, 171–182 (2010), the kernel version of the same method, with the author’s toolbox page.
• Scott, D. W. “On optimal and data-based histograms.” Biometrika 66, 605–610 (1979), the source of the \(n^{-1/3}\) width formula and the Gaussian reference rule.
• NumPy documentation, numpy.histogram_bin_edges, formulas and stated limits of the Sturges, Scott, Freedman–Diaconis and auto rules.
• Scargle, J. D., Norris, J. P., Jackson, B. and Chiang, J. “Studies in Astronomical Time Series Analysis. VI. Bayesian Block Representations.” The Astrophysical Journal 764, 167 (2013), variable-width segmentation of event data.
• Astropy documentation, astropy.stats.bayesian_blocks, an implementation of Bayesian blocks with an events fitness.