Visualizza wikitesto di 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 [[Stabilità numerica|instabilità numerica]] così come [[overflow aritmetico]] quando vengono trattati grandi valori.

== Algoritmo banale ==
Una formula per calcolare la varianza di un'intera [[Popolazione statistica|popolazione]] di dimensione ''N'' è:

:<math>\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}. \!</math>

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

:<math>s^2 = \displaystyle\frac {\sum_{i=1}^n x_i^2 - (\sum_{i=1}^n x_i)^2/n}{n-1}. \!</math>

Pertanto, un algoritmo banale per calcolare la stima della varianza è dato da:
<syntaxhighlight lang="python">
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)
</syntaxhighlight>

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

Poiché <math>SommaQuad</math> e <math>\frac{Somma^2}{n}</math> possono essere numeri molto simili, la [[Cancellazione numerica|cancellazione]] può far sì che la [[Cifra significativa|precisione]] del risultato sia molto inferiore alla precisione intrinseca dell'aritmetica in [[Numero in virgola mobile|virgola mobile]] usata per effettuare il calcolo.  Così, questo algoritmo non deve essere usato nella pratica.<ref name="Einarsson2005">{{Cita libro|autore=Bo Einarsson|titolo=Accuracy and Reliability in Scientific Computing|url=http://books.google.com/books?id=8hrDV5EbrEsC|accesso=17 febbraio 2013|data=1º agosto 2005|editore=SIAM|isbn=978-0-89871-584-2|p=47}}</ref><ref name="Chan1983">{{Cita pubblicazione|url=http://cpsc.yale.edu/sites/default/files/files/tr222.pdf|autore=T.F.Chan, G.H. Golub and R.J. LeVeque|titolo="Algorithms for computing the sample variance: Analysis and recommendations", The American Statistician, 37|pp=242-247|anno=1983}}</ref> 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 è [[Invarianza (matematica)|invariante]] rispetto ai cambiamenti di changes in un [[parametro di locazione]]

:<math>\operatorname{Var}(X-k)=\operatorname{Var}(X).</math>

con <math>k</math> costante qualunque, il che conduce alla nuova formula

:<math>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}. \!</math>

e più vicino è <math>K</math> 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 <math>(x_i - K)</math> 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.<ref name="Chan1983" />

Se prendiamo solo il primo campione come <math>K</math>, l'algoritmo può essere scritto nel [[Python|linguaggio di programmazione Python]] come segue

<syntaxhighlight lang="python">
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
</syntaxhighlight>
questa formula facilitata così il calcolo incrementale, che può essere espresso come segue 
<syntaxhighlight lang="python">
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)
</syntaxhighlight>

== Algoritmo a due passi ==
Un approccio alternativo, usando una formula diversa per la varianza, consiste prima nel calcolare la media campione,
:<math>\bar x = \displaystyle \frac {\sum_{j=1}^n x_j}{n}</math>,
e poi calcolare la somma dei quadrati delle differenze dalla media,
:<math>\mathrm{varianza} = s^2 = \displaystyle\frac {\sum_{i=1}^n (x_i - \bar x)^2}{n-1} \!</math>,
dove ''s'' è la deviazione standard.  Ciò è dato dal seguente [[codice sorgente]]:

<syntaxhighlight lang="python">
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
</syntaxhighlight>

Questo algoritmo è numericamente stabile se ''n'' è piccolo.<ref name="Einarsson2005"/><ref>{{Cita libro|nome=Nicholas |cognome=Higham |titolo=Accuracy and Stability of Numerical Algorithms (2 ed) (Problem 1.10)|editore=SIAM|anno=2002}}</ref> 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 [[Algoritmo di sommatoria di Kahan|somma compensata]] possono essere usate per ridurre questo errore.

=== Variante compensata ===
La versione con somma compensata dell'algoritmo visto sopra è:<ref name=":0" /><!--Nota da voce inglese: Da dove viene questo algoritmo? Non è una forma normale per una somma di Kahan.-->

<syntaxhighlight lang="python">
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
</syntaxhighlight>

== Algoritmo in linea ==

Spesso è utile essere in grado di calcolare la varianza in un passo singolo, elaborando ogni valore <math>x_i</math> 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 (statistica)|media]] e la varianza (stimata) della sequenza, per un elemento addizionale <math>x_{\mathrm{new}}</math>. Qui, <math>\bar x_n </math> denota la media campione dei primi <math>n</math> campioni <math>(x_1,\, \ldots ,\,x_n)</math>, <math>s^2_n</math> la loro varianza campione e <math>\sigma^2_n</math> la varianza della loro popolazione.

