взвешенная медиана в пакете spatstat - PullRequest
Новая
       89

взвешенная медиана в пакете spatstat

1 голос
/ 23 февраля 2020

Функция weighted.median() в пакете spatstat возвращает "10,5", когда я передаю равномерно взвешенные оценки 10, 11 и 12. Я ожидал ответ "11" (что вывод stats::median() и matrixStats::weightedMedian()).

Концепция взвешенной медианы не очень естественна для меня. Вывод неправильный или я неправильно понимаю назначение функции? 

x <- c(10, 11, 12)
w <- c( 1,  1,  1)

spatstat::weighted.median(x, w)
#> [1] 10.5
spatstat::weighted.quantile(x, w, probs = .5)
#>  50% 
#> 10.5


matrixStats::weightedMedian(x, w)
#> [1] 11
median(x)
#> [1] 11

Создано в 2020-02-23 с помощью пакета Представить (v0.3.0)

Ответы [ 2 ]

2 голосов
/ 24 февраля 2020

Здесь есть более фундаментальная проблема определения квантиля (включая медиану) в небольших конечных выборках.

В файле справки для базовой функции R quantile.default говорится, что существует аргумент type с 7 различными вариантами, которые будут давать разные ответы. Они тщательно описаны в прекрасной статье Роба Хиндмана, цитируемой в файле справки. По умолчанию для quantile.default установлено type=7.

Алгоритм в spatstat::weighted.quantile выполняет аналог type=4 (согласно его справочному файлу); то есть совокупная функция распределения F(x) линейно интерполируется, а затем вычисляется обратная функция. Этот алгоритм правильно реализован в коде spatstat.

Взвешенная медиана в другом упомянутом вами пакете вычисляет другое определение взвешенной медианы.

Спасибо большое много для привлечения нашего внимания к этому примеру. Это может побудить нас расширить реализацию spatstat::weighted.median, чтобы охватить другие типы.

Кстати, отчеты об ошибках для пакета CRAN должны публиковаться на странице отчетов об ошибках пакета, как показано в CRAN. Мне просто повезло, что я увидела этот пост. Но большое спасибо вам обоим за то, что вы обнаружили этот вопрос.

1 голос
/ 24 февраля 2020

Я считаю, что это недостаток в пакете, и я объясню почему.

Во-первых, weighted.median на самом деле просто вызывает weighted.quantile с вектором probs, установленным в 0.5. Но если вы позвоните weighted.quantile со своими данными, вы получите очень странные результаты: 

weighted.quantile(x, w)
#>    0%   25%   50%   75%  100% 
#> 10.00 10.00 10.50 11.25 12.00

Это неправильно.

Если вы посмотрите на тело этой функции, используя body(weighted.quantile), и следуйте логике c до конца, кажется, есть проблема с тем, как веса нормализуются в строке 10 в переменную с именем Fx. Для правильной работы нормализованные веса должны иметь вектор такой же длины, что и x, но начиная с 0 и заканчивая 1, с интервалом между ними, пропорциональным весам.

Но если вы посмотрите как это фактически вычисляется: 

body(weighted.quantile)[[10]]
#> Fx <- cumsum(w)/sum(w)

Вы можете видеть, что он не начинается с 0. В вашем случае, первый элемент будет 0.3333.

Итак, чтобы показать, что это случай, давайте напишем над этой строкой с правильным выражением. (Сначала нам нужно разблокировать привязку, чтобы получить доступ к функции) 

unlockBinding("weighted.quantile", asNamespace("spatstat"))
body(weighted.quantile)[[10]] <- substitute(Fx <- (cumsum(w) - min(w))/(sum(w) - min(w)))

Теперь мы получаем правильный результат для взвешенных квантилей (включая правильную медиану) 

weighted.quantile(x, w)
#>   0%  25%  50%  75% 100% 
#> 10.0 10.5 11.0 11.5 12.0
...