Алгоритмы генерации случайных величин с различными распределениями

Введение

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

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

 

1. Общие подходы к генерированию случайных величин

 

Основной составляющей, необходимой  для каждого метода генерирования случайных величин из любого распределения или случайного процесса, является источник независимых и одинаково распределенных случайных величин с распределением U(0,1) (Равномерно распределенная случайная величина на отрезке [0,1]). Существуют несколько альтернативных алгоритмов, которые можно использовать для генерирования случайных величин из заданного распределения: метод обратного преобразования, метод композиций, метод свёртки и метод принятия-отклонения.

 

1.1. Метод обратного преобразования

Предположим нам необходимо генерировать случайную величину X, являющуюся непрерывной и имеющую функцию распределения F-непрерывную и строго возрастающую, когда 0<F(x)<1. (если < и 0< ≤ <1, то < ). Примем что, - это обратная функция F. Тогда алгоритм для генерирования случайной величины X c функцией распределения F будет следующим:

1. Генерируем - (0,1)
2. Возвращаем X=

 всегда будет определено, поскольку 0≤U≤1. А областью функции F является [0,1]. На рисунке 1 этот алгоритм изображен графически, на нем случайная величина, соответствующая этой функции распределения, может принимать либо положительные, либо отрицательные значения; это зависит от конкретного значения U. На рисунке 1 случайное число в результате дает положительную конкретную случайную величину  , тогда как случайное число в результате дает случайную величину

рис. 1

Если величина X имеет нужную функцию распределения F, для любого вещественного числа х справедлива формула

Пример 1. Пусть Х имеет экспоненциальное распределение со средним β. Функция распределения

 

Задаем u=F(x), чтобы найти , и решаем уравнение для х, желая получить (В этом  случае можно использовать U вместо 1-U,  поскольку 1-U и U имеют одно и тоже распределение U(0,1) ) . Следовательно, получаем величину

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

1. Генерируем U – U(0,1)
2. Определяем положительное наименьшее число I, для которого   и возвращаем .

На рис.2 показан метод, с помощью  которого в данном случае мы генерируем  .

Рис.2

Чтобы убедиться в достоверности  метода обратного преобразования для  дискретных величин, необходимо доказать, что  для всех значений i. Для i=1 тогда и только тогда мы получем ,когда , так как мы расположили величины в порядке возрастания. Так как U – U(0,1), то , что и требовалось доказать. Для алгоритм устанавливает тогда  и только тогда, когда , поскольку значение i, выбираемое алгоритмом, является положительным наименьшим числом, для которого . Более того поскольку U – U(0,1) и ,то

 

 

1.2. Метод композиции.

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

  В частности, допускаем, что  для всех х функцию F(x) можно записать как

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

1. Генерируем положительное целое число J, для которого
2. Возвращаем Х с функцией распределения

Шаг 1 можно рассматривать как  выбор функции распределения  с вероятностью , его можно выполнить, например, с помощью метода обратного преобразования для дискретных величин. При условии что J=j, генерирование на шаге 2 должно выполняться независимо от J. Определившись со значением J, генерируемым на шаге 1, можем убедиться, что величина X, возвращаемая алгоритмом, будет иметь функцию распределения F:

 

 

 

Можно также дать геометрическую интерпретацию  метода композиции. Например, для непрерывной  случайной величины Х с плотностью распределения f следует разделить площадь под функцией f на участки площадью соответствующие декомпозиции f ее представлением выпуклой комбинацией. После этого шаг 1 можно рассматривать как выбор участка, а шаг 2 – как генерирование величины из распределения, соответствующего выбранному участку.

  Пример 2.  Распределение Лапласа имеет плотность для всех вещественных значений х. График плотности экспоненциального распределения изображен на рис.3 По графику можно увидеть, что f(x) представляет собой две плотночти экспоненциального распределения(за исключением нормализующего множителя 0,5), расположившиеся вплотную, что предполагает использование метода композиции. На самом деле плотность распределения можно выразить так

, где  обозначает индикаторную функцию множества А, определенную формулой

Таким образом, f(x) является выпуклой комбинацией и , обе функции представляют собой функции плотности, а . Следовательно, мы можем генерировать величину /X с плотностью распределения f с помощью метода композиции. Вначале генерируем и как независимые и одинаково распределенные величины с распределением U(0,1). Если . Возвращаем

Рис.3

 

 

 

 

1.3. Метод свертки.

