Algoritmi per il calcolo della varianza

Gli algoritmi per il calcolo della varianza giocano un ruolo molto importante nella statistica computazionale. Una difficoltà chiave nel progetto di un buon algoritmo per questo problema è che le formule per la varianza possono includere somme di quadrati, che possono condurre a instabilità numerica così come overflow aritmetico quando vengono trattati grandi valori.

Una formula per calcolare la varianza di un'intera popolazione di dimensione N è:

${\displaystyle \sigma ^{2}={\bar {(x^{2})}}-{\bar {x}}^{2}=\displaystyle {\frac {\sum _{i=1}^{N}x_{i}^{2}-(\sum _{i=1}^{N}x_{i})^{2}/N}{N}}.\!}$

Usando la correzione di Bessel per calcolare una stima priva di bias della varianza della popolazione da un campione finito di n osservazioni, la formula è:

${\displaystyle s^{2}=\displaystyle {\frac {\sum _{i=1}^{n}x_{i}^{2}-(\sum _{i=1}^{n}x_{i})^{2}/n}{n-1}}.\!}$

Pertanto, un algoritmo banale per calcolare la stima della varianza è dato da:

def stimaVarianza(x):    n = length(x)    somma = 0    sommaQuad = 0    for xi in x:       somma = somma + xi       sommaQuad = sommaQuad + xi^2    return (sommaQuad - somma^2 / n) / (n - 1) 

Questo algoritmo può facilmente essere adattato al calcolo della varianza di una popolazione finita: basta dividere per N invece di n − 1 nell'ultima linea.

Poiché ${\displaystyle SommaQuad}$ e ${\displaystyle {\frac {Somma^{2}}{n}}}$ possono essere numeri molto simili, la cancellazione può far sì che la precisione del risultato sia molto inferiore alla precisione intrinseca dell'aritmetica in virgola mobile usata per effettuare il calcolo. Così, questo algoritmo non deve essere usato nella pratica.^[1][2] Questo è particolarmente svantaggioso se la deviazione standard è piccola rispetto alla media. Tuttavia, l'algoritmo può essere migliorato adottando il metodo della media ipotizzata.

Possiamo usare una proprietà della varianza per evitare la cancellazione numerica in questa formula, poiché la varianza è invariante rispetto ai cambiamenti di changes in un parametro di locazione

${\displaystyle \operatorname {Var} (X-k)=\operatorname {Var} (X).}$

con ${\displaystyle k}$ costante qualunque, il che conduce alla nuova formula

${\displaystyle s^{2}=\displaystyle {\frac {\sum _{i=1}^{n}(x_{i}-K)^{2}-\left[\sum _{i=1}^{n}(x_{i}-K)\right]^{2}/n}{n-1}}.\!}$

e più vicino è ${\displaystyle K}$ al valore medio, più accurato il risultato sarà, ma soltanto lo scegliere un valore compreso nell'intervallo dei campioni garantirà la stabilitata desiderata. Se i valori ${\displaystyle (x_{i}-K)}$ sono piccoli, allora non ci sono problemi con la somma dei loro quadrati; al contrario, se essi sono grandi, ciò significa necessariamente che la varianza è anch'essa grande. In ogni caso, il secondo termine nella formula è sempre più piccolo del primo, quindi nessuna cancellazione può avvenire.^[2]

Se prendiamo solo il primo campione come ${\displaystyle K}$, l'algoritmo può essere scritto nel linguaggio di programmazione Python come segue

def shifted_data_variance(data):    if len(data) == 0:       return 0    K = data[0]    n = 0    sum_ = 0    sum_sqr = 0    for x in data:       n = n + 1       sum_ += x - K       sum_sqr += (x - K) * (x - K)    variance = (sum_sqr - (sum_ * sum_)/n)/(n - 1)    # use n instead of (n-1) if want to compute the exact variance of the given data    # use (n-1) if data are samples of a larger population    return variance 

questa formula facilitata così il calcolo incrementale, che può essere espresso come segue

K = 0 n = 0 ex = 0 ex2 = 0   def add_variable(x):     if (n == 0):       K = x     n = n + 1     ex += x - K     ex2 += (x - K) * (x - K)   def remove_variable(x):     n = n - 1     ex -= (x - K)     ex2 -= (x - K) * (x - K)   def get_meanvalue():    return K + ex / n  def get_variance():     return (ex2 - (ex*ex)/n) / (n-1) 