:<math>\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} \!</math>

:<math>s^2_n = \frac{(n-2)}{(n-1)} \, s^2_{n-1} + \frac{(x_n - \bar x_{n-1})^2}{n}, \quad n>1 </math>

:<math>\sigma^2_n = \frac{(n-1) \, \sigma^2_{n-1} + (x_n - \bar x_{n-1})(x_n - \bar x_{n})}{n}.</math>

Queste formule soffrono di instabilità numerica. Una quantità migliore per l'aggiornamento è la somma dei quadrati delle differenze dalla media corrente <math>\textstyle\sum_{i=1}^n (x_i - \bar x_n)^2</math>, qui indicate come <math>M_{2,n}</math>:

:<math>M_{2,n}\! = M_{2,n-1} + (x_n - \bar x_{n-1})(x_n - \bar x_n)</math>
:<math>s^2_n = \frac{M_{2,n}}{n-1}</math>
:<math>\sigma^2_n = \frac{M_{2,n}}{n}</math>

Un algoritmo numericamente stabile per la varianza campione è indicato sotto. Esso calcola anche la media.
Questo algoritmo è dovuto a Knuth,<ref>[[Donald E. Knuth]] (1998). ''[[The Art of Computer Programming]]'', volume 2: ''Seminumerical Algorithms'', 3rd edn., p. 232. Boston: Addison-Wesley.</ref> che cita Welford,<ref>B. P. Welford (1962).[https://www.jstor.org/stable/1266577 "Note on a method for calculating corrected sums of squares and products"]. ''Technometrics'' 4(3):419–420.</ref> ed è stato accuratamente analizzato.<ref>Chan, Tony F.; Golub, Gene H.; LeVeque, Randall J. (1983). Algorithms for Computing the Sample Variance: Analysis and Recommendations. The American Statistician 37, 242-247. https://www.jstor.org/stable/2683386</ref><ref>Ling, Robert F. (1974). Comparison of Several Algorithms for Computing Sample Means and Variances. Journal of the American Statistical Association, Vol. 69, No. 348, 859-866. {{doi|10.2307/2286154}}</ref> È anche comune denotare <math>M_k = \bar x_k</math> e <math>S_k = M_{2,k}</math>.<ref>http://www.johndcook.com/standard_deviation.html</ref>

<syntaxhighlight lang="python">
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)
</syntaxhighlight>

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'[[Algoritmi per il calcolo della varianza#Algoritmo parallelo|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 <math>n</math> con la somma dei pesi vista finora.  West (1979)<ref>D. H. D. West (1979). ''Communications of the ACM'', 22, 9, 532-535: ''Updating Mean and Variance Estimates: An Improved Method''</ref> suggerisce questo algoritmo incrementale:

<syntaxhighlight lang="python">
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)
</syntaxhighlight>

== Algoritmo parallelo ==
Chan e altri.<ref name=":0">{{Cita pubblicazione|cognome1= Chan     |nome1= Tony F.      |wkautore1= Tony F. Chan|cognome2= Golub    |nome2= Gene H.      |wkautore2= Gene H. Golub|cognome3= LeVeque  |nome3= Randall J.|contributo= Updating Formulae and a Pairwise Algorithm for Computing Sample Variances.|titolo= Technical Report STAN-CS-79-773|editore= Department of Computer Science, Stanford University|anno= 1979|urlcontributo=http://i.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf }}</ref> Notare che l'algoritmo in linea indicato sopra è un caso particolare di algoritmo che funziona per ogni partizione del campione <math>X</math> negli insiemi <math>X_A</math>, <math>X_B</math>:
:<math>\delta\! = \bar x_B - \bar x_A</math>
:<math>\bar x_X = \bar x_A + \delta\cdot\frac{n_B}{n_X}</math>
:<math>M_{2,X} = M_{2,A} + M_{2,B} + \delta^2\cdot\frac{n_A n_B}{n_X}</math>.
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 <math>n_A \approx n_B</math> ed entrambi sono grandi, poiché l'errore numerico in <math>\bar x_B - \bar x_A</math> non è ridotto in modo da essere nel caso <math>n_B = 1</math>. In tali casi, è preferibile che <math>\bar x_X = \frac{n_A \bar x_A + n_B \bar x_B}{n_A + n_B}</math>.

== Esempio ==
Supponiamo che tutte le operazioni in virgola mobile usino l'aritmetica con lo standard [[IEEE 754#Precisione doppia (64 bit)|IEEE 754 con precisione doppia]]. Consideriamo il campione <math>(4,\,7,\,13,\,16)</math> 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  <math>(10^8+4,\,10^8+7,\,10^8+13,\,10^8+16)</math>, 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 <math>29.33333333333333</math> 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 <math>(10^9+4,\, 10^9+7,\,10^9+13,\,10^9+16)</math>. Nuovamente, la varianza stimata della popolazione, di valore 30, è calcolata correttamente dall'algoritmo "a due passi", ma l'algoritmo "banale" ora calcola <math>-170.66666666666666</math>. 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<ref>{{Cita pubblicazione|cognome=Terriberry|nome=Timothy B.|anno=2007|titolo=Computing Higher-Order Moments Online|url=https://people.xiph.org/~tterribe/notes/homs.html|accesso=29 aprile 2019|urlarchivio=https://web.archive.org/web/20190408192147/https://people.xiph.org/~tterribe/notes/homs.html|urlmorto=sì}}</ref> estende le formule di Chan al calcolo del terzo e quarto momento centrale, necessarie per esempio quando si stimano [[Simmetria (statistica)|asimmetria]] e [[curtosi]]:
:<math>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_AM_{2,B} - n_BM_{2,A}}{n_X}</math>
:<math>\begin{align}
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_AM_{3,B} - n_BM_{3,A}}{n_X} \\
\end{align}</math>

