Бутстрап, или прикладная статистика почти без формул

BootstrapВ институтах студентов учат интегрировать аналитически, а потом обнаруживается, что на практике интегралы почти все считают численными методами. Ну или по крайней мере проверяют таким образом аналитическое решение.

В статистике тоже есть нечестный метод, который позволяет получить примерный ответ на многие практические вопросы без анализа, грубой компьютерной силой: бутстрап (англ. bootstrap). Придумал и опубликовал его в 1979 году Брэдли Эфрон.

Допустим, есть у нас интернет-магазин, где мы торгуем всякой всячиной и привлекаем клиентов разными способами. Понятное дело, что мы постоянно что-то тестируем — расположение картинок и кнопок на странице, рекламный текст, AdWords, баннеры на сайтах партнёров и так далее. И вот свежие результаты — в тестовой группе из 893 пришедших у нас что-то купили 34, а в контрольной группе из 923 пришедших что-то купили 28.

Возникает вопрос — идти к начальству и говорить «в тестовой группе конверсия 3.81%, в контрольной группе 3.03%, налицо улучшение на 26%, где моя премия?» или продолжать сбор данных, потому что разница в 6 человек — ещё не статистика?

Эту задачу несложно решить аналитически. Видим две случайные величины (процент конверсии в тестовой и контрольной группе). При большом количестве наблюдений биномиальное распределение похоже на нормальное. Нас интересует разность. Нормальное распределение бесконечно делимо, вычитаем матожидания и складываем дисперсии, получаем матожидание 34/893-28/923 = 0.77% и дисперсию (34/893)*(1-34/893)/893+(28/923)*(1-28/923)/923. Стандартное отклонение равно корню из дисперсии, в нашем случае 0.85%. Истинное значение с 95% вероятностью лежит в пределах плюс-минус двух стандартных отклонений от матожидания, то есть между -0.93% и 2.48%.

Так что премия пока не светит, надо продолжать собирать данные.

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

Алгоритм:


  1. Выбираем наугад одно наблюдение из имеющихся.
  2. Повторяем пункт 1 столько раз, сколько у нас есть наблюдений. При этом некоторые из них мы выберем несколько раз, некотороые не выберем вообще — это нормально.
  3. Считаем интересующие нас метрики по этой новой выборке. Запоминаем результат.
  4. Повторяем пункты 1-3 много раз. Например, 10 тысяч. Можно меньше, но точность будет хуже. Можно больше, но долго будет считать.


Теперь у нас есть распределение, на которое мы можем посмотреть или что-то по нему посчитать. Например, доверительный интервал, медиану или стандартное отклонение.

Обратите внимание на то, что мы не делаем никаких предположений о распределении чего-либо. Распределения могут быть несимметричными, с толстыми хвостами и вообще иметь причудливую форму. Алгоритм от этого не меняется.

Чудес, правда, не бывает. Если у распределения нет матожидания (такие встречаются), бутстрап его не найдёт. Ну то есть он найдёт матожидание выборки, но не генеральной совокупности. То же касается ситуации, когда выборка нерепрезентативна или просто маленькая.

Бутстрап прост в имплементации. Приведённый ниже пример написан на C (проще не бывает) и кроме ввода-вывода использует только две библиотечные функции: генератор псевдослучайных чисел и сортировку. Разберём его по частям.

Пролог особых комментариев не требует. В C нет bool, данные будем хранить в int.

#include <stdio.h>
#include <stdlib.h>

typedef int Data_t;
#define ARRAY_SIZE(x) sizeof(x)/sizeof(x[0])


Функция, которая на лету делает выборку и считает по ней процент конверсии.
Правильнее было бы использовать более точный алгоритм вычисления среднего значения,
но для нашего примера это не важно.

static double bootstrap(const Data_t* data, unsigned n)
{
    unsigned i;
    double sum = 0;
    for (i = 0; i < n; i++) {
        sum += data[rand() % n];
    }
    return sum / n;
}


Функция сравнения для сортировки результатов

static int compare(const void* a, const void* b)
{
    if (*(double*)a > *(double*)b) return 1;
    if (*(double*)a < *(double*)b) return -1;
    return 0;
}


Исходные данные

