Esercizio - Campionamento di riserva su un flusso di cinque elementi

Teoria: Algoritmi di streaming - campionamento di riserva e conteggio degli elementi distintiModello a flusso (un passaggio, memoria molto minore di n); campionamento di riserva (reservoir sampling) con dimostrazione per induzione; contatore di Morris (stima non distorta 2^X-1, varianza n(n-1)/2); elementi distinti: idea di Flajolet-Martin e stima del k-esimo minimo (k-1)/v_k; Count-Min sketch con garanzia (epsilon, delta) e dimensionamento w=e/epsilon, d=ln(1/delta).Algoritmi di streaming - campionamento di riserva e conteggio degli elementi distinti → (reservoir sampling), Linearità del valore atteso e variabili indicatriciLa linearità $E[X+Y]=E[X]+E[Y]$ vale sempre, anche senza indipendenza; con le variabili indicatrici ($E[I_A]=P(A)$) trasforma un conteggio in una somma di probabilità. Applicazioni: attesa del primo successo $1/p$, collezionista di figurine $nH_n$, punti fissi di una permutazione (media $1$), record (media $H_n$), inversioni ($n(n-1)/4$), attesa condizionata, identità di Wald $E[S_N]=E[N]E[X]$.Linearità del valore atteso e variabili indicatrici →.

Testo. Il campionamento di riserva con $k=2$ legge un flusso $x_1,\dots,x_5$: tiene $x_1,x_2$; per $i=3,4,5$ estrae $j$ uniforme in $\{1,\dots,i\}$ e, se $j\le2$, mette $x_i$ nello slot $j$.

1. Con $k=1$ e flusso $x_1,\dots,x_4$ mostrare che ogni elemento è il campione finale con probabilità $\frac14$ e simulare a mano la sequenza di estrazioni $j_2=2$, $j_3=1$, $j_4=3$.
2. Con $k=2$ calcolare, passo per passo, la probabilità con cui ogni elemento è nel campione dopo $i=3,4,5$ elementi, e la probabilità di ciascuna delle $\binom52=10$ coppie finali.
3. Un programmatore usa invece la regola: l'elemento $i>k$ entra con probabilità $\frac12$ e sostituisce uno slot scelto a caso. Calcolare le probabilità di inclusione per $k=2$, $N=5$. È uniforme?
4. Quante sostituzioni fa l'algoritmo in media su un flusso di $N$ elementi? Calcolare per $N=10^6$, $k=100$.

Soluzione

1. Caso $k=1$. L'elemento $x_i$ ($i\ge2$) sostituisce il campione se $j=1$, con probabilità $\frac1i$. $x_1$ sopravvive se i passi $2,3,4$ non lo sostituiscono: $\frac12\cdot\frac23\cdot\frac34=\frac14$. $x_2$ entra con $\frac12$ e sopravvive ai passi $3,4$: $\frac12\cdot\frac23\cdot\frac34=\frac14$. $x_3$: $\frac13\cdot\frac34=\frac14$. $x_4$: $\frac14$. Tutti $\frac14$ (in generale $\frac1j\prod_{t>j}\frac{t-1}t=\frac1N$).

Sequenza $j_2=2$, $j_3=1$, $j_4=3$. Passo $i=2$: $j=2>1$, si scarta $x_2$ e il campione resta $x_1$. Passo $3$: $j=1\le1$, $x_3$ sostituisce: campione $x_3$. Passo $4$: $j=3>1$, scarta $x_4$. Campione finale: $x_3$.

2. Caso $k=2$, $N=5$. Dopo $i=2$ il campione è $\{x_1,x_2\}$.