Для некоторых распределений искомую  случайную величину X можно выразить как сумму других независимых и одинаково распределенных случайных величин, которые проще генерировать, чем саму величину X. Допускается, что существуют независимые и одинаково распределенные случайные величины для заданного значения m, их сумма имеет то же распределение, что и величина X. Поэтому записываем  . Этот метод называется сверткой.

Алгоритм генерирования случайной  искомой величины X легко вывести интуитивно (пусть F будет функцией распределения величины X, а G – функцией распределения величины ).

1. Генерируем  как независимые и одинаково распределенные случайные величины с функцией распределения G каждая.
2. Возвращаем

       Чтобы доказать обоснованность этого алгоритма, используем то свойство, что X и  имеют одну и ту же функцию распределения – функцию F. Следовательно,  

Метод свертки, если им можно воспользоваться, очень прост для выполнения при  условии, что величины легко генерируются. Однако в зависимости от определенных параметров распределения величины X он может оказаться не самым эффективным методом генерирования. Например, генерирование случайной величины с распределением Эрланга m- го порядка с помощью метода свертки, когда m имеет большое значение, может потребовать много времени.

 

1.4. Метод принятия-отклонения

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

Для метода принятия-отклонения требуется  определить функцию  t, которая оценивает сверху плотность распределения f,  то есть для любых x. Теперь t не будет выступать плотностью распределения, поскольку

, а плотностью будет функция  r(x)=t(x)/c. Нам необходимо иметь возможность генерировать случайную величину Y c плотностью распределения r. Алгоритм:

1. Генерируем Y с плотностью распределения r.
2. Генерируем U-(0,1) независимо от Y
3. Если , возвращаем X=Y. В противном случае возвращаемся к шагу 1 и начинаем все сначала.

     Алгоритм продолжает  возвращаться к шагу 1 от тех  пор, пока на шагах 1 и 2 не  будет сгенерирована пара (Y,U), для которой , и тогда мы «примем» значение Y в качестве X.

 

2. Генерирование непрерывных случайных величин

 

2.1. Равномерное распределение (U(a,b))

Плотность
Функция распределения
Параметры	a и b – вещественные числа, для которых a<b; a – параметр положения, b – масштабный параметр
Область	[a,b]

 

рис.2.1 Плотность U(a,b)

:

  Функцию распределения случайной величины U(a,b) легко обратить, решив уравнение u=F(x), чтобы получить x для , 

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

1. Генерируем U-(0,1)
2. Возвращаем 

Если нужно генерировать много  значений X, константа b-a должна быть вычислена заранее.

 

2.2. Экспоненциальное распределение (expo(β))

Плотность
Функция распределения
Параметры	Масштабный  параметр β>0
Область	[0,∞)

Рис 2.2 Плотность expo(1)

Алгоритм:

1. Генерируем U-U(0,1)
2. Возвращаем

 

 

 

2.3. Гамма-распределение (Gamma(α))

Плотность
Функция распределения	Если α не является целым числом, конечной формы не существует. Если α-положительное целое число, тогда
Параметры	Параметры формы α>0, масштабный параметр β>0
Область	[0,∞)

Рис.2.3  Плотность gamma(α,1)

 

Случайные величины с гамма-распределением сложнее генерировать, нежели три  типа случайных величин, поскольку  функция распределения не имеет  для них конечной формы, для которой мы могли бы попробовать обратную функцию. Во-первых, при заданном X – gamma(α,β) для любого β>0 можно получить случайную величину X’,задав X’=βX, так что достаточно ограничить рассмотрение генерированием величин из распределения  gamma(α,1). Во-вторых, распределение gamma(1,1) является экспоненциальным со средним 1, то есть нам понадобится рассматривать только 0<α<1 и α>1. Поскольку алгоритмы для генерирования случайных величин с гамма-распределением, как правило, действительны только в одной из этих областей значений α, будем рассматривать отдельно.

Рассмотрим случай, когда 0<α<1. Алгоритм представляет собой метод принятия-отклонения с оцениваемой сверху функцией

Таким образом, , в результате получаем плотность как

Случайную величину Y с плотностью r(x) можно генерировать с помощью метода обратного преобразования; при этом функция распределения, соответствующая r, будет иметь вид

Ее можно обратить, чтобы получить

 

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