int main(int argc, char* argv[])
{
    Data_t test[893] = { 0 };
    Data_t control[923] = { 0 };
    unsigned i;
    for (i = 0; i < 34; i++) {
        test[i] = 1;
    }
    for (i = 0; i < 28; i++) {
        control[i] = 1;
    }


Инициализируем генератор псевдослучайных чисел параметром из командной строки. Если мы
всё сделали правильно, то результаты не должны сильно плавать при изменении этого
параметра.

    if (argc == 2) {
        srand(atoi(argv[1]));
    }


Сюда будем складывать результаты

    double t_minus_c[10000];


Главный цикл

    for (i = 0; i < ARRAY_SIZE(t_minus_c); i++) {
        t_minus_c[i] = bootstrap(test, ARRAY_SIZE(test)) 
                     - bootstrap(control, ARRAY_SIZE(control));
    }


Определяем 95% доверительный интервал: cортируем результаты, отбрасываем 2.5% снизу и столько же сверху, показываем результат.

    qsort(t_minus_c, ARRAY_SIZE(t_minus_c), sizeof(double), compare);
    printf("LCL=%g%%\n", 100. * t_minus_c[250]);
    printf("UCL=%g%%\n", 100. * t_minus_c[9750]);
    return 0;
}


Проверяем:

$ gcc -Wall -o bootstrap bootstrap.c
$ ./bootstrap
LCL=-0.891368%
UCL=2.43898%


Ещё пару раз с другими псевдослучайными числами:

$ ./bootstrap 42
LCL=-0.86589%
UCL=2.43898%
$ ./bootstrap 2013
LCL=-0.959673%
UCL=2.52184%


Похоже на теоретический результат (от -0.93% до 2.48%).

И зачем было огород городить?


У этой задачи есть простое аналитическое решение, но у многих реальных задач его или нет вообще, или оно есть, но очень сложное. Представьте, что вместо процента конверсии нас интересует отношение прибыли от клиента к затратам на его привлечение. Распределение такой метрики вряд ли будет нормальным, и формулы перестанут укладываться в пару строчек. А бутстрап будет работать точно так же, достаточно поменять Data_t на double и положить туда новые данные.

Первоисточники


  1. Efron B (1979). Bootstrap methods: Another look at the jackknife. Ann. Statist. 7 1–26
  2. Donald E. Knuth. Seminumerical Algorithms, volume 2 of The Art of Computer Programming, chapter 4.2.2, page 232. Addison-Wesley, Boston, third edition, 1998.
Поделиться публикацией

