Algoritmos “en línea” (iterador) para estimar la mediana estadística, modo, asimetría, curtosis?

¿Existe un algoritmo para estimar la mediana, el modo, la asimetría y / o la curtosis del conjunto de valores, pero eso NO requiere almacenar todos los valores en la memoria a la vez?

Me gustaría calcular las estadísticas básicas:

  • media: media aritmética
  • varianza: promedio de desviaciones al cuadrado de la media
  • desviación estándar: raíz cuadrada de la varianza
  • mediana: valor que separa la mitad más grande de los números de la mitad más pequeña
  • modo: valor más frecuente encontrado en el conjunto
  • asimetría: tl; Dr
  • curtosis: tl; Dr

La fórmula básica para calcular cualquiera de estos es la aritmética de la escuela primaria, y los conozco. También hay muchas bibliotecas de estadísticas que los implementan.

Mi problema es el gran número (miles de millones) de valores en los conjuntos que estoy manejando: trabajando en Python, no puedo hacer una lista o hash con miles de millones de elementos. Incluso si escribo esto en C, las matrices de mil millones de elementos no son demasiado prácticas.

Los datos no están ordenados Es producido al azar, sobre la marcha, por otros procesos. El tamaño de cada conjunto es muy variable, y los tamaños no se conocerán de antemano.

Ya he descubierto cómo manejar la media y la varianza bastante bien, iterando a través de cada valor en el conjunto en cualquier orden. (En realidad, en mi caso, los tomo en el orden en que se generan). Aquí está el algoritmo que estoy usando, cortesía de http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#On-line_algorithm :

  • Inicializa tres variables: count, sum y sum_of_squares
  • Para cada valor:
    • Recuento de incrementos
    • Agregue el valor a la sum.
    • Agregue el cuadrado del valor a sum_of_squares.
  • Divide sum por conteo, almacenando como la variable media.
  • Divida sum_of_squares por count, almacenando como la variable mean_of_squares.
  • Media cuadrada, almacenando como square_of_mean.
  • Reste square_of_mean de mean_of_squares, almacenándolo como varianza.
  • Media de salida y varianza.

Este algoritmo “en línea” tiene debilidades (por ejemplo, problemas de precisión ya que sum_of_squares crece rápidamente más grande que el rango entero o la precisión de flotación), pero básicamente me da lo que necesito, sin tener que almacenar cada valor en cada conjunto.

Pero no sé si existen técnicas similares para estimar las estadísticas adicionales (mediana, modo, asimetría, curtosis). Podría vivir con un estimador sesgado, o incluso un método que comprometa la precisión hasta cierto punto, siempre que la memoria requerida para procesar los valores de N sea sustancialmente menor que O (N).

Indicarme una biblioteca de estadísticas existente también ayudará si la biblioteca tiene funciones para calcular una o más de estas operaciones “en línea”.

Sesgo y Curtosis

Para los algoritmos en línea para Skewness y Kurtosis (a lo largo de las líneas de la varianza), consulte en la misma página wiki los algoritmos paralelos para estadísticas de momentos más altos.

Mediana

La mediana es difícil sin datos ordenados. Si sabe cuántos puntos de datos tiene, en teoría solo tiene que clasificar parcialmente, por ejemplo, utilizando un algoritmo de selección . Sin embargo, eso no ayuda demasiado con miles de millones de valores. Sugeriría usar conteos de frecuencia, ver la siguiente sección.

Mediana y Modo con Frecuencia Cuenta

Si se trata de números enteros, contaría las frecuencias , probablemente cortando los valores más altos y más bajos más allá de algún valor, donde estoy seguro de que ya no es relevante. Para flotantes (o demasiados enteros), probablemente crearía cubos / intervalos, y luego usaría el mismo enfoque que para los enteros. (Aproximado) modo y cálculo de la mediana que se pone fácil, basado en la tabla de frecuencias.

Variables aleatorias distribuidas normalmente

