goodness - ¿Por qué chisq.test ordena los datos en orden descendente antes de la suma?

chisq.test in r

Question

¿Por chisq.test función chisq.test en R ordena los datos antes de la suma en orden descendente ?

El código en cuestión es:

STATISTIC <- sum(sort((x - E)^2/E, decreasing = TRUE))

Si estuviera preocupado por la estabilidad numérica debido al uso de la aritmética flotante y quisiera utilizar un hackeo fácil de implementar, clasificaría los datos en orden creciente antes de la suma para evitar agregar un valor pequeño al valor grande en el acumulador (para evitar Recorte de los bits menos significativos en el resultado tanto como sea posible).

Busqué en el código fuente de la sum pero no explicó por qué pasar los datos en orden descendente a la sum() . ¿Qué extraño?

Un ejemplo:

x = matrix(1.1, 10001, 1)
x[1] = 10^16   # We have a vector with 10000*1.1 and 1*10^16
c(sum(sort(x, decreasing = TRUE)), sum(sort(x, decreasing = FALSE)))

El resultado:

10000000000010996 10000000000011000

Cuando ordenamos los datos en orden ascendente, obtenemos el resultado correcto. Si ordenamos los datos en orden descendente, obtenemos un resultado que está desactivado en 4.

Answer 1

EDITAR: El libro " Precisión y estabilidad de los algoritmos numéricos por Nicolas J. Higham " afirma que

"al sumar números no negativos por suma recursiva, el orden creciente es el mejor ordenamiento, en el sentido de tener el límite de error directo a priori más pequeño".

Gracias a @Lamia por compartir el libro en la sección de comentarios.

Este libro explica tres métodos de resumen, como las técnicas recursivas, de inserción y de pares. Cada técnica tiene sus propios méritos y deméritos en función de la magnitud del límite de error asociado a ellos, que se puede calcular haciendo un análisis sistemático de errores de la suma de los números de punto flotante.

En particular, los resultados de la suma de la técnica recursiva dependen de las estrategias de ordenación tales como el aumento, la disminución y el Psum (mire en el libro - página 82 - cuarto párrafo. También vea cómo funciona en el ejemplo que se encuentra al final de la página 82.) .

Al observar el código fuente de R para la función sum() que se puede obtener en summary.c informa que R implementa el método recursivo en su función sum() .

También el número de dígitos base en el significado de punto flotante es 53, que se puede obtener de

.Machine$double.digits
# [1] 53

Al establecer este número como los bits de precisión, podemos comparar la precisión de la operación de suma por base R y mpfr() de la biblioteca Rmpfr para diferentes estrategias de ordenación. Observe que el orden creciente produce resultados más cercanos a los que se ven en sumas conscientes de punto flotante, lo que corrobora la afirmación anterior de este libro.

Calculé la estadística de chi cuadrado usando los datos en bruto x .

library(''data.table'')
library(''Rmpfr'')
x1 = matrix(c( 10^16, rep(1.1, 10000)), 
            nrow = 10001, ncol = 1)
df1 <- data.frame(x = x1)
setDT(df1)
df1[, p := rep(1/length(x), length(x))]
s_x <- df1[, sum(x)]
df1[, E := s_x * p]
df1[, chi := ((x - E)^2/E)]

precBits <- .Machine$double.digits
x_chi <- data.frame( names = c("x_asc", "x_desc", "x_fp_asc", "x_fp_desc",
                               "chi_asc", "chi_desc", "chi_fp_asc", "chi_fp_desc"))
x_chi$vals <- c( ## x
  df1[order(x), format( sum(x), digits = 22)],
  df1[order(-x), format( sum(x), digits = 22)],
  df1[order(x), format( sum(mpfr(x, precBits = precBits)), digits = 22)],
  df1[order(-x), format( sum(mpfr(x, precBits = precBits)), digits = 22)],
  ## chi
  df1[order(chi), format( sum(chi), digits = 22)],
  df1[order(-chi), format( sum(chi), digits = 22)],
  df1[order(chi), format( sum(mpfr(chi, precBits = precBits)), digits = 22)],
  df1[order(-chi), format( sum(mpfr(chi, precBits = precBits)), digits = 22)])

x_chi
#         names                    vals
# 1       x_asc       10000000000011000
# 2      x_desc       10000000000010996
# 3    x_fp_asc 10000000000011000.00000
# 4   x_fp_desc 10000000000020000.00000
# 5     chi_asc    99999999999890014218
# 6    chi_desc    99999999999890030592
# 7  chi_fp_asc 99999999999890014208.00
# 8 chi_fp_desc 99999999999833554944.00

