superponer - ¿Qué algoritmo usa R para calcular la media?

superponer graficas en r

Question

Tengo curiosidad por saber qué función media del algoritmo R utiliza. ¿Hay alguna referencia a las propiedades numéricas de este algoritmo?

Encontré el siguiente código C en 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;

Parece hacer una media recta:

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

Luego agrega lo que asumo es una corrección numérica que parece ser la diferencia media con la media de los datos:

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

No he podido rastrear este algoritmo en ninguna parte (la media no es un gran término de búsqueda).

Cualquier ayuda sería muy apreciada.

Answer 1

Creo que el algoritmo R funciona de la siguiente manera.

El primer cálculo estándar de la media es efectivamente una estimación de la media algebraica, debido a errores de punto flotante (que se agrava a medida que la suma se aleja de los elementos que se acumulan).

La segunda pasada suma las diferencias de los elementos de la media estimada. No debería haber diferencia neta, ya que los valores a ambos lados de la media deberían compensarse, pero tenemos un error de punto flotante. Las diferencias de la media todavía tienen el potencial de error, pero deberían ser más pequeñas que la peor diferencia de potencial entre un elemento y la suma acumulada (al menos la media estimada vive en algún lugar dentro del rango de valores, mientras que la suma puede escapar de ella) . Dividir por N le da la diferencia promedio de la media, que luego utiliza para acercar su estimación inicial a la media real. Podría repetir esto para acercarse más y más, pero en algún punto el error de punto flotante al calcular la diferencia promedio de la media lo derrotará. Supongo que una pasada está lo suficientemente cerca.

Esto me lo explicó mi esposa.

No estoy seguro de cuál es la fuente del algoritmo, y no estoy seguro de cómo se compara esto con otros métodos, como la suma de Kahan. Supongo que tendré que hacer algunas pruebas.

Answer 2

No estoy seguro de qué algoritmo es este, pero Martin Maechler mencionó el método de actualización de West, 1979 en respuesta a PR#1228 , que fue implementado por Brian Ripley en R-2.3.0. No pude encontrar una referencia en el código fuente o en los registros de control de versión que enumeraron el algoritmo real utilizado. Se implementó en cov.c en la revisión 37389 y en summary.c en la revisión 37393.