c-i-median.R

---
title: "ДИ медианы"
author: "mz"
date: '`r Sys.time()`'
output:
  html_document: default
  pdf_document:
    keep_tex: yes
header-includes: \usepackage[T1,T2A]{fontenc}
---

```{r echo=FALSE,results='hide'}
suppressMessages(library(ggplot2))
suppressMessages(library(magrittr))
suppressMessages(library(dplyr))
suppressMessages(library(reshape2))
suppressMessages(library(knitr))
suppressMessages(library(nortest))
```

```{r setup, include=FALSE, results='hide'}
set.seed(2718281)
quant <- c(.025,.975)
# Это для моделирование выборок
# на основе которых будут рассчитываться ДИ
# различными способами
s.size          <- 1500
s.repetitions   <- 1250

# Это для бутстрепинга
bs.s.size          <- 1100
bs.s.repetitions   <- 100

gen.majority    <- rchisq(100000,8)
gen.median      <- median(gen.majority)
```

# Жили-были…

> «Космос велик — сообщает путеводитель по галактике во вступлении, страшно 
велик, вы просто не поверите на сколько умопомрачительно он велик. и т. д. 
Еще путеводитель сообщает, что набрав полную грудь воздуха, вы можете прожить в 
вакууме около 30 секунд. Однако, вероятность того, что в умопомрачительных 
просторах космоса вас подберет в течении 30 секунд проходящий мимо корабль, 
равняется одному к двум в степени 2079460347. По удивительному совпадению 
таков номер телефона некой квартиры в Лондоне, где Áртур познакомился на 
маскараде с очаровательной девушкой, которой не смог понравится.»
> --- <cite>[Д. Адамс «Автостопом по галактике» 0:25:40][adams1]</cite>

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

Для ясности освежим определение доверительного интервала, а называть им принято 
такой диапазон значений, построеный по данным единичной случайной выборки 
известного объёма, что при бесконечно повторении выборочного эксперимента 
процент случаев попадания значения исследуемого параметра (в нашем случае 
медианы) генеральной совокупности в рассчитаный интервал буден не менее 
заданной доли (обычно `r (quant[2]-quant[1])`) повторений.

Сразу оговоримся, что интерквáртильное расстояние включает 50 процентов 
наблюдений *выборки* и не соотносится с понятием доверительного интевала, 
как иногда утверждают.

