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.

Algoritmo banale

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

$\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 è:

$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é $SommaQuad$ e ${\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. Questo è particolarmente svantaggioso se la deviazione standard è piccola rispetto alla media. Tuttavia, l'algoritmo può essere migliorato adottando il metodo della media ipotizzata.

Calcolo con dati traslati

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

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

con $k$ costante qualunque, il che conduce alla nuova formula

$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 è $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 $(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.

Se prendiamo solo il primo campione come $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)

Algoritmo a due passi

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

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

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

$\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. 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.

Variante compensata

La versione con somma compensata dell'algoritmo visto sopra è:

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

Algoritmo in linea

Spesso è utile essere in grado di calcolare la varianza in un passo singolo, elaborando ogni valore $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 $x_{\mathrm {new} }$. Qui, ${\bar {x}}_{n}$ denota la media campione dei primi $n$ campioni $(x_{1},\,\ldots ,\,x_{n})$, $s_{n}^{2}$ la loro varianza campione e $\sigma _{n}^{2}$ la varianza della loro popolazione.

${\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}}\!$

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

$\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 $\textstyle \sum _{i=1}^{n}(x_{i}-{\bar {x}}_{n})^{2}$, qui indicate come $M_{2,n}$:

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

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

$\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, che cita Welford, ed è stato accuratamente analizzato. È anche comune denotare $M_{k}={\bar {x}}_{k}$ e $S_{k}=M_{2,k}$.

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.

Algoritmo incrementale ponderato

L'algoritmo può essere esteso per gestire pesi campione disuguali, sostituendo il semplice contatore $n$ con la somma dei pesi vista finora. West (1979) 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)

Algoritmo parallelo

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

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

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

$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 $n_{A}\approx n_{B}$ ed entrambi sono grandi, poiché l'errore numerico in ${\bar {x}}_{B}-{\bar {x}}_{A}$ non è ridotto in modo da essere nel caso $n_{B}=1$. In tali casi, è preferibile che ${\bar {x}}_{X}={\frac {n_{A}{\bar {x}}_{A}+n_{B}{\bar {x}}_{B}}{n_{A}+n_{B}}}$.

Esempio

Supponiamo che tutte le operazioni in virgola mobile usino l'aritmetica con lo standard IEEE 754 con precisione doppia. Consideriamo il campione $(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 $(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 $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 $(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 $-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.

Statistiche di ordine superiore

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

$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}}}$

${\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 $M_{k}$ sono nuovamente le somme di potenze di differenze dalla la media $\Sigma (x-{\overline {x}})^{k}$, avendo

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

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

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

$\delta \!=x-m$

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

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

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

$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 $\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 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 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:

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

dove $h(x_{k})$ e $H(x_{k})$ rappresentano la frequenza e la frequenza relativa per il pacchetto $x_{k}$ e $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) $n$ e i momenti centrali di $x(t)$ possono essere calcolati dal relativo istogramma:

$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}$

$\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 $^{(h)}$ indica che i momenti sono calcolati dall'istogramma. Per un pacchetto costante di larghezza $\Delta x_{k}=\Delta x$ queste due espressioni possono essere semplificate usando $I=A/\Delta x$:

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

$\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 , 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 $Q$ di momenti statistici sono noti: $(\gamma _{0,q},\mu _{q},\sigma _{q}^{2},\alpha _{3,q},\alpha _{4,q})\quad$ per $q=1,2,...,Q$, allora ogni $\gamma _{n}$ può essere espresso in termini dei momenti di origine zero (raw moments) equivalenti $n$ :

$\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 $\gamma _{0,q}$ si prende la durata della storia temporale $q^{th}$, o il numero di punti se $\Delta t$ è costante.

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

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

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

$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) ($m_{n}$) e il momento centrale $\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:

$\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$

Covarianza

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

$\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 $K_{x}$ e $K_{y}$ qualsiasi, si può scrivere:

$\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:

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

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

$\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 $\textstyle \sum x_{i}$ e $\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 $\textstyle C_{n}=\sum _{i=1}^{n}(x_{i}-{\bar {x}}_{n})(y_{i}-{\bar {y}}_{n})$:

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

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

$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 $\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 $\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

${\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:

$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}}}$.

Voci correlate

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

Collegamenti esterni

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