Un approccio alternativo, usando una formula diversa per la varianza, consiste prima nel calcolare la media campione,

${\displaystyle {\bar {x}}=\displaystyle {\frac {\sum _{j=1}^{n}x_{j}}{n}}}$,

e poi calcolare la somma dei quadrati delle differenze dalla media,

${\displaystyle \mathrm {varianza} =s^{2}=\displaystyle {\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})^{2}}{n-1}}\!}$,

dove s è la deviazione standard. Ciò è dato dal seguente codice sorgente:

def two_pass_variance(data):     n = 0     sum1 = 0     sum2 = 0          for x in data:         n += 1         sum1 += x          mean = sum1 / n      for x in data:         sum2 += (x - mean)*(x - mean)          variance = sum2 / (n - 1)     return variance 

Questo algoritmo è numericamente stabile se n è piccolo.^[1][3] Tuttavia, i risultati di entrambi questi algoritmi semplici ("banale" e "a due passi") possono dipendere non commutativamente dall'ordinamento dei dati e possono dare scarsi risultati per insiemi grandi di dati a causa di errori di arrotondamento ripetuti che si accumulano nelle somme. Tecniche quali la somma compensata possono essere usate per ridurre questo errore.

La versione con somma compensata dell'algoritmo visto sopra è:^[4]

def compensated_variance(data):     n = 0     sum1 = 0     for x in data:         n += 1         sum1 += x     mean = sum1/n           sum2 = 0     sum3 = 0     for x in data:         sum2 += (x - mean)**2         sum3 += (x - mean)     variance = (sum2 - sum3**2/n)/(n - 1)     return variance 

Spesso è utile essere in grado di calcolare la varianza in un passo singolo, elaborando ogni valore ${\displaystyle x_{i}}$ una sola volta; per esempio, quando i dati vengono raccolti senza abbastanza spazio per conservare tutti i valori, o quando i costi dell'accesso alla memoria sono preponderanti tra quelli del calcolo. Per un tale algoritmo in linea, è richiesta una relazione di ricorrenza tra quantità dalle quali le statistiche richieste possono essere calcolate in un modo numericamente stabile.

Le seguenti formule possono essere usate per aggiornare la media e la varianza (stimata) della sequenza, per un elemento addizionale ${\displaystyle x_{\mathrm {new} }}$. Qui, ${\displaystyle {\bar {x}}_{n}}$ denota la media campione dei primi ${\displaystyle n}$ campioni ${\displaystyle (x_{1},\,\ldots ,\,x_{n})}$, ${\displaystyle s_{n}^{2}}$ la loro varianza campione e ${\displaystyle \sigma _{n}^{2}}$ la varianza della loro popolazione.

${\displaystyle {\bar {x}}_{n}={\frac {(n-1)\,{\bar {x}}_{n-1}+x_{n}}{n}}={\bar {x}}_{n-1}+{\frac {x_{n}-{\bar {x}}_{n-1}}{n}}\!}$

${\displaystyle s_{n}^{2}={\frac {(n-2)}{(n-1)}}\,s_{n-1}^{2}+{\frac {(x_{n}-{\bar {x}}_{n-1})^{2}}{n}},\quad n>1}$

${\displaystyle \sigma _{n}^{2}={\frac {(n-1)\,\sigma _{n-1}^{2}+(x_{n}-{\bar {x}}_{n-1})(x_{n}-{\bar {x}}_{n})}{n}}.}$

Queste formule soffrono di instabilità numerica. Una quantità migliore per l'aggiornamento è la somma dei quadrati delle differenze dalla media corrente ${\displaystyle \textstyle \sum _{i=1}^{n}(x_{i}-{\bar {x}}_{n})^{2}}$, qui indicate come ${\displaystyle M_{2,n}}$:

${\displaystyle M_{2,n}\!=M_{2,n-1}+(x_{n}-{\bar {x}}_{n-1})(x_{n}-{\bar {x}}_{n})}$

${\displaystyle s_{n}^{2}={\frac {M_{2,n}}{n-1}}}$