Легкий гуглеж ведет на стандартные нынче 
[stackowerflow](http://stats.stackexchange.com), 
[старый форум по R](https://stat.ethz.ch), какой-то неведомый ресурс по 
[SPSS-тулам-скриптам](http://spsstools.net), 
[IBM'овская база знаний](http://www-01.ibm.com/support) и, что оказалось 
особенно полезным, на java-код одного околонаученого проекта на 
[github.com](http://github.com) с отличным комментариями.

Подходы к подсчету:

  * Наивное программное сэмплирование. Вычислительно емкий, но простой в 
  реализации и отладке метод, доступный в большинстве спец. пакетов, да и в 
  Excell'е даже. Будем использовать как арбитражный метод, потому как он 
  больше всего похож на само определение доверительного 
  интервала [[4][4],[5][5],[6][6]];
  * Чудовищный метод ~~с аналь~~ с «магическими» числами-константами 
  и SPSS'ом, являющийся репликой из какого-то [довольно древнего учебника][7], 
  который уже даже и скачать нельзя, но можно купить у шпрингера за 90 долларов 
  США, интересен только как история [[3][3],[8][8]];
  * Относительно аналитический метод через квáнтили биномиального распределения, 
  вот его-то и будем тестировать [[2][2],[5][5],[6][6]].

# [Das Model][3]

> «Каждое слово Лизы картинку меняет,<br/>
> Событий много, мозг не успевает»
>  --- <cite>[Ка4 «Лиза Овчинникова»][ka41]</cite>

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

## Генеральная совокупность

Для моделирования выборок нам понадобится некоторая генеральная совокупность, 
для наглядности в части асимметрии будет использована выборка объёмом $10^5$ 
из распределения $\chi^2$ с небольшим (8) числом степеней свободы. График 
плотности распределения не оставляет сомнений в его асиметричности:

```{r, echo=FALSE}
  ggplot(data.frame(x=gen.majority),aes(x))+geom_density()
```

Основные параметры полученной генеральной совокупности можно обобщить 
в следующей таблице:

```{r, echo=FALSE}
  kable(as.matrix(summary(gen.majority)))
```

## Выборочный метод

```{r echo=FALSE}
  sampling.start  <- Sys.time()
  samples         <- sapply(1:s.repetitions, 
                            function(i){sort(sample(gen.majority,s.size))})
  # для рисования ggplot'ом нужно перевести матрицу в таблиуц
  # а для правильно группировки нужно, чтобы melt давал правильные имена группам
  # проще всего этого достичь, если сразу поименовать строки и столбцы 
  # их номерами
  dimnames(samples) <- list(1:s.size,1:s.repetitions)
  melted.samples    <- melt(samples,varnames = c('sample.num','sample.serie'))
  sampling.duration <- as.numeric(Sys.time() - sampling.start) %>% 
    round(2)
  
  samples.plotting.start <- Sys.time()
  samples.dens.plot <- melted.samples %>% 
    ggplot(aes(value,group=sample.serie)) + geom_density()
  samples.plotting.duration <- as.numeric(Sys.time() - samples.plotting.start) %>% 
    round(2)
```

Сгенерируем `r s.repetitions` случайных выборок по `r s.size` значений от 
условной генеральной совокупности из предыдущего раздела. Это занимает 
приблизительно `r sampling.duration` секунды. И за еще `r samples.plotting.duration` 
можно построить график всех диаграмм плотностей распределения для полученных 
выборок.

```{r, echo=FALSE}
  samples.dens.plot
```

## Наивный подсчет через ресэмплинг

```{r echo=FALSE}
  sampled.medians <- melted.samples %>% 
    group_by(sample.serie) %>% 
    summarise(median = median(value)) %>% 
    ungroup() %>% 
    arrange(median) %>% 
    select(median) 
  medians.dens.plot <- sampled.medians %>% 
    ggplot(aes(median))+geom_density()  
  
  # Бутстрепинг-процедура
  get.median.measure <- function(sample.src){
    sapply(1:bs.s.repetitions, 
           function(i){sample(sample.src,
                              size = bs.s.size,
                              replace = T)}) %>% 
      apply(2, median) %>% 
      sort()
  }
  
  bs.start <- Sys.time()
  bs.medians <- apply(samples, 2, get.median.measure)
  bs.duration  <- Sys.time() - bs.start
  
  dimnames(bs.medians) <- list(1:bs.s.repetitions,1:s.repetitions)
  bs.medians.plot.start <- Sys.time()
  bs.medians.plot <- bs.medians %>% 
    t %>% 
    melt(varnames = c('sample.num','sample.serie')) %>% 
    ggplot(aes(value,group=sample.serie)) + geom_density()
  bs.medians.plot.duration <- Sys.time() - bs.medians.plot.start
  
  c.i.s <- bs.medians %>% 
    apply(2,
          function(z){quantile(z,quant)})
  miss.c.i <- c.i.s %>%   
    apply(2,
          function(y){y[1] > gen.median  | y[2] < gen.median })
```

Бóльшая часть работы уже выполнена, осталось только набутстрепать нужное
количество и подсчитать разброс `r (quant[2]-quant[1])*100`% медиан для всех 
выборок. Занимает это `r round(bs.duration,2)` c. и 
еще `r round(bs.medians.plot.duration,2)` с. занимает построение графика 
плотностей вероятностей распределений медиан набутстрепаных в количестве 
`r bs.s.repetitions` для каждой из `r s.repetitions` случайных выборок 
значений генеральной совокупности. 

```{r echo=FALSE}
  bs.medians.plot
```

Плотность распределения медиан выглядит следующим образом:

```{r echo=FALSE}
  medians.dens.plot
```

```{r echo=FALSE}
lf <- lillie.test(sampled.medians$median)
```

На вид нормально, тест Лиллиефорса (Шапиро-тест не может работать с выборками 
мощнее 5000) говорит: p-value=`r round(lf$p.value,3)`, то есть оно не мене
нормально, чем нормальное распределение.

Границы квáнтилей `r (quant[2]-quant[1])*100`%-го интервала (это не ДИ!) равны: 

```{r echo=FALSE}
  quantile(sampled.medians$median,quant)
```

Что касается интервалов, получаемых бутстрепингом по выборкам, то они содержат 
в себе медиану генеральной совокупности 
в `r round(100-sum(miss.c.i)/length(miss.c.i)*100,2)` процентах случаев.

## Аналитические варианты

Этот вариант вычислительно значительно легче, идея проста: нужно получить 
квáнтили биномиального распределения с параметром количества повторов равным 
числу наблюдений (в нашем случае — `r s.size`) и параметров вероятности 0.5.

```{r echo=FALSE}

  bounds    <- qbinom(c(.025,.975), s.size, 0.5)
  c.i.s  <- samples %>% 
    apply(2, function(x){x[bounds]}) 
  miss.c.i <- c.i.s %>% 
    apply(2, function(y){y[1] > gen.median  | y[2] < gen.median})

```

Всего несколько секунд, но результат по качетсву не отличается, процент 
попадания медианы генеральной совокупности в построенный интервал — 
`r round(100-sum(miss.c.i)/length(miss.c.i)*100,2)` процента.

Есть еще вариант из статьи [Д. Олива][5], в которой он показывает, что вполне
можно построить доверительный интервал к медиане методом похожим на используемы
при построении ДИ к среднему, то есть через квáнтили t-распределения (того 
которое Стьюдента, который на самом деле У. Госсет).

Содержание метода просто:

  1. Рассчитать выборочную медиану ($M_{S}$);
  1. Рассчитать границы индексов в упорядоченной выборке:
      
      * Нижнюю как $n_{L} = \lfloor{n \over 2}\rfloor - \lceil{\sqrt{n \over 4}}\rceil$;
      * Верхнюю как $n_{U} = n - n_{L}$;
      
  1. Рассчитать стандартную ошибку медианы как $SE_{M} = .5 \times (S_{n_{U}} - S_{n_{L}+1})$
  1. Взять квáнтили t-распределения с числом степеней свободы на единицу меньшим
     чем количество элементов между вычисленными ранее границами — $t_{n_{U} - n_{L} - 1}^{.025,.975}$;
  1. Вычесть из выборочной медианы произведение стандартной ошибки на величину 
     нижнего квáнтиля для получения нижней границы ДИ, и сложить медиану с 
     произведением верхнего квáнтиля на стандартную ошибку для получения верхней.

```{r }
  olive.median.ci <- function(samples){
    s.samples <- sort(samples)
    n         <- length(s.samples)
    m         <- median(s.samples)
    nL <- floor(n/2) - ceiling(sqrt(n/4))
    nU <- n - nL
    SE <- .5 * (s.samples[nU] - s.samples[nL+1])
    return(m + SE * qt(.975,nU-nL-1)*c(-1,1))
  }

  c.i.s <- samples %>% 
    apply(2, olive.median.ci)
  
  miss.c.i <- c.i.s %>% 
    apply(2, function(y){y[1] > gen.median  | y[2] < gen.median})
  miss.rate <- 100-sum(miss.c.i)/length(miss.c.i)*100
```

Вычислнные интервалы накрыл «истинную» медиану в `r round(miss.rate,2)` 
процентах случаев.

Так же в R есть способ использовать процедуру теста Уилкоксона, но её 
интерпретация мне, да и [не только мне][9], остается недоступной. Тем не 
менее просто для справки приведу пример:

```{r}
  wilcox.test(samples[,1],conf.level=0.95, conf.int = TRUE)
```

```{r echo=FALSE}
  wt.start <- Sys.time()
  c.i.s  <- samples %>% 
  apply(2,
        function(x){wilcox.test(x,
                                conf.level=0.95, 
                                conf.int = TRUE)$conf.int})
  
  miss.c.i <- c.i.s %>% 
    apply(2, function(y){y[1] > gen.median  | y[2] < gen.median})
  miss.rate <- 100-sum(miss.c.i)/length(miss.c.i)*100
  wt.duration <- Sys.time() - wt.start
```

`r round(Sys.time() - wt.start,2)` секунды и `r round(miss.rate,2)` процента 
«попаданий» в медиану. Что требует дополнительного исследования. А в остальном 
Дизраэли бы мог гордиться.

Воспроизвести исследование можно из исходных кодов на 
[gist](https://gist.github.com/aa989190f363e46d/3d9a2380cbec14738c48398d814d976f).

[adams1]:http://rutracker.org/forum/viewtopic.php?t=4059570
[ka41]:https://youtu.be/AenPw0AvpK8?t=184
[1]:https://www.youtube.com/watch?v=OQIYEPe6DWY
[2]:https://github.com/magro/elliptic-benchmark/blob/master/src/main/java/bb/science/Bootstrap.java
[3]:http://spsstools.net/ru/syntax/234/
[4]:http://www-01.ibm.com/support/docview.wss?uid=swg21476966
[5]:http://exploringdatablog.blogspot.com.by/2012/04/david-olives-median-confidence-interval.html
[6]:http://stats.stackexchange.com/questions/21103/confidence-interval-for-median
[7]:https://www.amazon.com/Mathematics-Handbook-Science-Engineering-Lennart/dp/3540211411
[8]:folk.ntnu.no/slyderse/medstat/KLMED8004/Ikkeparamet.eng.11.09.pdf
[9]:https://stat.ethz.ch/pipermail/r-help/2011-May/277910.html