Streaming Estimation¶

streaming ¶

Tail index estimation over a stream, without holding the sample.

The Hill estimator needs the top k order statistics. That sounds like it rules out a streaming version, since order statistics are a global property of the sample -- but only the top k + 1 values ever matter, and which values those are can be maintained incrementally. A min-heap of size k + 1 keeps them in O(k) memory and O(log k) time per observation, whatever the length of the stream.

That gives an exact result, not an approximation. :class:StreamingTailIndex holds the same numbers a batch estimator would sort out of the full sample, so it returns the same estimate to the last bit -- which the test suite asserts rather than assumes. Nothing is traded away except the ability to ask a different k later.

Two semantics, and the difference matters:

:class:StreamingTailIndex The whole stream, in O(k) memory. Every observation ever seen contributes, so a change in the tail is diluted by everything before it.

:class:WindowedTailIndex The most recent window observations, in O(window) memory. Responds to a change in the tail, at a cost in memory and variance.

The memory difference is not an implementation shortcoming to be fixed later. When the largest value in a window expires, the new maximum can be any of the remaining ones, so a window cannot be summarised more compactly than by keeping it. Anything claiming to track windowed order statistics in O(k) is either approximating or wrong.

StreamingTailIndex ¶

StreamingTailIndex(k)

Tail index estimators over an unbounded stream, in O(k) memory.

Holds the top k + 1 observations and nothing else. The estimate is exact: it is the same number a batch estimator would produce from the whole sample, because the estimators depend on the sample only through those values.

Parameters:

Name	Type	Description	Default
`k`	`int`	Number of top order statistics the estimators use.	required

Raises:

Type	Description
`ValueError`	If `k` is less than two, which no estimator here accepts.

Examples:

>>> from heavytails import Pareto
>>> from heavytails.tail_index import hill_estimator
>>> data = Pareto(alpha=2.0, xm=1.0).rvs(20000, seed=1)
>>> stream = StreamingTailIndex(k=500)
>>> stream.extend(data)
>>> stream.hill() == hill_estimator(data, k=500)
True

Source code in heavytails/streaming.py

def __init__(self, k: int) -> None:
    if not isinstance(k, int) or k < 2:
        raise ValueError("k must be an integer of at least 2.")
    self._k = k
    self._top = TopK(k + 1)

k `property` ¶

Number of top order statistics used.

n_seen `property` ¶

n_seen

How many observations have passed through.

ready `property` ¶

ready

Whether enough observations have been seen to estimate anything.

threshold `property` ¶

threshold

The order statistic the estimators measure exceedances above.

extend ¶

extend(values)

Add many observations.

Source code in heavytails/streaming.py

def extend(self, values: Iterable[float]) -> None:
    """Add many observations."""
    self._top.extend(values)

hill ¶

hill()

The Hill estimate of the extreme-value index.

Returns:

Type	Description
`float`	`gamma`, equal to `1 / alpha` for a Pareto tail.

Raises:

Type	Description
`ValueError`	If fewer than `k + 1` observations have been seen, or if the retained values are not positive.

Source code in heavytails/streaming.py

def hill(self) -> float:
    """
    The Hill estimate of the extreme-value index.

    Returns:
        ``gamma``, equal to ``1 / alpha`` for a Pareto tail.

    Raises:
        ValueError: If fewer than ``k + 1`` observations have been seen, or
            if the retained values are not positive.
    """
    return _hill(self._ready(), self._k)

moment ¶

moment()

The Dekkers-Einmahl-de Haan moment estimate.

Unlike the Hill estimator this does not assume the tail is heavy, so it is the one to reach for when the sign of gamma is in doubt.

Returns:

Type	Description
`tuple[float, float]`	`(gamma, alpha)`.

Raises:

Type	Description
`ValueError`	If there is not yet enough data, or it is not positive.

Source code in heavytails/streaming.py

def moment(self) -> tuple[float, float]:
    """
    The Dekkers-Einmahl-de Haan moment estimate.

    Unlike the Hill estimator this does not assume the tail is heavy, so it
    is the one to reach for when the sign of ``gamma`` is in doubt.

    Returns:
        ``(gamma, alpha)``.

    Raises:
        ValueError: If there is not yet enough data, or it is not positive.
    """
    return _moment(self._ready(), self._k)

update ¶

update(value)

Add one observation.

Source code in heavytails/streaming.py

def update(self, value: float) -> None:
    """Add one observation."""
    self._top.push(value)

TopK ¶

TopK(k)

The k largest values of a stream, in O(k) memory.

A min-heap holds the retained values, so its root is the smallest of them and is what a new observation has to beat. That is the whole trick: the comparison that decides whether to keep a value is against the smallest kept, not against anything in the discarded majority.

Parameters:

Name	Type	Description	Default
`k`	`int`	How many values to retain, at least one.	required

Raises:

Type	Description
`ValueError`	If `k` is not a positive integer.

Examples:

>>> top = TopK(3)
>>> top.extend([5.0, 1.0, 9.0, 3.0, 7.0])
>>> top.descending()
[9.0, 7.0, 5.0]
>>> top.n_seen, len(top)
(5, 3)

Source code in heavytails/streaming.py

def __init__(self, k: int) -> None:
    if not isinstance(k, int) or k < 1:
        raise ValueError("k must be a positive integer.")
    self._k = k
    self._heap: list[float] = []
    self._seen = 0

capacity `property` ¶

capacity

How many values are retained once the stream is long enough.