1. Генерируем и пусть . Если P>1, переходим к шагу 3. В противном случае переходим к шагу 2.
2. Допускаем , генерируем . Если , возвращаем X=Y. В противном случае возвращаемся к шагу 1.
3. Допускаем , генерируем . Если . Возвращаем X=Y. В противном случае переходим к шагу 1.

    Далее рассмотрим случай, когда  , который представляет собой модифицированный метод принятия-отклонения. Время выполнения этого алгоритма связано с , он работает быстрее по мере увеличения значения .(Модификация заключается в добавлении быстрой предварительной проверки на принятие). Для того чтобы получить оцениваемую сверху функцию t(x), сначала определим . Затем определяем t(x)=cr(x),где

Функция распределения, соответствующая  плотности r(x), составляет

Ее легко инвертировать, чтобы  получить

Алгоритм:

1. Генерируем  и как независимые и одинаково распределенные случайные величины.
2. Допускаем .
3. Если , возвращаем X=Y. В противном случае, переходим к шагу
4. Если , возвращаем X=Y. В противном случае переходим к шагу 1.

Шаг 3 – это добавленная предварительная  оценка, в случае прохождения которой  можно избежать вычисления логарифма  на шаге 4 обычного метода принятия-отклонения.

 

 

 

 

2.4. Распределение Вейбулла (Weibull(α,β))

 

Плотность
Функция распределения
Параметры	Параметр формы  , масштабный параметр
Область	[0,∞)

Рис. 2.4 Плотности Weibull(α,1):

 

Функцию плотности распределения  Вейбулла легко обратить, чтобы найти

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

1. Генерируем U-U(0,1)
2. Возвращаем

 

Снова используем тот факт, что 1-U и U имеют одно и то же распределение U(0,1). Этот алгоритм можно также обосновать, отметив, что если Y имеет экспоненциальное распределение со средним , то

 

 

 

 

 

 

 

2.5. Нормальное распределение (N(

))

Плотность
Функция распределения
Параметры	Параметр положения  , масштабный параметр
Область	(-∞,∞)

 

Рис.2.5 Плотность N(0,1)

 

Прежде всего, стоит обратить внимание, что при можно получить

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

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

1. Пусть задана функция распределения 

   и 

тогда , где - равномерно распределена на отрезке [0,1], будет распределена как N(0,1).

 

2.  Пусть две случайные величины независимы и равномерно распределены на [0,1] , тогда и   будут независимыми гауссовскими случайными величинами с параметрами , =0 и , =1.

 

3. Данный метод основан на центральной предельной теореме. Пусть задана последовательность независимых, равномерно распределенных случайных величин . Тогда случайная величина при будет сходиться к , где . Так как , а , то , которая в свою очередь будет распределена как N(0,1).

 

 

 

 

 

 

Плотность
Функция распределения
Параметры	Параметр положения  , масштабный параметр
Область	[0,∞)

2.6. Логнормальное распределение (LN( ))

Рис.2.6 Плотность

 

Специальное свойство логнормального распределения – если , то - используется чтобы получить следующий алгоритм

1. Генерируем
2. Возвращаем 

Параметры не являются средним и дисперсией распределения LN( ). В действительности,  если мы и мы определяем , то получается, что . Следовательно, если мы намерены генерировать логнормально распределенную случайную величину X с заданным средним и дисперсий , необходимо сначала определить  в исчислении . Перед тем как генерировать промежуточную нормально распределенную случайную величину Y. Получаем следующие формулы:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Плотность
Функция распределения
Параметры	Параметры формы
Область	[0,1]

2.7. Бета-распределение beta( )

 

Рис.2.7 Плотность beta (

)  а-г – при разных значениях

 

 

Мы можем получить beta( ) в интервале [a,b] для a<b, установив X’=a+(b-a)X, где X’ – beta( ) в интервале [0,1]. Так что достаточно рассмотреть последний случай.

Общий метод генерирования случайных  величин с распределением beta( ) для любых значений основывается на том факте, что если и являются независимыми величинами, то - beta( ). Таким образом, получаем следующий алгоритм:

1. Независимо генерируем и
2. Возвращаем

Этот метод достаточно удобен, поскольку он, по сути, выполняется для всех значений α>0 при условии, что есть генераторы случайных чисел с распределением Г( ). Эффективность его выполнения зависит от скорости работы выбранных генераторов случайных величин с гамма-распределением.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2.8. Распределение Пирсона типа V (PT5(α,β))

Плотность
Функция распределения
Параметры	Параметры формы
Область	[0,∞)

 

Рис.2.8 Плотность PT5(α,1)

, когда 1/X – gamma( ). Таким образом, мы получаем следующий алгоритм, основанный на специальном свойстве этого распределения.

1. Генерируем Y - gamma( ).
2. Возвращаем X=1/Y

 

 

 

 

 