Qui, gli <math>M_k</math> sono nuovamente le somme di potenze di differenze dalla la media <math>\Sigma(x - \overline{x})^k</math>, avendo
:asimmetria: <math>g_1 = \frac{\sqrt{n} M_3}{M_2^{3/2}},</math>
:cutrosi: <math>g_2 = \frac{n M_4}{M_2^2}-3.</math>
Per il caso incrementale (cioè, <math>B = \{x\}</math>), ciò si riduce a:
:<math>\delta\! = x - m</math>
:<math>m' = m + \frac{\delta}{n}</math>
:<math>M_2' = M_2 + \delta^2 \frac{ n-1}{n}</math>

:<math>M_3' = M_3 + \delta^3 \frac{ (n - 1) (n - 2)}{n^2} - \frac{3\delta M_2}{n}</math>

:<math>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}</math>

Conservando il valore <math>\delta / n</math>, è 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 è:
<syntaxhighlight lang="python">
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
</syntaxhighlight>

Pébay<ref>{{Cita pubblicazione|cognome=Pébay|nome=Philippe|anno=2008|contributo=Formulas for Robust, One-Pass Parallel Computation of Covariances and Arbitrary-Order Statistical Moments|titolo=Technical Report SAND2008-6212|editore=Sandia National Laboratories|url=http://infoserve.sandia.gov/sand_doc/2008/086212.pdf|urlmorto=sì|accesso=17 giugno 2016|dataarchivio=9 ottobre 2022|urlarchivio=https://ghostarchive.org/archive/20221009/http://infoserve.sandia.gov/sand_doc/2008/086212.pdf}}</ref>
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 (probabilità)|covarianza]].

