Публикация носит характер описания тропинки, выводящей к эффективной алгоритмизации методов вычисления доверительного интервала (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), "]")

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


В чем еще существенный плюс таких приближенных вычислений? Да можно просто перенести весь этот код к данным (которых, на самом деле, на несколько порядков больше), живущих в Clickhouse. Такие алгоритмы перекладываются на SQL в два клика. При этом счет времени получения расчетных показателей в продуктиве пойдет на миллисекунды.


P.S.


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

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