${\displaystyle \sigma _{n}^{2}={\frac {M_{2,n}}{n}}}$

Un algoritmo numericamente stabile per la varianza campione è indicato sotto. Esso calcola anche la media. Questo algoritmo è dovuto a Knuth,^[5] che cita Welford,^[6] ed è stato accuratamente analizzato.^[7][8] È anche comune denotare ${\displaystyle M_{k}={\bar {x}}_{k}}$ e ${\displaystyle S_{k}=M_{2,k}}$.^[9]

def online_variance(data):     n = 0     mean = 0.0     M2 = 0.0           for x in data:         n += 1         delta = x - mean         mean += delta/n         M2 += delta*(x - mean)      if n < 2:         return float('nan')     else:         return M2 / (n - 1) 

Questo algoritmo è molto meno incline a perdita di precisione a causa della cancellazione numerica, ma potrebbe non essere efficiente a causa dell'operazione di divisione all'interno del ciclo. Per un algoritmo a due passi particolarmente robusto per calcolare la varianza, si può prima calcolare e sottrarre una stima della media e poi utilizzare questo algoritmo sui residui.

Sotto, l'algoritmo parallelo illustra come unire più insiemi di statistiche calcolate in linea.

L'algoritmo può essere esteso per gestire pesi campione disuguali, sostituendo il semplice contatore ${\displaystyle n}$ con la somma dei pesi vista finora. West (1979)^[10] suggerisce questo algoritmo incrementale:

def weighted_incremental_variance(dataWeightPairs):     sumweight = 0     mean = 0     M2 = 0      for x, weight in dataWeightPairs:  # Alternatively "for x, weight in zip(data, weights):"         temp = weight + sumweight         delta = x - mean         R = delta * weight / temp         mean += R         M2 += sumweight * delta * R  # Alternatively, "M2 = M2 + weight * delta * (x−mean)"         sumweight = temp      variance_n = M2/sumweight     variance = variance_n * len(dataWeightPairs)/(len(dataWeightPairs) - 1) 

Chan e altri.^[4] Notare che l'algoritmo in linea indicato sopra è un caso particolare di algoritmo che funziona per ogni partizione del campione ${\displaystyle X}$ negli insiemi ${\displaystyle X_{A}}$, ${\displaystyle X_{B}}$:

${\displaystyle \delta \!={\bar {x}}_{B}-{\bar {x}}_{A}}$

${\displaystyle {\bar {x}}_{X}={\bar {x}}_{A}+\delta \cdot {\frac {n_{B}}{n_{X}}}}$

${\displaystyle M_{2,X}=M_{2,A}+M_{2,B}+\delta ^{2}\cdot {\frac {n_{A}n_{B}}{n_{X}}}}$.

Ciò può essere utile quando, per esempio, più unità di elaborazione possono essere assegnate a parti discrete di input.

Il metodo di Chan per stimare la media è instabile numericamente quando ${\displaystyle n_{A}\approx n_{B}}$ ed entrambi sono grandi, poiché l'errore numerico in ${\displaystyle {\bar {x}}_{B}-{\bar {x}}_{A}}$ non è ridotto in modo da essere nel caso ${\displaystyle n_{B}=1}$. In tali casi, è preferibile che ${\displaystyle {\bar {x}}_{X}={\frac {n_{A}{\bar {x}}_{A}+n_{B}{\bar {x}}_{B}}{n_{A}+n_{B}}}}$.