Al mirar el código fuente de la función de edit(chisq.test) informa que no hay ninguna operación de clasificación involucrada en él.

Además, como se señaló en la sección de comentarios, no está relacionado con el signo (+ ve o -ve) de los valores de los datos sin procesar utilizados en la función chisq.test() . Esta función no acepta valores negativos, por lo que escupirá el error al detener la función con este mensaje "all entries of ''x'' must be nonnegative and finite" .

set.seed(2L)
chisq.test(c(rnorm(10, 0, 1)))
# Error in chisq.test(c(rnorm(10, 0, 1))) : 
#   all entries of ''x'' must be nonnegative and finite

La diferencia en los valores al sumar números de punto flotante se refiere a la aritmética de doble precisión. Vea la demostración a continuación. Cuando la precisión de los números de punto flotante se mantiene en 200 dígitos utilizando la función mpfr() disponible en el paquete Rmpfr , la operación de suma da resultados idénticos independientemente del orden del vector x1 o x2 . Sin embargo, se observan valores desiguales cuando no se mantiene una precisión de punto flotante.

Sin precisión de FP:

x1 = matrix(c( 10^16, rep(1.1, 10000)), 
            nrow = 10001, ncol = 1)
## reverse
x2 = matrix(c( rep(1.1, 10000), 10^16 ), 
            nrow = 10001, ncol = 1)

c( format(sum(x1), digits = 22), 
   format(sum(x2), digits = 22))
# [1] "10000000000010996" "10000000000011000"

FP Precisión mantenida:

library(''Rmpfr'')
##
prec <- 200
x1 = matrix(c( mpfr( 10^16, precBits = prec),
              rep( mpfr(1.1, precBits = prec), 10000)), 
           nrow = 10001, ncol = 1)

## reverse
x2 = matrix(c( rep(mpfr(1.1, precBits = prec), 10000), 
              mpfr( 10^16, precBits = prec) ), 
           nrow = 10001, ncol = 1)
c( sum(x1), sum(x2))
# 2 ''mpfr'' numbers of precision  200   bits 
# [1] 10000000000011000.000000000000888178419700125232338905334472656
# [2] 10000000000011000.000000000000888178419700125232338905334472656

El número de punto flotante positivo más pequeño en la base R se puede obtener a partir del código siguiente y cualquier número menor que este número se truncará en la base R, lo que produce resultados variables de la operación de suma.

.Machine$double.eps
# [1] 2.220446e-16

Comparación de funciones aritméticas de doble precisión e inconscientes de la función chisq.test() .

Se chisq.test() parte relevante de chisq.test() y se hace una nueva función chisq.test2() . En el interior, verá opciones de comparación antes y después de aplicar el mpfr() doble precisión de 250 dígitos utilizando la función mpfr() para la estadística de chi cuadrado. Puede ver resultados idénticos para la función de reconocimiento de punto flotante, pero no para los datos sin procesar.

# modified chi square function:
chisq.test2 <- function (x, precBits) 
{
  if (is.matrix(x)) {
    if (min(dim(x)) == 1L) 
      x <- as.vector(x)
  }

  #before fp precision
  p = rep(1/length(x), length(x))
  n <- sum(x)
  E <- n * p

  # after fp precision
  x1 <- mpfr(x, precBits = precBits)
  p1 = rep(1/length(x1), length(x1))
  n1 <- sum(x1)
  E1 <- n1 * p1

  # chisquare statistic
  STATISTIC <- c(format(sum((x - E)^2/E), digits=22),           # before decreasing
                 format(sum(sort((x - E)^2/E, decreasing = FALSE)), digits=22), # before increasing
                 sum((x1 - E1)^2/E1),                           # after decreasing 
                 sum(sort((x1 - E1)^2/E1, decreasing = FALSE))) # after increasing

  return(STATISTIC)
}

# data
x1 = matrix(c( 10^16, rep(1.1, 10000)), 
            nrow = 10001, ncol = 1)

chisq.test2(x = x1, precBits=250)

Salida:

# [[1]]  # before fp decreasing
# [1] "99999999999890030592"
# 
# [[2]]  # before fp increasing
# [1] "99999999999890014218"
# 
# [[3]]  # after fp decreasing 
# ''mpfr1'' 99999999999889972569.502489584522352514811399898444554440067408531548230046685
# 
# [[4]]  # after fp increasing
# ''mpfr1'' 99999999999889972569.502489584522352514811399898444554440067408531548230251906