• Passo $3$ ($j\in\{1,2,3\}$ ciascuno con probabilità $\frac13$): $j=3$ nulla, campione $\{1,2\}$; $j=1$: $\{3,2\}$; $j=2$: $\{1,3\}$. Ogni elemento compare in due casi su tre: $P=\frac23=\frac ki$ con $k=2$, $i=3$.
• Passo $4$ ($j\in\{1,..,4\}$): $j=3,4$ lasciano il campione, $j=1,2$ inseriscono $x_4$. Le $6$ coppie possibili hanno ciascuna probabilità $\frac16$: $\{1,2\}$ (da $\{1,2\}$ con $j\ge3$: $\frac13\cdot\frac24=\frac16$), $\{1,3\},\{2,3\}$ (da $\{1,3\}$ e $\{3,2\}$ con $j\ge3$: $\frac13\cdot\frac12=\frac16$ ciascuna), e le tre coppie con $x_4$: $\{4,2\}$, $\{1,4\}$ e $\{3,4\}$, ciascuna con probabilità $\frac16$. Inclusione: $\frac36=\frac12=\frac24$ per ogni elemento.
• Passo $5$ ($j\in\{1,\dots,5\}$): $x_5$ entra con probabilità $\frac25$; ogni vecchio elemento nel campione ($\frac12$) è espulso con probabilità $\frac15$ (il suo slot è esattamente uno dei $5$ valori di $j$), quindi resta con $\frac12\cdot\frac45=\frac25$. Le $10$ coppie finali hanno probabilità $\boxed{\frac1{10}}$ ciascuna: l'enumerazione di tutti gli esiti (le $13$ configurazioni ordinate per slot si ripartiscono sulle $10$ coppie) lo conferma.

3. L'algoritmo "entra con $\frac12$". Con $k=2$ i nuovi elementi entrano con probabilità $\frac12$, nello slot $1$ o $2$ con probabilità $\frac14$ ciascuno. Un elemento già presente resta con probabilità $\frac12+\frac14=\frac34$ a ogni passo (non entra nulla, oppure entra nell'altro slot). Il calcolo esatto dà $$P(x_1\text{ presente})=P(x_2)=\frac{27}{64}=\bigl(\tfrac34\bigr)^3,\quad P(x_3)=\frac9{32},\quad P(x_4)=\frac38,\quad P(x_5)=\frac12.$$ Non è uniforme: gli ultimi elementi sono favoriti ($\frac12$ contro $\frac{27}{64}=0{,}42$; per un flusso lungo i primi sparirebbero come $(\frac34)^N$). Il fattore $\frac ki$ che decresce con $i$ è essenziale: il $j$ uniforme in $\{1..i\}$ rende sempre meno probabile l'ingresso di un elemento nuovo, compensando il fatto che gli elementi vecchi hanno già avuto più occasioni di essere espulsi.

4. Numero di sostituzioni. Il passo $i>k$ sostituisce con probabilità $\frac ki$, quindi per la linearità $$E[\#\text{sostituzioni}]=\sum_{i=k+1}^N\frac ki=k(H_N-H_k).$$ Per $N=10^6$, $k=100$: $100\,(14{,}3927-5{,}1874)=\boxed{920{,}5}$ sostituzioni su $10^6$ elementi: costo di memoria $k=100$ e circa $920$ scritture. Per $N=10^5$ la formula dà $690{,}3$ e $100$ simulazioni danno media $687{,}4$ (la deviazione standard teorica è circa $25$).

Codice di verifica

from fractions import Fraction
import random

def passo(stati, i, k):
    nuovo = {}
    for s, p in stati.items():
        for j in range(1, i + 1):
            t = s
            if j <= k:
                t = list(s); t[j - 1] = i; t = tuple(t)       # x_i prende lo slot j
            nuovo[t] = nuovo.get(t, 0) + p * Fraction(1, i)
    return nuovo

stati = {(1, 2): Fraction(1)}
for i in (3, 4, 5):
    stati = passo(stati, i, 2)
    print(i, {x: str(sum(p for s, p in stati.items() if x in s)) for x in range(1, i + 1)})
coppie = {}
for s, p in stati.items():
    coppie[frozenset(s)] = coppie.get(frozenset(s), 0) + p
print(len(stati), sorted(set(map(str, coppie.values()))), len(coppie))

random.seed(1)
def sostituzioni(N, k):
    return sum(random.randrange(1, i + 1) <= k for i in range(k + 1, N + 1))
r = [sostituzioni(10 ** 5, 100) for _ in range(100)]
print(sum(r) / len(r), 100 * (sum(1 / i for i in range(101, 10 ** 5 + 1))))

Risultato: le probabilità di inclusione 2/3, 1/2, 2/5 (per ogni elemento), 13 ['1/10'] 10, e infine 687.43 690.2768... (media simulata e valore teorico delle sostituzioni).