Si se distribuye normalmente, utilizaría la media de muestra de la población, la varianza , la asimetría y la curtosis como estimadores de máxima verosimilitud para un pequeño subconjunto. Los algoritmos (en línea) para calcularlos, usted ya lo hace ahora. Por ejemplo, lea en un par de cientos de miles de millones de puntos de datos, hasta que su error de estimación sea lo suficientemente pequeño. Solo asegúrese de elegir al azar de su conjunto (por ejemplo, que no introduzca un sesgo al elegir los primeros 100’000 valores). El mismo enfoque también se puede usar para estimar el modo y la mediana para el caso normal (para ambos, la media de la muestra es un estimador).

Más comentarios

Todos los algoritmos anteriores se pueden ejecutar en paralelo (incluidos muchos algoritmos de clasificación y selección, por ejemplo, QuickSort y QuickSelect), si esto ayuda.

Siempre he supuesto (a excepción de la sección sobre la distribución normal) que hablamos de momentos de muestra, mediana y modo, no estimadores de momentos teóricos dada una distribución conocida.

En general, el muestreo de los datos (es decir, solo mirar un subconjunto) debería ser bastante exitoso dada la cantidad de datos, siempre y cuando todas las observaciones sean realizaciones de la misma variable aleatoria (tengan las mismas distribuciones) y los momentos, modo y la mediana realmente existe para esta distribución. La última advertencia no es inocua. Por ejemplo, la media (y todos los momentos superiores) para la distribución de Cauchy no existe. En este caso, la media muestral de un subconjunto “pequeño” podría estar masivamente fuera de la media muestral de la muestra completa.

Utilizo estos estimadores medianos y promedios incrementales / recursivos, que usan almacenamiento constante:

mean += eta * (sample - mean) median += eta * sgn(sample - median) 

donde eta es un parámetro de velocidad de aprendizaje pequeño (p. ej., 0.001) y sgn () es la función de signo que devuelve uno de {-1, 0, 1}. (Use una constante eta si los datos no son estacionarios y desea rastrear los cambios a lo largo del tiempo; de lo contrario, para las fonts estacionarias puede usar algo como eta = 1 / n para el estimador medio, donde n es el número de muestras vistas hasta ahora … desafortunadamente, esto no parece funcionar para el estimador mediano.)

Este tipo de estimador de promedios incrementales parece usarse en todas partes, por ejemplo, en reglas de aprendizaje de redes neuronales no supervisadas, pero la versión mediana parece ser mucho menos común, a pesar de sus beneficios (solidez a valores atípicos). Parece que la versión mediana podría usarse como reemplazo del estimador promedio en muchas aplicaciones.

Me encantaría ver un estimador de modo incremental de una forma similar …

ACTUALIZAR

