Подтвердить что ты не робот

Почему данные сортировки chisq.test в порядке убывания перед суммированием

Почему функция chisq.test в R сортирует данные перед суммированием в порядке убывания?

Данный код:

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

Если меня беспокоило числовую стабильность из-за использования арифметики float и захотело использовать какой-то простой для развертывания взлома, я бы сортировал данные в порядке возрастания перед суммированием, чтобы не добавлять в накопитель крошечное значение в большое значение (в чтобы избежать максимально возможного обрезки наименее значимых бит в результате).

Я просмотрел исходный код sum, но он не объяснил, почему передавать данные в порядке убывания на sum(). Что мне не хватает?

Пример:

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)))

Результат:

10000000000010996 10000000000011000

Когда мы сортируем данные в порядке возрастания, получаем правильный результат. Если мы сортируем данные в порядке убывания, мы получаем результат, который выключен на 4.

4b9b3361

ОТВЕТЫ

Ответ 1

EDIT: Книга " Точность и стабильность численных алгоритмов Николя Дж. Хайама" гласит, что

"при суммировании неотрицательных чисел рекурсивным суммированием увеличение упорядочивания - лучший порядок, в смысле наличия наименьшая априорная ошибка ошибки."

Спасибо @Lamia за обмен книгой в разделе комментариев.

В этих книгах объясняются три метода суммирования, такие как рекурсивные, вставные и попарные методы. Каждый метод имеет свои достоинства и недостатки, основанные на величине связанной с ними привязки ошибок, которые могут быть вычислены путем систематического анализа ошибок суммирования чисел с плавающей запятой.

Примечательно, что результат суммирования из рекурсивного метода зависит от стратегий упорядочения, таких как увеличение, уменьшение и Psum (смотрите в книге - стр. 82 - 4-й абзац, а также посмотрите, как это работает в примере, приведенном внизу страницы 82.).

Глядя на исходный код R для функции sum(), которую можно получить из summary.c, сообщает, что R реализует рекурсивный метод в своем sum().

Также число базовых цифр в значении с плавающей запятой равно 53, которое можно получить из

.Machine$double.digits
# [1] 53

Установив это число как прецизионные биты, мы можем сравнить точность операции суммирования по базе R и mpfr() с библиотекой Rmpfr для разных стратегий упорядочения. Обратите внимание, что увеличение порядка дает результаты, близкие к тем, которые наблюдаются при суммировании с плавающей точкой, что подтверждает вышеупомянутое утверждение из этой книги.

Я вычислил статистику хи-квадратов с использованием исходных данных 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

Взгляд на исходный код функции edit(chisq.test) сообщает, что в нем нет операции сортировки.

Кроме того, как указано в разделе комментариев, оно не связано со знаком (+ ve или -ve) значений необработанных данных, используемых в функции chisq.test(). Эта функция не принимает отрицательные значения, поэтому она будет выплевывать ошибку, останавливая функцию с помощью этого сообщения "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

Разница в значениях при суммировании чисел с плавающей запятой связана с арифметикой с двойной точностью. См. Демонстрацию ниже. Когда точность чисел с плавающей запятой поддерживается на 200 цифр с помощью функции mpfr(), доступной в пакете Rmpfr, операция суммирования дает одинаковые результаты независимо от порядка вектора x1 или x2. Однако неравные значения наблюдаются, когда точность с плавающей точкой не поддерживается.

Нет точности 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:

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

Наименьшее положительное число с плавающей запятой в базе R может быть получено из приведенного ниже кода, и любое число, меньшее этого числа, будет усечено в базе R, что дает различные результаты из операции суммирования.

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

Сравнение функций арифметики с двойной точностью и неважных функций chisq.test().

Выделяется соответствующая часть chisq.test(), и с ней создается новая функция chisq.test2(). Внутри вы увидите опции для сравнения до и после применения двухзначной осведомленности о значении 250 цифр с использованием функции mpfr() для статистики хи-квадратов. Вы можете видеть, что идентичные результаты видны для функции с плавающей точкой, но не для исходных данных.

# 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)

Вывод:

# [[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