Комментарии 23

    +7
    Это тот редкий случай, когда обосновать законность применения метода сложнее, чем применить сам метод. И его самый важный недостаток, как и у всех методов Монте-Карло, точность растет медленнее, чем количество наблюдений.
      +1
      То же самое можно сказать про классические параметрические методы. Я ни разу не видел в настоящих данных нормальное распределение, поэтому строго говоря всякие t-тесты и p-значения от них неприменимы (что мало кому мешает их использовать на практике). Точность у параметрических методов тоже растёт медленнее, чем количество наблюдений.

      Другое дело что если предположения, лежащие в основе параметрических методов, выполняются, тогда их результат в некотором смысле оптимален. Про Монте-Карло такое сказать нельзя, по крайней мере при конечном количестве итераций. А при бесконечном количестве итераций получается сферический конь в вакууме.
        +1
        1. Причиной того, что Вы не видели в настоящих данных нормального распределения может быть то, что Вы ожидаете теоретическое нормальное распределение (вместо емпирического). Посмотрите на гистограммы выборок, полученных генератором случайных нормально распределенных чисел при разных размерах выборок. Вы должны сравнивать именно с такими гистограммами, а не с теоретической кривой. При малых размерах выборки гистограммы нормального распределения достаточно сильно отличаются.
        2. Не согласен что параметрические тесты неприменимы: при умеренных отклонениях от нормального распределения они всегда дают правильный ответ. Но согласен, что их часто используют в ситуациях когда они становятся ошибочны.
        3. Понятия точности для статистического теста не определено. Есть понятие мощности теста. Мощность — вероятность, что тест отклонит нулевую гипотезу, если она на самом деле не верна. Мощность параметрических тестов гораздо выше чем непараметрических и возрастает (в приближении) пропорционально квадрату количества наблюдений.
        4. Число возможных сочетаний с повторениями конечно. При использовании в бутстрапе всех таких сочетаний результат оптимален. Конечно, в примере, данном в статье, число сочетаний с повторениями слишком огромно, но при малых размерах выборки вполне достижимо.
          0
          Извиняюсь, в пункте 3 должно быть "… пропорционально квадратному корню с количества наблюдений"
            +3
            1. Нет, дело не в шуме из-за размера выборки (у меня они большие), дело именно в кривизне распределений. Как в правиле 20/80, по которому 20% людей выпивают 80% пива. Некоторые распределения можно приблизить логнормальным (количество денег в банке) или всякой экзотикой вроде распределения Васичека (потери), но нет гарантии, что это правильно отражает суть того, что на самом деле происходит. Даже такие вроде бы нормальные величины как вес или рост не следуют нормальному распределению, вопрос только в том, какой должен быть размер выборки, чтобы это понять.

            2. Параметрические тесты неприменимы с математической точки зрения. Они действительно дают правильные ответы даже тогда, когда условия строго говоря нарушены, так же как бутстрап даёт правильные ответы в случаях, для которых его оптимальность не доказана.

            3. Мощность параметрических тестов гораздо выше, но они опираются на предположения, которые не обязательно выполняются. При больших размерах выборки (10^5 и больше) t-тест почти всегда находит статистически значимую разницу между тестовой и контрольной группой в A/A тесте, просто потому что они выбраны хоть немного неслучайно.
              +3
              Поддерживаю и согласен со всем. 1. В разных сферах — разная природа данных. Их природа и определяет форму распределения. Но следует учитывать, что в статистических тестах нормальным должно быть распределение остатков, а не исходной выборки (хоть они и связаны). У меня при анализе медицинских данных около 30% случаев имеют распределение остатков достаточно близкое к нормальному.
              3. Поэтому размер эффекта более важен нежели р-величины.
            +1
            Хм, я вот регулярно сталкиваюсь с нормальным распределением в настоящих данных (психологические измерения), но тут всё зависит от того, как мы определяем нормальность. Критерии Шапиро-Уилкса, Колмогорова-Смирнова и «я на глазок прикинул, мамой клянусь — гауссиана» — дают порой довольно разные результаты, тем более если размеры выборок невелики.

            Когда я преподавал матметоды я обычно просил студентов проверить нормальность равномерного распределения из 30 случаев (на 10, кажется, интервалах). Ну так вот, нет оснований считать, что такая выборка извлечена из совокупности с ненормальным распределением — если мы пользуемся формальным критерием.

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

              Но тут даже не в тестах дело, а просто в сути данных. Вот допустим меряю я атмосферное давление, отвалился у меня датчик и три месяца в данных одни нули, пока датчик обратно не привинтил. Вероятность нулей три месяца подряд по нормальному распределению очень низкая, а в реальных данных они есть.

              Или то же потребление пива. Много людей не пьют вообще, у пьющих есть горб где-то посередине, а потом — длинный тонкий хвост. Опять ненормальное распределение. Можно нули отдельно моделировать, а пьющих преобразовать, но как именно? Предположения какие-то надо делать по поводу формы хвоста, от них результат зависеть будет.

              Одна из особенностей бутстрапа как раз в том, что он не требует делать такие предположения.
                +1
                Действительно, тесты на нормальность при малых размерах выборки не смогут отбросить гипотезу о нормальном распределении даже если ним там и не пахнет, а при больших размерах даже при пустяковых отклонениях дадут р<0,05. Поэтому их использование справедливо только для небольшого интервала размеров (где то посередине). Так вот, на глаз — а именно инспекция квантиль-квантильного нормального графика — как раз самое то. А нули при измерении давления — это выбросы, от них нужно избавлятся.
                  0
                  Избавляться от выбросов чревато. Недавно вот финансовый кризис был, не в последнюю очередь потому, что финансовые моделисты наделали предположений о форме распределений и избавились от выбросов. Нассим Талеб про это хорошо пишет.

                  Одно дело если я метеорологией занимаюсь, тогда спокойно могу нули выкинуть. Если атмосферное давление действительно равно нулю, то мои исследования никому уже не нужны.

                  А если я инженер, который датчик давления делает, для меня эти выбросы — весьма важные данные.
                  0
                  Тут проблема-то в том, что если мы в принципе оцениваем «среднее», хоть бутстрапом, хоть как угодно, мы используем этот параметр для оценки данных, а вот осмысленность его использования — зависит от распределения. Поэтому всё равно, чтобы сравнивать две выборки по какому-то такому вот параметру, нужно, чтобы их распределения относились к _одному типу_. Не обязательно — нормальному.

                  Про пьющих пиво — это распределение Парето. Меня вообще очень интересует вопрос применимости различных статистических техник к этому распределению, так как я занимаюсь вещами на грани с лингвистикой, где это распределение является куда более естественным, чем нормальное (см. закон Ципфа).
                    0
                    Согласен. А для регрессии с Парето-распределенной зависимой переменной можно использовать достаточно универсальный инструмент — GAMLSS (имплементирован на R и Java).
            +2
            Берем на вооружение!
              +2
              В астрономии мы опыты повторять не можем, а статистику навести хочется, потому пользуем бутстрап.
                +2
                В C нет bool, данные будем хранить в int.

                Не принципиально, но вообще-то в C есть тип bool, начиная с C99.
                  +1
                  Благодаря КДПВ узнал, что bootstrap означает «шнурок», как-то раньше в голову не приходило. Спасибо!
                    +3
                    bootstrap — это не шнурок, а петля на задней стороне ботинок.
                      0
                      Да, действительно. Не разобрался.
                      0
                      В данном контексте — вытягивание себя из болота за шнурки, по примеру барона Мюнхаузена.
                      +2
                      По поводу реализации: стандартный библиотечный rand() может иметь весьма короткий период повторения последовательности: в лучшем случае около 10^9, а в худшем — 32767! Для любых реализаций методов Монте-Карло рекомендуется использовать известные алгоритмы генерации псевдослучайных чисел, такие как Mersenne-Twister.
                        0
                        В примере rand() вызывается около 2*10^7 раз. У стандартного библиотечного rand() из glibc период повторения последовательности 2*10^9. Там вроде сдвиговый регистр с обратной связью используется, отсюда и цикл длиной 2^31.

                        То есть в это ограничение мы пока не упёрлись, но на настоящих данных если размер выборки окажется кратен длине цикла, то действительно можно попасть.
                          +1
                          В MS CRT период 32767.
                        0
                        В качестве тренировки реализовал то же самое на Haskell:

                        -- bootstrap.hs
                        module Main where
                        
                        import Data.List
                        import Control.Monad
                        import Control.Monad.Trans.State
                        import System.Environment
                        import System.Random
                        import Text.Printf
                        
                        randomsN :: (RandomGen g, Random a) => (g -> (a, g)) -> Int -> g -> ([a], g)
                        randomsN rnd n = runState (replicateM n (state rnd))
                        
                        bootstrap :: (RandomGen g) => (Int, Int) -> g -> (Double, g)
                        bootstrap (sN, tN) gen =
                        	(successCount / realToFrac tN, newGen)
                        	where 
                        		rnd = randomR (1, tN)
                        		(rs, newGen) = randomsN rnd tN gen
                        		successCount = genericLength $ filter (<=sN) rs
                        	
                        minusMetric :: RandomGen g => (Int, Int) -> (Int, Int) -> g -> (Double, g)
                        minusMetric test control g =
                        	(m1 - m2, g2)
                        	where
                        		(m1,g1) = bootstrap test g
                        		(m2,g2) = bootstrap control g1
                        
                        createGen :: [String] -> IO StdGen
                        createGen [] = getStdGen
                        createGen [seed]= return $ mkStdGen $ read seed
                        
                        main :: IO ()
                        main = do
                        	gen <- getArgs >>= createGen
                        	let 
                        		metric = minusMetric (34, 893) (28, 923)
                        		(ms,_) = randomsN metric 10000 gen
                        		os = sort ms
                        		lcl = (os !! 250) * 100
                        		ucl = (os !! 9750) * 100
                        	putStrLn $ printf "LCL=%g%%" lcl
                        	putStrLn $ printf "UCL=%g%%" ucl
                        


                        Компилируем:
                        $ ghc --make -O2 bootstrap.hs
                        


                        Запускаем:
                        $ ./bootstrap.exe
                        LCL=-0.8913676736965855%
                        UCL=2.4389770442796324%
                        
                        $ ./bootstrap.exe 42
                        LCL=-0.8840882316900799%
                        UCL=2.5182016381170995%
                        
                        $ ./bootstrap.exe 2013
                        LCL=-0.8768087896835746%
                        UCL=2.4535359282926428%
                        
                        $ ./bootstrap.exe 2013
                        LCL=-0.8768087896835746%
                        UCL=2.4535359282926428%

                        Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                        Самое читаемое