Supponiamo che tutte le operazioni in virgola mobile usino l'aritmetica con lo standard IEEE 754 con precisione doppia. Consideriamo il campione ${\displaystyle (4,\,7,\,13,\,16)}$ da una popolazione infinita. Sulla base di questo campione, la media stimata della popolazione è 10 e la stima priva di bias della varianza della popolazione è 30. Entrambi gli algoritmi, quello "banale" e quello "a due passi", calcolano questi valori correttamente. Ora, invece, consideriamo il campione ${\displaystyle (10^{8}+4,\,10^{8}+7,\,10^{8}+13,\,10^{8}+16)}$, che dà luogo alla stessa varianza stimata come per il primo campione. L'algoritmo "a due passi" calcola correttamente questa varianza stimata, ma l'algoritmo "banale" restituisce ${\displaystyle 29.33333333333333}$ invece di 30. Mentre questa perdita di precisione può essere tollerabile e considerata come un difetto minore dell'algoritmo "a due passi", è facile trovare dati che rivelano un grave difetto nell'algoritmo "banale": Prendiamo il campione ${\displaystyle (10^{9}+4,\,10^{9}+7,\,10^{9}+13,\,10^{9}+16)}$. Nuovamente, la varianza stimata della popolazione, di valore 30, è calcolata correttamente dall'algoritmo "a due passi", ma l'algoritmo "banale" ora calcola ${\displaystyle -170.66666666666666}$. Ciò costituisce un serio problema con l'algoritmo "banale" ed è dovuto alla cancellazione numerica nella sottrazione di due numeri simili nella fase finale dell'algoritmo.

Terriberry^[11] estende le formule di Chan al calcolo del terzo e quarto momento centrale, necessarie per esempio quando si stimano asimmetria e curtosi:

${\displaystyle M_{3,X}=M_{3,A}+M_{3,B}+\delta ^{3}{\frac {n_{A}n_{B}(n_{A}-n_{B})}{n_{X}^{2}}}+3\delta {\frac {n_{A}M_{2,B}-n_{B}M_{2,A}}{n_{X}}}}$

${\displaystyle {\begin{aligned}M_{4,X}=M_{4,A}+M_{4,B}&+\delta ^{4}{\frac {n_{A}n_{B}\left(n_{A}^{2}-n_{A}n_{B}+n_{B}^{2}\right)}{n_{X}^{3}}}\\&+6\delta ^{2}{\frac {n_{A}^{2}M_{2,B}+n_{B}^{2}M_{2,A}}{n_{X}^{2}}}+4\delta {\frac {n_{A}M_{3,B}-n_{B}M_{3,A}}{n_{X}}}\\\end{aligned}}}$

Qui, gli ${\displaystyle M_{k}}$ sono nuovamente le somme di potenze di differenze dalla la media ${\displaystyle \Sigma (x-{\overline {x}})^{k}}$, avendo

asimmetria: ${\displaystyle g_{1}={\frac {{\sqrt {n}}M_{3}}{M_{2}^{3/2}}},}$

cutrosi: ${\displaystyle g_{2}={\frac {nM_{4}}{M_{2}^{2}}}-3.}$

Per il caso incrementale (cioè, ${\displaystyle B=\{x\}}$), ciò si riduce a:

${\displaystyle \delta \!=x-m}$

