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 è:
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 è:
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é e 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
con costante qualunque, il che conduce alla nuova formula
e più vicino è 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 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 , 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 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, , e poi calcolare la somma dei quadrati delle differenze dalla media, , 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 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 . Qui, denota la media campione dei primi campioni , la loro varianza campione e la varianza della loro popolazione.
Queste formule soffrono di instabilità numerica. Una quantità migliore per l'aggiornamento è la somma dei quadrati delle differenze dalla media corrente , qui indicate come :
Un algoritmo numericamente stabile per la varianza campione è indicato sotto. Esso calcola anche la media. Questo algoritmo è dovuto a Knuth,Donald E. Knuth (1998). The Art of Computer Programming, volume 2: Seminumerical Algorithms, 3rd edn., p. 232. Boston: Addison-Wesley. che cita Welford,B. P. Welford (1962)."Note on a method for calculating corrected sums of squares and products". Technometrics 4(3):419–420. ed è stato accuratamente analizzato.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/2683386Ling, 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. È anche comune denotare e .http://www.johndcook.com/standard_deviation.html
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 con la somma dei pesi vista finora. West (1979)D. H. D. West (1979). Communications of the ACM, 22, 9, 532-535: Updating Mean and Variance Estimates: An Improved Method 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 negli insiemi , :
. 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 ed entrambi sono grandi, poiché l'errore numerico in non è ridotto in modo da essere nel caso . In tali casi, è preferibile che .
Esempio Supponiamo che tutte le operazioni in virgola mobile usino l'aritmetica con lo standard IEEE 754 con precisione doppia. Consideriamo il campione 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 , 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 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 . Nuovamente, la varianza stimata della popolazione, di valore 30, è calcolata correttamente dall'algoritmo "a due passi", ma l'algoritmo "banale" ora calcola . 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:
Qui, gli sono nuovamente le somme di potenze di differenze dalla la media , avendo asimmetria: cutrosi: Per il caso incrementale (cioè, ), ciò si riduce a:
Conservando il valore , è 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:
dove e rappresentano la frequenza e la frequenza relativa per il pacchetto e è l'area totale dell'istogramma. Dopo questa normalizzazione, i momenti di origine zero (raw moments) e i momenti centrali di possono essere calcolati dal relativo istogramma: