Недавно меня попросили написать отзыв на автореферат кандидатской диссертации, в которой обсуждалось моделирование случайных величин с использованием Python и C++. Я разбираюсь в моделировании, но не в программировании. Обсуждая работу, я поинтересовался у соискателя, почему он выбрал эти инструменты и не рассматривал ли Excel. Он ответил, что в их среде Excel не используется. «А жаль», — подумал я. Особенно учитывая, что в работе выборки не превышали сотни элементов. Excel легко справляется даже с миллионом и имеет десятки встроенных функций для таких целей. 

В этой статье в блоге ЛАНИТ я покажу, как с помощью Excel можно эффективно генерировать случайные величины различных распределений и почему этот инструмент не стоит недооценивать. 

С появлением в Excel 2019 новых функций массива генерация случайных чисел существенно упростилась. В Excel также есть десятки статистических функций, возвращающих значения разнообразных распределений, а современные ПК быстро создадут выборку даже в миллион значений. Так почему же не использовать этот широко распространенный инструмент?

В основе моделирования любых распределений лежит равномерное распределение. С него и начнем.

Равномерное распределение

Распределение вероятностей называют равномерным, если на интервале [a, b], которому принадлежат все возможные значения случайной величины X, плотность распределения сохраняет постоянное значение.

Плотность равномерного распределения:

Наиболее интересен случай, когда . Такое равномерное распределение называют стандартным (рис. 1а):

Интегрируя плотность вероятности (1.2), получаем функцию распределения (рис. 1б):

Процедура в Excel для каждого распределения будет включать:

• генерацию массива случайных чисел,

• подсчет случайных чисел, попадающих в карманы (они же диапазоны),

• построение кривой (или гистограммы) плотности вероятности,

• нахождение среднего значения и стандартного отклонения и сравнение с теоретическими значениями для верификации процесса генерации.

Чтобы разыграть (сгенерировать) непрерывную стандартную равномерно распределенную случайную величину в Excel используем формулу:

  

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

Для построения гистограммы нужно выбрать интервалы конечной ширины (рис. 2б). Вероятность попадания величины (в результате испытания) в интервал, принадлежащий интервалу , равна его длине:

Например, разбивая интервал на 10 равных частей мы получим вероятность случайной величины Х попасть в один из отрезков (рис. 2в). Если разбить интервал на 20 равных частей столбцы гистограммы выстроятся на высоте 5%. Вот почему диаграммы плотности вероятности иногда строят без указания чисел по оси ординат.

В приложенном архиве вы найдете Excel-файлы, использованные для моделирования и построения диаграмм. Для уменьшения объема файлов я сохранил их со значением .

Параметры распределения и статистики выборки

В статистике генеральную совокупность описывают параметрами, обозначаемыми греческими буквами: математическое ожидание , среднеквадратичное отклонение . Выборки описывают статистиками, обозначаемыми латинскими буквами: среднее арифметическое (или просто среднее), стандартное отклонение .

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

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

Матожидание равномерного распределения:

Матожидание стандартного равномерного распределения:

Дисперсия равномерного распределения:

Среднеквадратичное отклонение стандартного равномерного распределения:

Теоретические значения параметров генеральной совокупности и статистики выборки для стандартного равномерного распределения приведены в верхней правой части рис. 2. Мы будем проверять адекватность выборки, сравнивая статистики с параметрами генеральной совокупности.

Статистические функции Excel

В Excel в разделе статистических функций представлено два десятка функций непрерывных распределений. Большинство из них парные, но не все:

Некоторые распределения представлены большим числом функций. Например, СТЬЮДЕНТ.РАСП, СТЬЮДЕНТ.РАСП.2Х, СТЬЮДЕНТ.РАСП.ПХ. Для наших целей достаточно одного представителя от семейства. Рассмотрим структуру функций на примере пары НОРМ.СТ.РАСП и НОРМ.СТ.ОБР.

Синтаксис:

=НОРМ.СТ.РАСП(;ЛОЖЬ) — возвращает плотность вероятности стандартного нормального распределения (рис. 4а);

=НОРМ.СТ.РАСП(;ИСТИНА) — возвращает [интегральную] функцию стандартного нормального распределения (рис. 4б);

=НОРМ.СТ.ОБР — возвращает обратное значение стандартного нормального распределения (рис. 4в). Обратите внимание на аргумент функции — , вероятность.

Аналогично устроены все функции в таблице на рис. 3. Функции с суффиксом РАСП позволяют строить плотность вероятности и функцию распределения, выбирая соответствующий последний аргумент, а функции с суффиксом ОБР — находить значение распределения по его вероятности. Последние используются для разыгрывания случайной величины.

Алгоритм прост: записываем формулу на основе функции с суффиксом ОБР и вместо одного значения вероятности р подставляем случайный массив десятичных значений от 0 до 1 (рис. 5а).

В таблице на рис. 3 представлены 8 функций с суффиксом ОБР, которые можно использовать для генерирования случайных величин, что называется «из коробки»: F.ОБР, БЕТА.ОБР, ГАММА.ОБР, ЛОГНОРМ.ОБР, НОРМ.ОБР, НОРМ.СТ.ОБР и СТЬЮДЕНТ.ОБР. 

НОРМ.ОБР, НОРМ.СТ.ОБР я описал выше. 

Оставшиеся шесть функций:

Если в Excel нет готовой функции для генерирования распределения, можно попробовать найти обратную функцию аналитически.

Метод обратной функции

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

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

Так как в интервале всех возможных значений функция распределения монотонно возрастает от 0 до 1, то в этом интервале существует, причем только одно, такое значение аргумента при котором функция распределения примет значение . Другими словами, уравнение имеет единственное решение:

где — функция, обратная функции . 

