SergBag 11 фев в 10:01

Как в Excel сгенерировать случайную величину произвольного распределения

Средний

10 мин

3.6K

Блог компании ГК ЛАНИТАлгоритмы*Математика*Программирование*

Туториал

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

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

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

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

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

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

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

$(1.1)f(x)= \begin{cases} 0 & \quad \text {при } х < a \\ \frac{1}{b-a} & \quad \text {при } a\leq x\leq b \\ 0 & \quad \text {при } x>b \end{cases}$

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

$(1.2) f(x)=\begin {cases} 0 & \quad \text {при } x<0 \\1 & \quad \text {при } 0\leq x\leq 1 \\ 0 & \quad \text {при } x>1 \end{cases}$

Изображение выглядит как линия, диаграмма, Прямоугольник, График Автоматически созданное описание — Рис. 1. Теоретические кривые для непрерывной случайной величины , распределенной равномерно в интервале *[0, 1]*: (а) плотность вероятности, (б) функция распределения

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

$(1.3)F_X(x)=P(X\leq x)=\int\limits_{-\infty}^0f(t)dt=\begin{cases} 0 & \quad \text {при } x<0 \\ x & \quad \text {при } 0\leq x\leq 1 \\ 1 & \quad \text {при } x>1 \end{cases}$

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

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

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

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

Рис. 2. Сгенерированное стандартное равномерное распределение: а) фрагмент массива; б) диапазоны и соответствующие им частоты; в) эмпирическая плотность вероятности;

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

$(1.4) P(c\leq X\leq d)=d-c$

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

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

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

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

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

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

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

$(1.5) \mu =\frac {a+b}{2}$

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

$(1.6)\mu =\frac{0+1}{2}=1/2$

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

$(1.7) \sigma^2=\frac{(b-a)^2}{12}$

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

$(1.8)\sigma =\sqrt{\frac{(1-0)^2}{12}}=1/\sqrt12$

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

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

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

Изображение выглядит как текст, снимок экрана, Шрифт, число Автоматически созданное описание — Рис. 3. Функции распределений в Excel

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

Синтаксис:

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

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

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

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

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

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

Рис. 5. Сгенерированное нормальное распределение: а) фрагмент массива; б) диапазоны и соответствующие им частоты; в) эмпирические плотности вероятности;

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

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

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

Изображение выглядит как диаграмма, График, линия, текст Автоматически созданное описание — Рис. 6. Генерирование случайных величин на основе функций Excel с суффиксом ОБР; эмпирические плотности вероятности,

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

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

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

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

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

$(2.1) x_i=F^{-1}(p_i)$

где $F^{-1}$ — функция, обратная функции .

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

$(2.2) F(x)=1-e^{-λx} (x>0)$

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

$(2.3) 1-e^{-λx_i}=p_i$

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

$(2.4) e^{-λx_i}=1-p_i$ или

Отсюда

$(2.5) x_i=-\frac1λln⁡(1-p_i)$

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

$(2.6)x_i=-\frac1λlnp_i$

Изображение выглядит как линия, диаграмма, График, Прямоугольник Автоматически созданное описание — Рис. 7. Теоретические кривые экспоненциального распределения: функция распределения ; обратная функция $F^{-1}(p)$

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

Изображение выглядит как текст, снимок экрана, Параллельный, линия Автоматически созданное описание — Рис. 8. Разыгрывание случайной экспоненциально распределенной величины в Excel; иногда для более точной визуализации распределения нужно поиграть с осью ОХ; здесь она сдвинута относительно карманов на 0,1

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

$(2.7) f(x)=\frac{1}{2\phi} e^{-\frac{|x-\theta|}{\phi}}, -∞<x<+∞$

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

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

$(2.8) f(x)=\frac12e^{-|x|}$

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

$(2.9) F_X(x)=\begin{cases} \frac12e^x, x\leq 0 \\ 1-\frac12e^{-x}, x\geq 0 \end{cases}$

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

$(2.10) F^{-1} (p)=\begin {cases} ln (2p), 0<p<\frac {1}{2} \\ ln(2-2p), \frac{1}{2}\leq p < 1 \end{cases}$

Изображение выглядит как текст, диаграмма, снимок экрана, линия Автоматически созданное описание — Рис. 9. Генерирование случайной величины, распределенной по Лапласу,

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

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

$(2.11) f(x) = \frac {1}{\pi(1+x^2)}; x_i=tg(\pi(p_i-0,5))$

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

$(2.12) f(x)= cx^{c-1}e^{-x^c};$ $x_i=\left\{-ln⁡(1-p_i)\right\}^{1/c}$ ,

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

$(2.13) f(x)= \frac{e^{-x}}{(1+e^{-x})^2}; x_i=ln⁡\frac{p_i}{1-p_i}$

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

$(2.14) f(x)= \frac{1}{πchx}; x_i=ln tg \frac{πp_i}{2}$

Изображение выглядит как текст, диаграмма, График, линия Автоматически созданное описание — Рис. 10. Некоторые аналитически представимые обратные распределения

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

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

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

Пусть есть случайная величина , имеющая плотность вероятности и функцию распределения , для которой не получается найти аналитическое выражение для $F^{-1}(p)$ . Подберем вспомогательную функцию с плотностью вероятности , для которой относительно легко разыграть случайную величину по вероятности , такую что

$(3.1) g(x)\geq f(x)≥0$

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

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

Изображение выглядит как линия, График, диаграмма, снимок экрана Автоматически созданное описание — Рис. 11. Разыгрываемая и вспомогательная распределенная равномерно

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

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

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

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

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

Изображение выглядит как текст, снимок экрана, число, документ Автоматически созданное описание — Рис. 12. Реализация метода исключения в Excel,

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

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

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

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

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

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

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

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

Изображение выглядит как диаграмма, График, текст Автоматически созданное описание — Рис. 13. Сумма независимых равномерных распределений быстро сходится к нормальному,

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

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

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

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

$(3.2) f(x|\alpha,x_{min})=\frac{\alpha x_{min}^\alpha}{x^{\alpha+1}}$

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

Стандартное распределение Парето определено для $x_{min} = 1$ на интервале :

$(3.3) f(x|\alpha), 1)=\alpha x^{-(\alpha+1)}$

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

$(3.4) f(x|2, 1)=2x^{-3}$

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

$(3.5) x_i=F^{-1}(p_i)=\frac{1}{\sqrt{1-p_i}}$

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

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

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

Изображение выглядит как диаграмма, График, дизайн Автоматически созданное описание — Рис. 15. Сходимость суммы стандартных распределений Парето с α = 2, n = 50k

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

Вывод

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

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

Теги:

Хабы: