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

в 7:01, , рубрики: Ланит, моделирование в Excel, непрерывные распределения, случайная величина, теория вероятностей, численные методы

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

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

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

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

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

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

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

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

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

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

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

Изображение выглядит как линия, диаграмма, Прямоугольник, График  Автоматически созданное описание

Рис. 1. Теоретические кривые для непрерывной случайной величины Х, распределенной равномерно в интервале [0, 1]: (а) плотность вероятности, (б) функция распределения

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

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

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

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

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

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

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

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

 =СЛМАССИВ(n;;0;1;ЛОЖЬ). 

Здесь n – размер выборки, Как в Excel сгенерировать случайную величину произвольного распределения - 10 – минимальное значение, 1 – максимальное значение, ЛОЖЬ – задает плотность вероятности (ИСТИНА задала бы функцию распределения). В результате получим вертикальный массив из n десятичных значений в диапазоне от Как в Excel сгенерировать случайную величину произвольного распределения - 15 до 1 (рис. 2а).

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

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

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

(1.4) P(cleq Xleq d)=d-c

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

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

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

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

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

В настоящей заметке параметры генеральной совокупности известны. Их мы используем в формулах 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Х, СТЬЮДЕНТ.РАСП.ПХ. Для наших целей достаточно одного представителя от семейства. Рассмотрим структуру функций на примере пары НОРМ.СТ.РАСП и НОРМ.СТ.ОБР.

Синтаксис:

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

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

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

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

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

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

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

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

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

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

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

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

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

Рис. 6. Генерирование случайных величин на основе функций Excel с суффиксом ОБР; эмпирические плотности вероятности, n=1M

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

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

Пусть требуется разыграть непрерывную случайную величину X, т. е. получить последовательность ее возможных значений x_i (i=1, 2, ..., n), зная функцию распределения F(х).

Если p_i — случайное число в диапазоне (0; 1), то возможное значение x_i разыгрываемой непрерывной случайной величины X с заданной функцией распределения F(х), соответствующее p_i, является корнем уравнения F(х_i)=p_i.

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

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

где F^{-1} — функция, обратная функции y=F(x)

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

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

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

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

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

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

Отсюда

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

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

(2.6)x_i=-frac1λlnp_i

Изображение выглядит как линия, диаграмма, График, Прямоугольник  Автоматически созданное описание

Рис. 7. Теоретические кривые экспоненциального распределения: функция распределения F(x); обратная функция F^{-1}(p)

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

Изображение выглядит как текст, снимок экрана, Параллельный, линия  Автоматически созданное описание

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

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

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

где θ — параметр сдвига (медиана распределения), -∞ < θ < +∞; ϕ – параметр масштаба (ширина распределения), ϕ > 0.

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

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

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

(2.9) F_X(x)=begin{cases} frac12e^x, xleq 0 \  1-frac12e^{-x}, xgeq 0 end{cases}

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

(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. Генерирование случайной величины, распределенной по Лапласу, n=1M

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

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

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

Стандартное распределение Вейбулла (x ≥ 0):

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

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

(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), метод фильтрации и алгоритм Неймана (по имени автора).

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

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

для всей области определения f(x)

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

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

Рис. 11. Разыгрываемая f(x) и вспомогательная распределенная равномерно g(x)

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

  1. Сгенерируйте кандидата x_i из плотности g(x).

  2. Сгенерируйте случайное число u_i из равномерного распределения U[0;1].

  3. Проверьте условие u_i ≤ f(x_i)/g(x_i). Если оно выполняется, примите x_i, если нет, отклоните x_i и перейдите к шагу 1.

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

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

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

Изображение выглядит как текст, снимок экрана, число, документ  Автоматически созданное описание

Рис. 12. Реализация метода исключения в Excel, n=1M
Пояснения модели Excel

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

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

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

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

При конечной дисперсии случайных величин X_i, распределение суммы этих случайных величин φ_m=X_1 + X_2 + … + X_m, нормированное на m, в пределе стремится к нормальному распределению.

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

Если случайная величина X_1 равномерно распределена на отрезке [0, 1], плотность вероятности φ_1(х) будет постоянной и равной 1. Поскольку я использовал 100 диапазонов, вероятность попасть в один из них равна 1/100 или 1% (φ_1 на рис. 13). Теперь добавим к X_1 другую случайную величину X_2, независимую и с таким же распределением. У суммы X_1 + X_2 распределение будет другим! График плотности вероятности для суммы φ_2(x) стал треугольным. Добавим еще одну переменную, и плотность вероятности φ_3 для распределения суммы X_1 + X_2 + X_3 станет колоколом. Нам достаточно 3–4 слагаемых, чтобы распределение φ_m приняло вид нормального:

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

Рис. 13. Сумма независимых равномерных распределений быстро сходится к нормальному, n=1M

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

 =СЛМАССИВ(n;;0;1;ЛОЖЬ)

Для генерации  S_2=X_1 + X_2 я использовал среднее двух массивов =(СЛМАССИВ(n;;0;1;ЛОЖЬ)+СЛМАССИВ(n;;0;1;ЛОЖЬ))/2, и т. д.

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

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

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

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

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

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

Для α=2 получаем:

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

А обратное распределение имеет вид:
Как в Excel сгенерировать случайную величину произвольного распределения - 198

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

Рассмотрим случайную величину φ_m=X_1 + X_2 + … + X_m, где каждое по отдельности X_i — независимая случайная величина, распределенная по (3.4). Посмотрим, с какой скоростью φ_m сходится к нормальному распределению при росте m. Для разыгрывания случайной величины φ_m используем сумму m обратных распределений (3.5), деленную на m.

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

=BYROW(СЛМАССИВ(n;m;0;1;ЛОЖЬ);LAMBDA(массив;СУММ(1/(КОРЕНЬ(1-массив)))/m))

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

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

Рис. 15. Сходимость суммы стандартных распределений Парето с α = 2, n = 50k

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

Вывод

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

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

Автор: SergBag

Источник

* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js