Quale algoritmo usa R per calcolare la media?

Question

Sono curioso di sapere quale algoritmo utilizza la funzione media di R. C'è qualche riferimento alle proprietà numeriche di questo algoritmo?

ho trovato il seguente codice C in summary.c: do_summary():

case REALSXP: 
PROTECT(ans = allocVector(REALSXP, 1)); 
for (i = 0; i < n; i++) s += REAL(x)[i]; 
s /= n; 
if(R_FINITE((double)s)) { 
    for (i = 0; i < n; i++) t += (REAL(x)[i] - s); 
    s += t/n; 
} 
REAL(ans)[0] = s; 
break; 

Sembra di fare un dritto up medio:

for (i = 0; i < n; i++) s += REAL(x)[i]; 
s /= n; 

Poi si aggiunge quello che presumo è un correzione numerica che sembra essere la differenza media dalla media dei dati:

for (i = 0; i < n; i++) t += (REAL(x)[i] - s); 
s += t/n; 

Non sono stato in grado di rintracciare s algoritmo giù ovunque (significa che non è un grande termine di ricerca).

Qualsiasi aiuto sarebbe molto apprezzato.

Answer 1

Non sono sicuro di quale algoritmo sia, ma Martin Maechler ha menzionato il metodo di aggiornamento di West, 1979 in risposta a PR#1228, che è stato implementato da Brian Ripley in R-2.3.0. Non sono riuscito a trovare un riferimento nel codice sorgente o nei registri di controllo della versione che elencavano l'algoritmo utilizzato. È stato implementato in cov.c nella revisione 37389 e in summary.c nella revisione 37393.

Answer 2

Credo che l'algoritmo R funzioni come segue.

Il primo calcolo standard della media è effettivamente una stima della media algebrica, a causa di errori in virgola mobile (che peggiora ulteriormente man mano che la somma si allontana dagli elementi accumulati).

Il secondo passaggio somma le differenze degli elementi dalla media stimata. Non ci dovrebbero essere differenze nette in quanto i valori su entrambi i lati della media dovrebbero essere bilanciati, ma abbiamo un errore in virgola mobile. Le differenze rispetto alla media hanno ancora il potenziale di errore, ma queste dovrebbero essere più piccole della peggiore differenza potenziale tra un elemento e la somma accumulata (almeno la media stimata vive da qualche parte all'interno dell'intervallo di valori, mentre la somma può sfuggire ad essa) . Dividere per N ti dà la differenza media dalla media, che poi usi per spingere la tua stima iniziale più vicina alla media vera. Potresti ripetere questo per avvicinarti sempre di più, ma ad un certo punto l'errore del punto in virgola mobile nel calcolare la differenza media dalla media ti sconfiggerà. Immagino che un passaggio sia abbastanza vicino.

Questo mi è stato spiegato da mia moglie.

Non sono sicuro di dove sia la fonte dell'algoritmo e non sono sicuro di come questo si colleghi ad altri metodi, come la sommatoria di Kahan. Credo che dovrò fare qualche test.