Медианы выборок. Доверительные интервалы и сравнение

от автора

Публикация носит характер описания тропинки, выводящей к эффективной алгоритмизации методов вычисления доверительного интервала (Confidence Interval = CI) для

  • медианы распределения;
  • разницы медиан двух распределений.

Задача сугубо практическая, в глубины математики погружаться можно, но это не самоцель, да и не всегда хватает баллона, чтобы добраться до дна.

Выборки по объему большие, 10^5 — 10^7 записей, ощутимо ассимметричные, с длинными хвостами, могут иметь несколько мод. В этом случае медианы более устойчивы к выбросам.

Применение классической статистики, например, критерия Уилкоксона-Манна-Уитни, для оценки разницы медиан на таких объемах не проходит. Да и очень много чего читать надо под звездочками, чтобы правильно применять эти критерии. Ведь этот критерий проверяет отнюдь не равенство медиан, да и для медиан он работает только при одинаковых формах распределений двух выборок. И т.д. и т.п.

Хвататься за молоток бутстрапа можно, но и с ним надо думать + на симуляцию требуется время и память.

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

CI для медиан

Отправная дискуссия на StackExchange "Confidence interval for median" выводит на статью David J. Olive "A Simple Confidence Interval for the Median", 2005 и весьма элегантный код, проще которого сложно что-то придумать:

x <- faithful$waiting sort(x)[qbinom(c(.025,.975), length(x), 0.5)]

Код взят отсюда.

CI для разницы медиан

Проверяем гипотезу о статистической неразличимости медиан двух различных выборок.
Отправные дискуссии на StackExchange "How to construct a 95% confidence interval of the difference between medians?" и
"Bootstrap hypothesis test for median of differences".

В последней, хоть вопрос шел об одном, но ответ приложен на нужный вопрос.

И еще 2 публикации:

Исходный код

# test from table 3 of b&p 2002 x1 <- c(77, 87, 88, 114, 151, 210, 219, 246, 253, 262, 296, 299, 306,         376, 428, 515, 666, 1310, 2611) x2 <- c(59, 106, 174, 207, 219, 237, 313, 365, 458, 497, 515, 529,         557, 615, 625, 645, 973, 1065, 3215)  # sort vectors x1 <- sort(x1) x2 <- sort(x2)  # get medians x1_mdn <- median(x1) x2_mdn <- median(x2)  # stuff to calculate variance of medians x1_n <- length(x1) x2_n <- length(x2)  x1_aj <- round((x1_n + 1) / 2 - x1_n ^ (1 / 2)) x2_aj <- round((x2_n + 1) / 2 - x2_n ^ (1 / 2))  z <- 1.855 # from table 1 of b&p 2002, see p. 376  # calculate variance x1_var <- ((x1[x1_n - x1_aj + 1] - x1[x1_aj]) / (2 * z)) ^ 2 x2_var <- ((x2[x2_n - x2_aj + 1] - x2[x2_aj]) / (2 * z)) ^ 2  # contrast coefficients, such that its median(d) - median(dg) x1_cj <- 1 x2_cj <- -1  # median difference mdn_diff <- x1_mdn * x1_cj + x2_mdn * x2_cj  # standard error mdn_diff_se <- (((x1_cj ^ 2) * x1_var) + ((x2_cj ^ 2) * x2_var)) ^ (1 / 2)  # 95% confidence interval lb <- mdn_diff - 1.96 * mdn_diff_se ub <- mdn_diff + 1.96 * mdn_diff_se  # within roundng error of p. 376 of b&p 2002 paste0(mdn_diff, " [", round(lb), ", ", round(ub), "]")

Код взят отсюда.

P.S.

  1. Публикация никак не претендует на какую-либо полноту или глубину. Просто подход к решению практической задачи. Если кто-то знает иные быстрые и элегантные способы расчета, буду благодарен комментариям.
  2. Интересные изыскания на примерно подобные вопросы были найдены здесь: "ДИ медианы".

Предыдущая публикация — «Карантин, онлайн-системы и data science. Кто думает об удержании клиентов?».

ссылка на оригинал статьи https://habr.com/ru/post/505834/


Комментарии

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *