Streaming Statistics and Online Learning

"In a world drowning in data, the ability to compute statistics in a single pass with bounded memory is not a luxury — it is a necessity." — Cormode & Muthukrishnan, 2005

Why Streaming Matters

Traditional batch algorithms require the entire dataset to be stored in memory. In modern applications:

• Network traffic analysis: Billions of packets per day; storing all flow records is infeasible
• Sensor networks: Continuous streams from thousands of IoT devices with limited battery and memory
• Financial tick data: Real-time computation of portfolio statistics as trades execute
• Social media analytics: Rolling computation of engagement metrics across millions of posts
• Database query processing: Approximate aggregate queries over massive tables

The streaming model forces us to rethink fundamental statistical operations: how do we compute a mean, a variance, or a quantile without storing the data?

Welford's Online Algorithm

Streaming Variance (Welford's Method)

Parallel Streaming: Chan's Algorithm

Count-Min Sketch

Query Estimation

Space Complexity

HyperLogLog

Algorithm

Bloom Filters

Online Gradient Descent

Multi-Armed Bandits

Upper Confidence Bound (UCB1)

Thompson Sampling

Python Implementation

import numpy as np
from collections import defaultdict
import hashlib

# --- Welford's Online Algorithm ---
class WelfordOnlineStatistics:
    """Track mean and variance using Welford's algorithm in one pass."""
    def __init__(self):
        self.n = 0
        self.mean = 0.0
        self.M2 = 0.0

    def update(self, x):
        self.n += 1
        delta = x - self.mean
        self.mean += delta / self.n
        delta2 = x - self.mean
        self.M2 += delta * delta2

    @property
    def variance(self):
        return self.M2 / (self.n - 1) if self.n > 1 else 0.0

    @property
    def std(self):
        return np.sqrt(self.variance)

    def __repr__(self):
        return f"Welford(n={self.n}, mean={self.mean:.4f}, var={self.variance:.4f})"

# Demonstration
np.random.seed(42)
data = np.random.normal(loc=5.0, scale=2.0, size=10000)

w = WelfordOnlineStatistics()
for x in data:
    w.update(x)

print(f"Streaming: {w}")
print(f"Batch mean: {data.mean():.4f}, Batch var: {data.var():.4f}")
# Output: Streaming mean/var matches batch exactly

# --- Count-Min Sketch ---
class CountMinSketch:
    def __init__(self, width, depth):
        self.width = width
        self.depth = depth
        self.table = np.zeros((depth, width), dtype=int)
        self.hashes = [self._make_hash(i) for i in range(depth)]

    def _make_hash(self, seed):
        def h(x):
            return int(hashlib.md5(f"{seed}:{x}".encode()).hexdigest(), 16) % self.width
        return h

    def add(self, item, count=1):
        for i in range(self.depth):
            self.table[i][self.hashes[i](item)] += count

    def estimate(self, item):
        return min(self.table[i][self.hashes[i](item)] for i in range(self.depth))

# Demonstration
cms = CountMinSketch(width=1000, depth=5)
items = np.random.choice(1000, size=100000, replace=True)
for item in items:
    cms.add(str(item))

print(f"\nCount-Min Sketch estimates for first 5 items:")
for i in range(5):
    print(f"  Item '{i}': true={np.sum(items==i)}, est={cms.estimate(str(i))}")

# --- HyperLogLog (simplified) ---
class HyperLogLog:
    def __init__(self, p=14):
        self.p = p
        self.m = 1 << p
        self.buckets = [0] * self.m

    def _hash(self, item):
        h = int(hashlib.sha256(str(item).encode()).hexdigest(), 16)
        return h

    def _leading_zeros(self, h, max_bits=64):
        if h == 0:
            return max_bits
        count = 0
        while (h & 1) == 0:
            count += 1
            h >>= 1
        return count

    def add(self, item):
        h = self._hash(item)
        bucket = h & (self.m - 1)
        w = h >> self.p
        lz = self._leading_zeros(w)
        self.buckets[bucket] = max(self.buckets[bucket], lz)

    def cardinality(self):
        alpha = 0.7213 / (1 + 1.079 / self.m)
        raw = alpha * self.m ** 2 / sum(2.0 ** (-b) for b in self.buckets)
        return raw

# Demonstration
hll = HyperLogLog(p=14)
true_n = 500000
for i in range(true_n):
    hll.add(i)

print(f"\nHyperLogLog:")
print(f"  True cardinality: {true_n}")
print(f"  Estimated: {hll.cardinality():.0f}")
print(f"  Error: {abs(hll.cardinality() - true_n) / true_n * 100:.2f}%")

# --- Bloom Filter ---
class BloomFilter:
    def __init__(self, size, num_hashes):
        self.size = size
        self.num_hashes = num_hashes
        self.bit_array = [False] * size
        self.hashes = [self._make_hash(i) for i in range(num_hashes)]

    def _make_hash(self, seed):
        def h(item):
            return int(hashlib.md5(f"{seed}:{item}".encode()).hexdigest(), 16) % self.size
        return h

    def add(self, item):
        for h in self.hashes:
            self.bit_array[h(item)] = True

    def __contains__(self, item):
        return all(self.bit_array[h(item)] for h in self.hashes)

# Demonstration
bf = BloomFilter(size=10000, num_hashes=7)
for i in range(1000):
    bf.add(f"item_{i}")

false_positives = sum(f"new_{i}" in bf for i in range(10000))
print(f"\nBloom Filter: {false_positives}/10000 false positives "
      f"(expected ~1%)")

# --- Multi-Armed Bandit (UCB1) ---
class UCB1:
    def __init__(self, k):
        self.k = k
        self.counts = np.zeros(k)
        self.values = np.zeros(k)
        self.total = 0

    def select_arm(self):
        # Pull each arm once first
        for i in range(self.k):
            if self.counts[i] == 0:
                return i
        ucb_values = self.values + np.sqrt(2 * np.log(self.total) / self.counts)
        return np.argmax(ucb_values)

    def update(self, arm, reward):
        self.counts[arm] += 1
        n = self.counts[arm]
        self.values[arm] += (reward - self.values[arm]) / n
        self.total += 1

# 10-armed Gaussian bandit test
np.random.seed(42)
true_means = np.random.normal(0, 1, 10)
optimal_arm = np.argmax(true_means)

agent = UCB1(10)
cumulative_regret = []

for t in range(1, 2001):
    arm = agent.select_arm()
    reward = np.random.normal(true_means[arm], 1.0)
    agent.update(arm, reward)
    regret = true_means[optimal_arm] - true_means[arm]
    cumulative_regret.append(regret)

print(f"\nUCB1 Bandit: final avg reward = {agent.values[optimal_arm]:.3f} "
      f"(true optimal = {true_means[optimal_arm]:.3f})")
print(f"Cumulative regret after 2000 rounds: {sum(cumulative_regret):.1f}")

Key Takeaways

Next Steps

• -> Big Data Statistics
• -> Statistics Meets Machine Learning
• -> Monte Carlo Simulation