Acabo de modificar el estimador mediano incremental para estimar cuantiles arbitrarios. En general, una función cuantil ( http://en.wikipedia.org/wiki/Quantile_function ) le indica el valor que divide los datos en dos fracciones: p y 1-p. A continuación, se estima este valor de forma incremental:

 quantile += eta * (sgn(sample - quantile) + 2.0 * p - 1.0) 

El valor p debe estar dentro de [0,1]. Esto esencialmente desplaza la salida simétrica {-1,0,1} de la función sgn () para inclinarse hacia un lado, dividiendo las muestras de datos en dos contenedores de tamaño desigual (las fracciones p y 1-p de los datos son menores que / mayores que la estimación del cuantil, respectivamente). Tenga en cuenta que para p = 0.5, esto se reduce al estimador mediano.

Implementé el algoritmo P-Square para el cálculo dynamic de cuantiles e histogtwigs sin almacenar observaciones en un módulo limpio de Python que escribí llamado LiveStats . Debería resolver su problema bastante efectivamente. La biblioteca admite todas las estadísticas que mencione, excepto el modo. Todavía no he encontrado una solución satisfactoria para la estimación de modo.

Ryan, me temo que no estás haciendo la media y la varianza correcta … Esto apareció hace unas semanas aquí . Y uno de los puntos fuertes de la versión en línea (que en realidad se conoce con el nombre del método de Welford) es el hecho de que es especialmente preciso y estable, consulte la discusión aquí . Uno de los puntos fuertes es el hecho de que no necesita almacenar la sum total o la sum total de cuadrados …

No puedo pensar en ningún enfoque en línea para el modo y la mediana, que parecen requerir considerar toda la lista a la vez. Pero bien podría ser que un enfoque similar al de la varianza y la media funcionará también para la asimetría y la curtosis …

El artículo de Wikipedia citado en la pregunta contiene las fórmulas para calcular la asimetría y la curtosis en línea.

Por modo, creo, no hay forma de hacerlo en línea. ¿Por qué? Suponga que todos los valores de su entrada son diferentes, además de la última que duplica una anterior. En este caso, debe recordar todos los valores ya vistos en la entrada para detectar que el último valor duplica un valor visto antes y lo convierte en el más frecuente.

Para la mediana, es casi lo mismo: hasta la última entrada no se sabe qué valor se convertirá en la mediana si todos los valores de entrada son diferentes porque podría ser antes o después de la mediana actual. Si conoce la longitud de la entrada, puede encontrar la mediana sin almacenar todos los valores en la memoria, pero aún tendrá que almacenar muchos de ellos (supongo que alrededor de la mitad) porque una secuencia de entrada incorrecta podría cambiar la mediana en gran medida en el la segunda mitad posiblemente haga cualquier valor desde la primera mitad de la mediana.

(Tenga en cuenta que solo me refiero al cálculo exacto).

Si tiene miles de millones de puntos de datos, entonces no es probable que necesite respuestas exactas, en lugar de respuestas cerradas. En general, si tiene miles de millones de puntos de datos, el proceso subyacente que los genera probablemente obedecerá a algún tipo de propiedad estadística de estacionariedad / ergodicidad / mezcla. También puede importar si espera que las distribuciones sean razonablemente continuas o no.

En estas circunstancias, existen algoritmos para on-line, poca memoria, estimación de cuantiles (la mediana es un caso especial de cuantil 0.5), así como modos, si no necesita respuestas exactas. Este es un campo activo de estadísticas.

Ejemplo de estimación de cuantiles: http://www.computer.org/portal/web/csdl/doi/10.1109/WSC.2006.323014

ejemplo de estimación de modo: Bickel DR. Estimadores robustos del modo y la asimetría de los datos continuos. Estadística computacional y análisis de datos. 2002; 39: 153 – 163. doi: 10.1016 / S0167-9473 (01) 00057-3.

Estos son campos activos de estadísticas computacionales. Está ingresando a los campos donde no existe un algoritmo exacto mejor, sino una diversidad de ellos (en realidad, estimadores estadísticos) que tienen diferentes propiedades, suposiciones y desempeño. Son las matemáticas experimentales. Probablemente haya cientos o miles de trabajos sobre el tema.

La última pregunta es si realmente necesita asimetría y curtosis por sí mismos, o más probablemente algunos otros parámetros que pueden ser más confiables para caracterizar la distribución de probabilidad (¡suponiendo que tenga una distribución de probabilidad!). ¿Estás esperando un gaussiano?

¿Tiene formas de limpiar / preprocesar los datos para que sea principalmente gaussiano? (por ejemplo, los montos de las transacciones financieras a menudo son algo gaussianos después de tomar los logaritmos). ¿Esperas desviaciones estándar finitas? ¿Esperas colas gordas? ¿Son las cantidades que te importan en las colas o en la masa?

Todo el mundo sigue diciendo que no se puede hacer el modo de una manera en línea, pero eso simplemente no es cierto. Aquí hay un artículo que describe un algoritmo para hacer exactamente este mismo problema inventado en 1982 por Michael E. Fischer y Steven L. Salzberg de la Universidad de Yale. Del artículo:

El algoritmo de búsqueda de la mayoría utiliza uno de sus registros para el almacenamiento temporal de un único elemento de la secuencia; este artículo es el candidato actual para el elemento de la mayoría. El segundo registro es un contador inicializado a 0. Para cada elemento de la secuencia, le pedimos al algoritmo que realice la siguiente rutina. Si el contador dice 0, instale el elemento de flujo actual como el nuevo candidato de la mayoría (desplazando cualquier otro elemento que pueda estar ya en el registro). Luego, si el elemento actual coincide con el candidato de la mayoría, incremente el contador; de lo contrario, disminuya el contador. En este punto del ciclo, si la parte del flujo vista hasta ahora tiene un elemento de mayoría, ese elemento está en el registro de candidato, y el contador tiene un valor mayor que 0. ¿Qué pasa si no hay un elemento de mayoría? Sin hacer una segunda pasada a través de los datos, lo cual no es posible en un entorno de flujo continuo, el algoritmo no siempre puede dar una respuesta inequívoca en esta circunstancia. Simplemente promete identificar correctamente el elemento de la mayoría, si es que hay alguno.

También se puede ampliar para encontrar la N superior con más memoria, pero esto debería resolverla para el modo.

En última instancia, si no tiene un conocimiento paramétrico a priori de la distribución, creo que debe almacenar todos los valores.

Dicho esto, a menos que esté tratando con algún tipo de situación patológica, el remedio (Rousseuw y Bassett 1990) bien puede ser lo suficientemente bueno para sus propósitos.

Simplemente implica calcular la mediana de lotes de medianas.

la mediana y el modo no se pueden calcular en línea usando solo el espacio constante disponible. Sin embargo, como la mediana y el modo son de todos modos más “descriptivos” que “cuantitativos”, puede estimarlos, por ejemplo, muestreando el conjunto de datos.

Si los datos se distribuyen normalmente en el largo plazo, entonces podrías usar tu media para estimar la mediana.

También puede estimar la mediana usando la siguiente técnica: establezca una estimación mediana M [i] para cada, digamos, 1,000,000 entradas en el flujo de datos de modo que M [0] sea la mediana del primer millón de entradas, M [1] el mediana del segundo millón de entradas, etc. Luego use la mediana de M [0] … M [k] como el estimador mediano. Esto, por supuesto, ahorra espacio, y puede controlar cuánto desea usar el espacio “ajustando” el parámetro 1,000,000. Esto también se puede generalizar recursivamente.

OK amigo prueba estos:

para c ++:

 double skew(double* v, unsigned long n){ double sigma = pow(svar(v, n), 0.5); double mu = avg(v, n); double* t; t = new double[n]; for(unsigned long i = 0; i < n; ++i){ t[i] = pow((v[i] - mu)/sigma, 3); } double ret = avg(t, n); delete [] t; return ret; } double kurt(double* v, double n){ double sigma = pow(svar(v, n), 0.5); double mu = avg(v, n); double* t; t = new double[n]; for(unsigned long i = 0; i < n; ++i){ t[i] = pow( ((v[i] - mu[i]) / sigma) , 4) - 3; } double ret = avg(t, n); delete [] t; return ret; } 

donde dice que ya puede calcular la varianza de la muestra (svar) y el promedio (prom), los señala a sus funciones para hacer eso.

Además, eche un vistazo a la aproximación de Pearson. en un conjunto de datos tan grande sería bastante similar. 3 (media - mediana) / desviación estándar tiene una mediana como max - min / 2

para el modo flotadores no tiene significado. uno típicamente los pegaría en contenedores de tamaño insignificante (como 1/100 * (máx. - mín.)).

Yo tendería a usar cubos, que podrían ser adaptables. El tamaño del cucharón debe ser la precisión que necesita. Luego, a medida que ingrese cada punto de datos, agregue uno al conteo del cubo correspondiente. Esto debería proporcionarle aproximaciones simples a la mediana y curtosis, contando cada cubo como su valor ponderado por su conteo.

El único problema podría ser la pérdida de resolución en punto flotante después de miles de millones de operaciones, es decir, ¡agregar una ya no cambia el valor! Para evitar esto, si el tamaño máximo de cubeta supera algún límite, puede quitar un número grande de todos los conteos.

 for j in range (1,M): y=np.zeros(M) # build the vector y y[0]=y0 #generate the white noise eps=npr.randn(M-1)*np.sqrt(var) #increment the y vector for k in range(1,T): y[k]=corr*y[k-1]+eps[k-1] yy[j]=y list.append(y)