2.9. Распределение Пирсона VI  (PT6(

))

Плотность
Функция распределения
Параметры	Параметры формы  , масштабный параметр β>0
Область	[0,∞)

 

 

 

 

 

 

Рис. 2.9а Плотность PT6(

)

a -

; б -

Рис. 2.9б Плотность PT6(

)

в -

; г -

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

1. Независимо генерируем  ,
2. Возвращаем 

 

 

 

 

 

2.10. Лог-логистическое распределение (LL(α,1))

Плотность
Функция распределения
Параметры	Параметр формы α>0, масштабный параметр β>0
Область	[0,∞)

 

 

 

 

 

 

 

 

 

 

 

 

Рис. 2.10 Плотность LL(α,1)

 

Функцию лог-логистического распределения легко обратить, чтобы получить

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

1. U – U(0,1);

2.

 

2.11. Связанное распределение Джонсона (JSB(

, ,a,b))

Плотность
Функция распределения
Параметры	Параметр  , масштабный параметр b-a (b>a) , параметр формы и
Область	[a,b]

 

 

Рис. 2.11 Плотность JSB(

, ,1) a,б – при разных значениях ,

 

Случайная величина X – JSB( , ,a,b),тогда и только тогда, когда               . Решим это уравнение для X в исчислении Z, чтобы получить следующий алгоритм, основанный на специальном свойстве.

1. Генерируем Z – N(0,1)
2. Пусть
3. Возвращаем 

 

 

 

 

 

 

 

 

 

 

 

2.12 Несвязанное распределение Джонсона (JSU(

)

Плотность
Функция распределения
Параметры	Параметр положения  , масштабный параметр , параметр формы и
Область	( )

 

Рис.2.12 Плотность JSU(

)

a – при

, б – при

 

Случайная величина , тогда и только тогда, когда

Решим это уравнение для X в исчислении Z, чтобы получить следующий алгоритм, основанный на специальном свойстве распределения.

1. Генерируем  Z – N(0,1)
2. Допускаем, что
3. Возвращаем

 

 

 

2.13. Треугольное распределение (triang(a,b,c))

Плотность
Функция распределения
Параметры	Вещественные числа a,b,c для которых a<c<b; a – параметр положения, b-a – масштабный параметр, c  - параметр формы
Область	[a,b]

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Рис.2.13 Плотноcть triang(a,b,c)

 

Если X – triag[0,1,(c-a)/(b-a)], X’=a+(b-a)X . Поэтому можно ограничиться рассмотрением случайных величин с распределением triag(0,1,c), где 0<c<1. Функцию треугольного распрделения легко обратить, чтобы получить для

 

 

 

Сформулируем следующий алгоритм обратного преобразования для генерирования    X – triag(0,1,c):

1. Генерируем U – U(0,1)
2. Если , возвращаем . В противном случае возвращаем

 

14. Прочие распределения

 

1. Распределение Эрланга m-го порядка.

   Случайная величина X c распределением Эрланга m-го порядка и средним значением β может быть определена как сумма m независимых и одинаково экспоненциально распределенных случайных величин с общим средним β/m. Значит, чтобы генерировать случайную величину X, можно генерировать  как независимые и одинаково распределенные случайные величины со средним β/m, а затем возвратить 

 

Рис.2.14 Плотность Erlang(1)

2. Распределение .

Если Y – N(0,1), то будет иметь распределение «хи-квадрат» с одной степенью свободы.

      Если  являются независимыми и одинаково распределенными величинами с распределением , то

Рис.2.15 Плотность  Hi?(4)

3. Распределение Стьюдента

Если Y – N(0,1), , а Y и Z являются независимыми, то считается, что имеет распределение Стьюдента с k степенями свободы.

Рис. 2.16 Плотность

  k=1.

 

 

3. Генерирование дискретных случайных величин

 

1. Распределение Бернулли (Bernulli(p))

Вероятностная мера
Распределение
Параметры
Область

Рис.3.1 Вероятностная мера Bernulli(p)

 

1. Генерируем U – U(0,1)
2. Если , возвращаем X=1. В противном случае возвращаем X=0.

 

2. Дискретное равномерное распределение (DU(i,j))

 

 

Вероятностная мера
Распределение
Параметры	i,j – целые числа, для которых , i – параметр положения, j-i – масштабный параметр
Область	{i,i+1…j}

 

Рис. 3.2 Плотность DU(i,j)

1. Генерируем U – U(0,1)
2. Возвращаем

 

 

3. Биномиальное распределение (Bin(t,p))