Публикация носит характер описания тропинки, выводящей к эффективной алгоритмизации методов вычисления доверительного интервала (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)]Код взят отсюда.
Про применение биномиального распределения:
- J.H. Zar, Biostatistical Analysis, Fifth edition 2010, PEARSON, ISBN 10: 0-13-100846-3. "24.9 CONFIDENCE INTERVAL FOR A POPULATION MEDIAN" (pages 548-549)
- J.H. Zar, Biostatistical Analysis, Fifth edition 2014, PEARSON, ISBN 10: 1-292-02404-6. "9 CONFIDENCE INTERVAL FOR A POPULATION MEDIAN" (pages 584-585)
CI для разницы медиан
Проверяем гипотезу о статистической неразличимости медиан двух различных выборок.
Отправные дискуссии на StackExchange "How to construct a 95% confidence interval of the difference between medians?" и "Bootstrap hypothesis test for median of differences".
В последней, хоть вопрос шел об одном, но ответ приложен на нужный вопрос.
Плюс еще 2 публикации:
- Price R.M., Bonett D.G. "Distribution-free confidence intervals for difference and ratio of medians". J Stat Comput Simul 72, 2002
- Price R.M., Bonett D.G. "Statistical inference for a linear function of medians: Confidence intervals, hypothesis testing, and sample size requirements
# 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), "]")Код взят отсюда.
В чем еще существенный плюс таких приближенных вычислений? Да можно просто перенести весь этот код к данным (которых, на самом деле, на несколько порядков больше), живущих в Clickhouse. Такие алгоритмы перекладываются на SQL в два клика. При этом счет времени получения расчетных показателей в продуктиве пойдет на миллисекунды.
P.S.
- Публикация никак не претендует на какую-либо полноту или глубину. Просто подход к решению практической задачи. Если кто-то знает иные быстрые и элегантные способы расчета, буду благодарен комментариям.
- Простые и внятные рекомендации на заданные вопросы очень трудно найти в рафинированном виде. Много встречается общих рассуждений, либо ошибочных утверждений. Интересные русскоязычные изыскания на примерно подобные вопросы с небольшим обзором "литературы" были найдены только здесь: "ДИ медианы".
Предыдущая публикация — «Карантин, онлайн-системы и data science. Кто думает об удержании клиентов?».
