Además flotante punto: los problemas de pérdida de precisión

Question

En pocas palabras: ¿cómo puedo ejecutar a+b tal que cualquier pérdida de precisión debido al truncamiento es lejos de cero en lugar de hacia el cero

La larga historia

Estoy calculando la suma de una larga serie de valores de punto flotante para el propósito de calcular la media muestral y la varianza del conjunto. Desde Var (X) = E (X ² ) - E (X) ² , basta para mantener funcionando el recuento de todos los números, el suma de todos los números hasta el momento, y la suma de los cuadrados de todos los números hasta el momento.

Hasta aquí todo bien.

Sin embargo, es absolutamente necesario que la E (X ² )> E (X) ² , que debido a ISN precisión de punto flotante' t siempre el caso. En pseudo-código, el problema es el siguiente:

int count;
double sum, sumOfSquares;
...
double value = <current-value>;
double sqrVal = value*value; 

count++;
sum += value; //slightly rounded down since value is truncated to fit into sum
sumOfSquares += sqrVal; //rounded down MORE since the order-of-magnitude 
//difference between sqrVal and sumOfSquares is twice that between value and sum;

Para las secuencias variables, esto no es un gran problema - se termina ligeramente por debajo de estimación de la varianza, pero a menudo no es un gran problema. Sin embargo, para los conjuntos de constantes o casi constantes con un no-media cero, puede significar que E (X ² ) 2 , lo que resulta en una variación negativa computarizada, lo que viola las expectativas de consumo de código.

Ahora, sé de Kahan La suma, que no es una solución atractiva. En primer lugar, hace que el código susceptible a los caprichos de optimización (en función de parámetros de optimización, el código puede o no puede presentar este problema), y en segundo lugar, el problema no es realmente debido a la precisión - que es bueno lo suficiente - que es por adición introduce sistemática de error hacia cero. Si pudiera ejecutar la línea

sumOfSquares += sqrVal;

de tal manera que se garantice que sqrVal se redondea hacia arriba, no hacia abajo, en la precisión de SumOfSquares, tendría una solución numéricamente razonable. Pero, ¿cómo puedo lograr eso?

_{Editar:? finalizados pregunta - ¿Por qué pulsando Entrar en la lista extensible gota en el campo de la etiqueta presentar la pregunta de todos modos}

Answer 1

Hay otro algoritmo de paso único que reorganiza el cálculo un poco. En pseudocódigo:

n = 0
mean = 0
M2 = 0

for x in data:
    n = n + 1
    delta = x - mean
    mean = mean + delta/n
    M2 = M2 + delta*(x - mean)  # This expression uses the new value of mean

variance_n = M2/n         # Sample variance
variance = M2/(n - 1)     # Unbiased estimate of population variance

(Fuente: http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance )

Esto parece mejor comportamiento con respecto a los problemas que señaló con el algoritmo habitual.

Answer 2

IEEE proporciona cuatro modos de redondeo, (hacia -inf, hacia + inf, hacia 0, tonearest). Hacia + inf es lo que parece querer. No hay control estándar en C90 o C ++. C99 añade el <fenv.h> encabezado que también está presente como una extensión de alguna C90 y aplicación C ++. Para respetar la norma C99, que tendría que escribir algo como:

#include <fenv.h>
#pragma STDC FENV_ACCESS ON

int old_round_mode = fegetround();
int set_round_ok = fesetround(FE_UPWARD);
assert(set_round_ok == 0);
...
int set_round_ok = fesetround(old_round_mode);
assert(set_round_ok == 0);

Es bien sabido que el algoritmo utiliza es numéricamente inestable y tiene un problema de precisión. Es mejor para la precisión de hacer dos pasadas sobre los datos.

Answer 3

Si usted no se preocupe por la precisión, pero sólo de una variación negativa, ¿por qué no sólo tiene que hacer V(x) = Max(0, E(X^2) - E(X)^2)