${\displaystyle m'=m+{\frac {\delta }{n}}}$

${\displaystyle M_{2}'=M_{2}+\delta ^{2}{\frac {n-1}{n}}}$

${\displaystyle M_{3}'=M_{3}+\delta ^{3}{\frac {(n-1)(n-2)}{n^{2}}}-{\frac {3\delta M_{2}}{n}}}$

${\displaystyle M_{4}'=M_{4}+{\frac {\delta ^{4}(n-1)(n^{2}-3n+3)}{n^{3}}}+{\frac {6\delta ^{2}M_{2}}{n^{2}}}-{\frac {4\delta M_{3}}{n}}}$

Conservando il valore ${\displaystyle \delta /n}$, è necessaria una sola operazione di divisione e le statistiche di ordine superiore possono così essere calcolate per poco costo incrementale.

Un esempio dell'algoritmo in linea per la curtosi implementato come descritto è:

def online_kurtosis(data):     n = 0     mean = 0     M2 = 0     M3 = 0     M4 = 0      for x in data:         n1 = n         n = n + 1         delta = x - mean         delta_n = delta / n         delta_n2 = delta_n * delta_n         term1 = delta * delta_n * n1         mean = mean + delta_n         M4 = M4 + term1 * delta_n2 * (n*n - 3*n + 3) + 6 * delta_n2 * M2 - 4 * delta_n * M3         M3 = M3 + term1 * delta_n * (n - 2) - 3 * delta_n * M2         M2 = M2 + term1      kurtosis = (n*M4) / (M2*M2) - 3     return kurtosis 

Pébay^[12] estende ulteriormente questi risultati per momenti centrali di ordine arbitrario, per il caso incrementale e per il caso di coppie. Si possono trovare formule simili per la covarianza.

Choi e Sweetman ^[13] offrono due metodi alternativi per calcolare l'asimmetria e la curtosi, ognuno dei quali può far risparmiare requisiti di memoria del computer e tempo di occupazione della CPU in certe applicazioni. Il primo approccio consiste nel calcolare i momenti statistici separando i dati in pacchetti e poi calcolare i momenti dalla geometria dell'istogramma risultante, che effettivamente diventa un algoritmo a un passo per momenti superiori. Un vantaggio è che il calcolo dei momenti statistici può essere effettuato con accuratezza arbitraria tale che i calcoli possono essere portati, per esempio, alla precisione del formato di memorizzazione dei dati o alla capacità originale dell'hardware. Un istogramma relativo ad una variabile casuale può essere costruito nel modo convenzionale: l'intervallo di valori potenziali è diviso in pacchetti e il numero di occorrenze all'interno di ogni pacchetto è contato e graficato in modo tale che l'area di ogni rettangolo sia uguale alla porzione di valori campione in questo pacchetto:

${\displaystyle H(x_{k})={\frac {h(x_{k})}{A}}}$

dove ${\displaystyle h(x_{k})}$ e ${\displaystyle H(x_{k})}$ rappresentano la frequenza e la frequenza relativa per il pacchetto ${\displaystyle x_{k}}$ e ${\displaystyle A=\sum _{k=1}^{K}h(x_{k})\,\Delta x_{k}}$ è l'area totale dell'istogramma. Dopo questa normalizzazione, i momenti di origine zero (raw moments) ${\displaystyle n}$ e i momenti centrali di ${\displaystyle x(t)}$ possono essere calcolati dal relativo istogramma:

${\displaystyle m_{n}^{(h)}=\sum _{k=1}^{K}x_{k}^{n}\,H(x_{k})\Delta x_{k}={\frac {1}{A}}\sum _{k=1}^{K}x_{k}^{n}\,h(x_{k})\Delta x_{k}}$

${\displaystyle \theta _{n}^{(h)}=\sum _{k=1}^{K}{\Big (}x_{k}-m_{1}^{(h)}{\Big )}^{n}\,H(x_{k})\Delta x_{k}={\frac {1}{A}}\sum _{k=1}^{K}{\Big (}x_{k}-m_{1}^{(h)}{\Big )}^{n}\,h(x_{k})\Delta x_{k}}$

dove la sovrascritta ${\displaystyle ^{(h)}}$ indica che i momenti sono calcolati dall'istogramma. Per un pacchetto costante di larghezza ${\displaystyle \Delta x_{k}=\Delta x}$ queste due espressioni possono essere semplificate usando ${\displaystyle I=A/\Delta x}$:

${\displaystyle m_{n}^{(h)}={\frac {1}{I}}{\sum _{k=1}^{K}x_{k}^{n}\,h(x_{k})}}$

${\displaystyle \theta _{n}^{(h)}={\frac {1}{I}}{\sum _{k=1}^{K}{\Big (}x_{k}-m_{1}^{(h)}{\Big )}^{n}\,h(x_{k})}}$

Il secondo approccio, da Choi and Sweetman ^[13], consiste in una metodologia analitica per combinare momenti statistici da segmenti individuali di una storia temporale tale che i momenti complessivi risultanti siano quelli di una storia temporale completa. Questa metodologia potrebbe essere utilizzata per il calcolo parallelo di momenti statistici con successiva combinazione di quei momenti, o per la combinazione di momenti statistici calcolati in tempi sequenziali.

Se gli insiemi ${\displaystyle Q}$ di momenti statistici sono noti: ${\displaystyle (\gamma _{0,q},\mu _{q},\sigma _{q}^{2},\alpha _{3,q},\alpha _{4,q})\quad }$ per ${\displaystyle q=1,2,...,Q}$, allora ogni ${\displaystyle \gamma _{n}}$ può essere espresso in termini dei momenti di origine zero (raw moments) equivalenti ${\displaystyle n}$ :

${\displaystyle \gamma _{n,q}=m_{n,q}\gamma _{0,q}\qquad \quad {\textrm {per}}\quad n=1,2,3,4\quad {\text{ e }}\quad q=1,2,\dots ,Q}$

dove, generalmente, come ${\displaystyle \gamma _{0,q}}$ si prende la durata della storia temporale ${\displaystyle q^{th}}$, o il numero di punti se ${\displaystyle \Delta t}$ è costante.

Il vantaggio di esprimere i momenti statistici in termini di ${\displaystyle \gamma }$ è che gli insiemi ${\displaystyle Q}$ possono essere combinati per addizione e non c'è un limite superiore per il valore di ${\displaystyle Q}$.

${\displaystyle \gamma _{n,c}=\sum _{q=1}^{Q}\gamma _{n,q}\quad \quad {\textrm {for}}\quad n=0,1,2,3,4}$

dove la sovrascritta ${\displaystyle _{c}}$ rappresenta la storia temporale ${\displaystyle \gamma }$ concatenata o combinata. Questi valori combinati di ${\displaystyle \gamma }$ possono essere trasformati inversamente in momenti di origine zero (raw moments) rappresentando la storia temporale concatenata completa

${\displaystyle m_{n,c}={\frac {\gamma _{n,c}}{\gamma _{0,c}}}\quad {\textrm {for}}\quad n=1,2,3,4}$

relazioni note tra i momenti di origine zero (raw moments) (${\displaystyle m_{n}}$) e il momento centrale ${\displaystyle \left(\theta _{n}=E[(x-\mu )^{n}]\right)}$ sono allora usate per calcolare i momenti centrali della storia temporale concatenata. Infine, i momenti statistici della storia concatenata vengono calcolati dai momenti centrali:

${\displaystyle \mu _{c}=m_{1,c}\ \ \ \ \ \sigma _{c}^{2}=\theta _{2,c}\ \ \ \ \ \alpha _{3,c}={\frac {\theta _{3,c}}{\sigma _{c}^{3}}}\ \ \ \ \ \alpha _{4,c}={\frac {\theta _{4,c}}{\sigma _{c}^{4}}}-3}$

Si possono utilizzare algoritmi molto simili per calcolare la covarianza. L'algoritmo "banale" è:

${\displaystyle \operatorname {Cov} (X,Y)=\displaystyle {\frac {\sum _{i=1}^{n}x_{i}y_{i}-(\sum _{i=1}^{n}x_{i})(\sum _{i=1}^{n}y_{i})/n}{n}}.\!}$

Per l'algoritmo scritto sopra, si può utilizzare il seguente codice Python:

def naive_covariance(data1, data2):     n = len(data1)     sum12 = 0     sum1 = sum(data1)     sum2 = sum(data2)      for i in range(n):         sum12 += data1[i]*data2[i]      covariance = (sum12 - sum1*sum2 / n) / n     return covariance 

Come per la variaznza, la covarianza di due variabili casuali è anche invariante per traslazione, allora, dati due valori costanti ${\displaystyle K_{x}}$ e ${\displaystyle K_{y}}$ qualsiasi, si può scrivere:

${\displaystyle \operatorname {Cov} (X,Y)=\operatorname {Cov} (X-k_{x},Y-k_{y})=\displaystyle {\frac {\sum _{i=1}^{n}(x_{i}-K_{x})(y_{i}-K_{y})-(\sum _{i=1}^{n}(x_{i}-K_{x}))(\sum _{i=1}^{n}(y_{i}-K_{y}))/n}{n}}.\!}$

e, nuovamente, scegliendo un valore all'interno dell'intervallo dei valori stabilizzerà la formula contro la cancellazione numerica in modo da renderlo più robusto nei confronti delle grandi somme. Prendendo il primo valore in ogni insieme di dati, l'algoritmo può essere scritto come segue:

def shifted_data_covariance(dataX, dataY):    n = len(dataX)    if (n < 2):      return 0    Kx = dataX[0]    Ky = dataY[0]    Ex = 0    Ey = 0    Exy = 0    for i in range(n):       Ex += dataX[i] - Kx       Ey += dataY[i] - Ky       Exy += (dataX[i] - Kx) * (dataY[i] - Ky)    return (Exy - Ex * Ey / n) / n 

L'algoritmo a due passi prima calcola le medie campione e dopo la covarianza:

${\displaystyle {\bar {x}}=\displaystyle \sum _{i=1}^{n}x_{i}/n}$

${\displaystyle {\bar {y}}=\displaystyle \sum _{i=1}^{n}y_{i}/n}$

${\displaystyle \operatorname {Cov} (X,Y)=\displaystyle {\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}{n}}.\!}$

L'algoritmo a due passi può essere scritto come segue:

def two_pass_covariance(data1, data2):     n = len(data1)      mean1 = sum(data1) / n     mean2 = sum(data2) / n      covariance = 0      for i in range(n):         a = data1[i] - mean1         b = data2[i] - mean2         covariance += a*b / n      return covariance 

Una versione compensata un po' più accurata esegue l'algoritmo "banale" completo sui residui. Le somme finali ${\displaystyle \textstyle \sum x_{i}}$ e ${\displaystyle \textstyle \sum y_{i}}$ dovrebbero essere zero, ma il secondo passaggio compensa qualsiasi piccolo errore.

Una lieve modifica dell'algoritmo in linea per il calcolo della varianza produce un algoritmo in linea per la covarianza:

def online_covariance(data1, data2):     mean1 = mean2 = 0.     M12 = 0.     n = len(data1)     for i in range(n):         delta1 = (data1[i] - mean1) / (i + 1)         mean1 += delta1         delta2 = (data2[i] - mean2) / (i + 1)         mean2 += delta2         M12 += i * delta1 * delta2 - M12 / (i + 1)     return n / (n - 1.) * M12 

Esiste un algoritmo stabile a un passo, simile a quello visto sopra, che calcola il co-momento ${\displaystyle \textstyle C_{n}=\sum _{i=1}^{n}(x_{i}-{\bar {x}}_{n})(y_{i}-{\bar {y}}_{n})}$:

${\displaystyle {\bar {x}}_{n}={\bar {x}}_{n-1}+{\frac {x_{n}-{\bar {x}}_{n-1}}{n}}\!}$

${\displaystyle {\bar {y}}_{n}={\bar {y}}_{n-1}+{\frac {y_{n}-{\bar {y}}_{n-1}}{n}}\!}$

${\displaystyle C_{n}=C_{n-1}+(x_{n}-{\bar {x}}_{n})(y_{n}-{\bar {y}}_{n-1})=C_{n-1}+(y_{n}-{\bar {y}}_{n})(x_{n}-{\bar {x}}_{n-1})}$

L'asimmetria apparente in quest'ultima equazione è dovuta al fatto che ${\displaystyle \textstyle (x_{n}-{\bar {x}}_{n})={\frac {n-1}{n}}(x_{n}-{\bar {x}}_{n-1})}$, così entrambi i termini di aggiornamento sono pari a ${\displaystyle \textstyle {\frac {n-1}{n}}(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n-1})}$. Un'accuratezza ancora maggiore può essere ottenuta calcolando prima le medie e dopo usando l'algoritmo stabile a un passo sui residui.

Così, possiamo calcolare la covarianza come

${\displaystyle {\begin{aligned}\operatorname {Cov} _{N}(X,Y)={\frac {C_{N}}{N}}&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+(x_{n}-{\bar {x}}_{n})(y_{n}-{\bar {y}}_{n-1})}{N}}\\&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+(y_{n}-{\bar {y}}_{n})(x_{n}-{\bar {x}}_{n-1})}{N}}\\&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+{\frac {N-1}{N}}(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n-1})}{N}}.\end{aligned}}}$

Allo stesso modo, vi è una formula per combinare le covarianze di due insiemi che può essere utilizzata per parallelizzare calcolo:

${\displaystyle C_{X}=C_{A}+C_{B}+({\bar {x}}_{A}-{\bar {x}}_{B})({\bar {y}}_{A}-{\bar {y}}_{B})\cdot {\frac {n_{A}n_{B}}{n_{X}}}}$.

• Formula computazionale per la varianza
• Algoritmo di sommatoria di Kahan

• (EN) Eric W. Weisstein, Algoritmi per il calcolo della varianza, su MathWorld, Wolfram Research.

Portale Statistica: accedi alle voci di Wikipedia che trattano di statistica