Choi e Sweetman
<ref name="Choi2010">
{{Cita pubblicazione|cognome1     = Choi|nome1        = Muenkeun|cognome2     = Sweetman|nome2        = Bert|anno = 2010|titolo = Efficient Calculation of Statistical Moments for Structural Health Monitoring|url = http://www.rms-group.org/RMS_Papers/TAMUG_Papers/MK/Efficient_Moments_2010.pdf|urlmorto = sì|accesso = 17 giugno 2016|urlarchivio = https://web.archive.org/web/20160303181457/http://www.rms-group.org/RMS_Papers/TAMUG_Papers/MK/Efficient_Moments_2010.pdf}}
</ref>
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:

: <math> H(x_k)=\frac{h(x_k)}{A}</math>

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

: <math>
 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
</math>

: <math>
 \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
</math>

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

: <math>
 m_n^{(h)}= \frac{1}{I} {\sum_{k=1}^{K} x_k^n \, h(x_k)}
</math>

: <math>
 \theta_n^{(h)}= \frac{1}{I}{\sum_{k=1}^{K} \Big(x_k-m_1^{(h)}\Big)^n \, h(x_k)}
</math>

Il secondo approccio, da Choi and Sweetman
<ref name="Choi2010" />,
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 <math>Q</math> di momenti statistici sono noti:
<math>(\gamma_{0,q},\mu_{q},\sigma^2_{q},\alpha_{3,q},\alpha_{4,q})
\quad </math> per <math>q=1,2,...,Q </math>, allora ogni <math>\gamma_n</math> può essere espresso in termini dei momenti di origine zero (raw moments) equivalenti <math>n</math> :

: <math>
 \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
</math>

dove, generalmente, come <math>\gamma_{0,q}</math> si prende la durata della storia temporale <math>q^{th}</math>, o il numero di punti se <math>\Delta t</math> è costante.

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

: <math>
 \gamma_{n,c}= \sum_{q=1}^{Q}\gamma_{n,q} \quad \quad \textrm{for} \quad n=0,1,2,3,4
</math>

dove la sovrascritta <math>_c</math> rappresenta la storia temporale <math>\gamma</math> concatenata o combinata. Questi valori combinati di
<math>\gamma</math> possono essere trasformati inversamente in momenti di origine zero (raw moments) rappresentando la storia temporale concatenata completa

: <math>
 m_{n,c}=\frac{\gamma_{n,c}}{\gamma_{0,c}} \quad \textrm{for} \quad n=1,2,3,4
</math>

relazioni note tra i momenti di origine zero (raw moments) (<math>m_n</math>) e il momento centrale <math>\left( \theta_n = E[(x-\mu)^n]\right)</math>
sono allora usate per calcolare i momenti centrali della storia temporale concatenata. Infine, i momenti statistici della storia concatenata vengono calcolati dai momenti centrali:

: <math>
 \mu_c=m_{1,c}
 \ \ \ \ \ \sigma^2_c=\theta_{2,c}
 \ \ \ \ \ \alpha_{3,c}=\frac{\theta_{3,c}}{\sigma_c^3}
 \ \ \ \ \ \alpha_{4,c}={\frac{\theta_{4,c}}{\sigma_c^4}}-3
</math>

== Covarianza ==
Si possono utilizzare algoritmi molto simili per calcolare la [[covarianza (probabilità)|covarianza]]. L'algoritmo "banale" è:
:<math>\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}. \!</math>

Per l'algoritmo scritto sopra, si può utilizzare il seguente codice Python:
<syntaxhighlight lang="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
</syntaxhighlight>

Come per la variaznza, la covarianza di due variabili casuali è anche invariante per traslazione, allora, dati due valori costanti <math>K_x</math> e <math>K_y</math> qualsiasi, si può scrivere:

:<math>\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}. \!</math>

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:

<syntaxhighlight lang="python">
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
</syntaxhighlight>

L'algoritmo a due passi prima calcola le medie campione e dopo la covarianza:
:<math>\bar x = \displaystyle \sum_{i=1}^n x_i/n</math>
:<math>\bar y = \displaystyle \sum_{i=1}^n y_i/n</math>
:<math>\operatorname{Cov}(X,Y) = \displaystyle\frac {\sum_{i=1}^n (x_i - \bar x)(y_i - \bar y)}{n}. \!</math>

L'algoritmo a due passi può essere scritto come segue:
<syntaxhighlight lang="python">
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
</syntaxhighlight>

Una versione compensata un po' più accurata esegue l'algoritmo "banale" completo sui residui. Le somme finali <math>\textstyle\sum x_i</math> e <math>\textstyle\sum y_i</math> 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:

<syntaxhighlight lang="python">
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
</syntaxhighlight>
Esiste un algoritmo stabile a un passo, simile a quello visto sopra, che calcola il co-momento <math>\textstyle C_n = \sum_{i=1}^n (x_i - \bar x_n)(y_i - \bar y_n)</math>:
:<math>\bar x_n = \bar x_{n-1} + \frac{x_n - \bar x_{n-1}}{n} \!</math>
:<math>\bar y_n = \bar y_{n-1} + \frac{y_n - \bar y_{n-1}}{n} \!</math>
:<math>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})</math>
L'asimmetria apparente in quest'ultima equazione è dovuta al fatto che <math>\textstyle (x_n - \bar x_n) = \frac{n-1}{n}(x_n - \bar x_{n-1})</math>, così entrambi i termini di aggiornamento sono pari a <math>\textstyle \frac{n-1}{n}(x_n - \bar x_{n-1})(y_n - \bar y_{n-1})</math>. 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
:<math>\begin{align}
\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{align}</math>

Allo stesso modo, vi è una formula per combinare le covarianze di due insiemi che può essere utilizzata per parallelizzare calcolo:
:<math>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}</math>.

== Note ==
<references />

== Voci correlate ==
*[[Formula computazionale per la varianza]]
*[[Algoritmo di sommatoria di Kahan]]

== Collegamenti esterni ==
* {{Collegamenti esterni}}

{{portale|statistica}}

[[Categoria:Statistica computazionale]]