Например, для стандартного нормального распределения функции и представлены на рис. 4б и 4в. Не всегда (более того, лишь изредка) уравнение удается решить в явном виде относительно . Начнем с примера, когда случайная величина распределена по показательному закону, заданному функцией распределения:

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

Решим это уравнение относительно :

или

Отсюда

Случайное число заключено в интервале ; следовательно, число также случайное и принадлежит интервалу . Другими словами, величины и распределены одинаково. Поэтому для отыскания хi можно воспользоваться более простой формулой:

Это и есть наш генератор случайной величины, распределенной по экспоненте:

При распределении Лапласа (двойного экспоненциального) плотность вероятности случайной величины имеет вид :

где θ — параметр сдвига (медиана распределения), – параметр масштаба (ширина распределения), .

Распределение Лапласа с параметром называют стандартным. Его плотность вероятности:

Интегрируя , получаем функцию распределения:

Решая уравнения относительно , получаем обратную функцию:

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

Стандартное распределение Коши:

Стандартное распределение Вейбулла :

, 

Стандартное логистическое распределение:

Стандартное распределение Чампернауна:

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

Метод исключения

Также встречаются названия метод выборки с отклонением (Rejection Sampling), метод отбраковки, метод принятия-отклонения (Acceptance-Rejection Method), метод фильтрации и алгоритм Неймана (по имени автора).

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

для всей области определения . 

Для простоты в качестве можно взять функцию, распределенную равномерно на области определения и проходящую через максимум функции . Это легко сделать, если область определения ограничена. Если же область определения не ограничена, нужно задать границы принудительно.

Алгоритм метода включает три шага.

1. Сгенерируйте кандидата из плотности .

2. Сгенерируйте случайное число из равномерного распределения .

3. Проверьте условие . Если оно выполняется, примите , если нет, отклоните и перейдите к шагу 1.

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

Теоретическое обоснование метода исключения и много чего еще интересного по теме есть в Хабра-заметке Александра Самарина «Генераторы непрерывно распределенных случайных величин».

Разберем алгоритм в Excel на примере стандартного нормального распределения.

Пояснения модели Excel

В ячейках В2:В3 задаю границы области определения и , значения и . В5:В6 — проверяю среднее значение и стандартное отклонение массива значений , полученных методом исключения. В7:В8 — проверяю среднее значение и стандартное отклонение массива значений , полученных с помощью функции НОРМ.СТ.ОБР(). D2 — разыгрываю аргумент функций и , случайное число , распределенное равномерно на интервале от до . E2 — вычисляю значение для , где — стандартное нормальное распределение F2 — вычисляю значение , равное максимуму функции . G2 — разыгрываю случайное число , распределенное равномерно на интервале . H2 – если , возвращаю , иначе «нет». I2— разыгрываю случайное число , распределенное по стандартному нормальному распределению, используя функцию Excel «из коробки» НОРМ.СТ.ОБР(). K2 — задаю карманы (диапазоны) для построения диаграммы. L2 — в каждом диапазоне подсчитываю число разыгранных значений по методу исключения. M2 — в каждом диапазоне подсчитываю число разыгранных значений с помощью функции НОРМ.СТ.ОБР().

Напоследок рассмотрим интересную практическую задачу разыгрывания случайной величины.

Скорость сходимости центральной предельной теоремы

Год назад прочитал книгу Нассима Талеба «Статистические последствия жирных хвостов». Книга о математике, лежащей в основе историй Талеба, рассказанных в его предыдущих эссе. В частности автор показывает, что распределения с толстыми хвостами довольно медленно сходятся к нормальному в соответствии с центральной предельной теоремой (ЦПТ):

При конечной дисперсии случайных величин , распределение суммы этих случайных величин , нормированное на , в пределе стремится к нормальному распределению.

Быстрая сходимость равномерного распределения

Если случайная величина равномерно распределена на отрезке , плотность вероятности будет постоянной и равной . Поскольку я использовал диапазонов, вероятность попасть в один из них равна или 1% ( на рис. 13). Теперь добавим к другую случайную величину , независимую и с таким же распределением. У суммы распределение будет другим! График плотности вероятности для суммы стал треугольным. Добавим еще одну переменную, и плотность вероятности для распределения суммы станет колоколом. Нам достаточно 3–4 слагаемых, чтобы распределение приняло вид нормального:

В Excel я разыграл случайную равномерно распределенную на отрезке [0,1] величину формулой:

. 

Для генерации  я использовал среднее двух массивов =, и т. д.

Медленная сходимость распределения Парето

Плотность вероятности распределения Парето:

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

Стандартное распределение Парето определено для на интервале :

Для получаем:

А обратное распределение имеет вид:

Рассмотрим случайную величину , где каждое по отдельности — независимая случайная величина, распределенная по (3.4). Посмотрим, с какой скоростью сходится к нормальному распределению при росте . Для разыгрывания случайной величины используем сумму обратных распределений , деленную на .

В Excel я применил формулу (авторство принадлежит AlienSx с форума Планета Excel):

где — размер выборки,  — число слагаемых, каждое из которых соответствует , подробности во вложенном файле «15. ЦПТ. Парето, альфа 2.xlsx».

Центральная предельная теорема работает, но не так быстро, как ожидалось. Как говорит Насим Талеб: «Распределение Парето так и не приблизилось к гауссиане, хотя при это произойдет — если у вас хватит терпения и вы будете жить долго-долго».

Вывод

Итак, мы рассмотрели, как с помощью Excel моделировать (разыгрывать) случайные величины с различными законами распределения. Выяснили, что для этого доступны встроенные функции, формулы на основе обратных функций или метод исключения. Оказывается, что возможностей Excel в целом ряде прикладных задач будет вполне достаточно без необходимости прибегать к программированию. 

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