n_seen `property` ¶

n_seen

How many observations have been offered.

len ¶

__len__()

How many values are currently retained.

Source code in heavytails/streaming.py

def __len__(self) -> int:
    """How many values are currently retained."""
    return len(self._heap)

descending ¶

descending()

The retained values, largest first.

Source code in heavytails/streaming.py

def descending(self) -> list[float]:
    """The retained values, largest first."""
    return sorted(self._heap, reverse=True)

extend ¶

extend(values)

Offer many observations.

Source code in heavytails/streaming.py

def extend(self, values: Iterable[float]) -> None:
    """Offer many observations."""
    for value in values:
        self.push(value)

push ¶

push(value)

Offer one observation.

Source code in heavytails/streaming.py

def push(self, value: float) -> None:
    """Offer one observation."""
    self._seen += 1
    item = float(value)
    if len(self._heap) < self._k:
        heapq.heappush(self._heap, item)
    elif item > self._heap[0]:
        heapq.heapreplace(self._heap, item)

WindowedTailIndex ¶

WindowedTailIndex(window, k)

The same estimators over the most recent window observations.

The whole-stream version dilutes a change in the tail with everything that came before, which is the wrong behaviour for monitoring: a portfolio whose tail index has moved from 3 to 1.5 does not want an estimate averaging the two. This forgets.

Memory is O(window), not O(k), and that is inherent rather than a gap to close later. When the largest value in the window expires, the new largest can be any of the survivors, so nothing smaller than the window itself determines the answer.

Parameters:

Name	Type	Description	Default
`window`	`int`	How many recent observations to keep.	required
`k`	`int`	Number of top order statistics the estimators use, below `window`.	required

Raises:

Type	Description
`ValueError`	If `window` is not a positive integer, or `k` is not at least two and strictly below `window`.

Examples:

>>> from heavytails import Pareto
>>> monitor = WindowedTailIndex(window=5000, k=400)
>>> monitor.extend(Pareto(alpha=3.0, xm=1.0).rvs(5000, seed=1))
>>> round(monitor.hill(), 3)  # true gamma 0.333
0.33
>>> monitor.extend(Pareto(alpha=1.5, xm=1.0).rvs(5000, seed=2))
>>> round(monitor.hill(), 3)  # true gamma 0.667, the old regime gone
0.667

Source code in heavytails/streaming.py

def __init__(self, window: int, k: int) -> None:
    if not isinstance(window, int) or window < 2:
        raise ValueError("window must be an integer of at least 2.")
    if not isinstance(k, int) or k < 2:
        raise ValueError("k must be an integer of at least 2.")
    if k >= window:
        raise ValueError(f"k must be below window; got k={k}, window={window}")
    self._window = window
    self._k = k
    self._recent: deque[float] = deque()
    # The same values kept in ascending order, so the top k + 1 are a slice
    # rather than a sort on every query. Insertion and removal are linear
    # in the window, which for the sizes this is meant for is a memmove.
    self._ordered: list[float] = []
    self._seen = 0

k `property` ¶

Number of top order statistics used.

n_seen `property` ¶

n_seen

How many observations have passed through, including evicted ones.

ready `property` ¶

ready

Whether the window holds enough observations to estimate anything.

threshold `property` ¶

threshold

The order statistic the estimators measure exceedances above.

window `property` ¶

window

How many recent observations are kept.

extend ¶

extend(values)

Add many observations.

Source code in heavytails/streaming.py

def extend(self, values: Iterable[float]) -> None:
    """Add many observations."""
    for value in values:
        self.update(value)

hill ¶

hill()

The Hill estimate over the current window.

Source code in heavytails/streaming.py

def hill(self) -> float:
    """The Hill estimate over the current window."""
    return _hill(self._ready(), self._k)

moment ¶

moment()

The moment estimate over the current window.

Source code in heavytails/streaming.py

def moment(self) -> tuple[float, float]:
    """The moment estimate over the current window."""
    return _moment(self._ready(), self._k)

update ¶

update(value)

Add one observation, evicting the oldest if the window is full.

Source code in heavytails/streaming.py

def update(self, value: float) -> None:
    """Add one observation, evicting the oldest if the window is full."""
    item = float(value)
    self._seen += 1
    self._recent.append(item)
    bisect.insort(self._ordered, item)
    if len(self._recent) > self._window:
        oldest = self._recent.popleft()
        index = bisect.bisect_left(self._ordered, oldest)
        del self._ordered[index]

values ¶

values()

The current window contents, oldest first.

Source code in heavytails/streaming.py

def values(self) -> list[float]:
    """The current window contents, oldest first."""
    return list(self._recent)

Streaming Estimation¶

streaming ¶

StreamingTailIndex ¶

k property ¶

n_seen property ¶

ready property ¶

threshold property ¶

extend ¶

hill ¶

moment ¶

update ¶

TopK ¶

capacity property ¶

n_seen property ¶

__len__ ¶

descending ¶

extend ¶

push ¶

WindowedTailIndex ¶

k property ¶

n_seen property ¶

ready property ¶

threshold property ¶

window property ¶

extend ¶

hill ¶

moment ¶

update ¶

values ¶

k `property` ¶

n_seen `property` ¶

ready `property` ¶

threshold `property` ¶

capacity `property` ¶

n_seen `property` ¶

len ¶

k `property` ¶

n_seen `property` ¶

ready `property` ¶

threshold `property` ¶

window `property` ¶