Моделювання неперервних випадкових величин

Існує кілька методів моделювання значень неперервних випадкових величин з до­вільним законом розподілу на основі випадкових чисел, рівномірно розподілених у інтервалі [0, 1]: метод оберненої функції, метод відсіювання, наближені методи тощо.

4.7.1 Метод оберненої функції

Розглянемо метод моделювання випадкової величини, яка має функцію щільності ймовірностей f(x) і монотонно зростаючу функцію розподілу F(х) (рис. 4.9). Суть методу така. За допомогою генератора випадкових чисел генеруємо зна­чення випадкової величини r _i ·, якому відповідає точка на осі ординат. Значення випадкової величини х_i з функцією розподілу F (х) можемо одержати з рівняння F(x_i) = r _i.

Дійсно, якщо на осі ординат відкласти значення r _i випадкової величини, розподіленої рівномірно в інтервалі [0, 1], і на осі абсцис знайти значення х_i випадкової величини (рис. 4.6), при якомуF(x _i) = r _i, то випадкова величина X = F^-1 (r) буде мати функцію розподілу F(x).

Таким чином, послідовність випадкових чисел r₁, r ₂, r₃,... перетворюється на послідовність х_и х₂, х₃,..., яка має задану функцію щільності розподілу f(x).

Рисунок 4.6 – Використання методу оберненої функції для генерування неперервної випадкової величини

Звідси випливає загальний алгоритм моделювання випадкових неперервних величин, що мають задану функцію розподілу ймовірностей:

- генерується випадкове число r _i є [0, 1];

- обчислюється випадкове число х_i яке є розв'язком рівняння

Приклади застосування методу наведені нижче.

4.7.2 Рівномірний розподіл

У загальному випадку випадкова величина X є рівномірно розподіленою на від­різку [а, b], якщо її щільність розподілу ймовірностей має вигляд

Функцію розподілу ймовірностей можна знайти як

тобто

Математичне сподівання та дисперсія випадкової величини X визначаються як

Для моделювання випадкової рівномірно розподіленої на відрізку [а, b] величини можна скористатись методом оберненої функції. Обчислимо функцію роз­поділу випадкової величини та прирівняємо її до значення r_i.

Звідси знаходимо значення випадкової величини з функцією розподілу f (x):

Цю формулу також можна отримати, якщо виконати лінійне перетворення інтервалу [0,1] у відрізок [а, b]. Для цього потрібно змінити масштаб функції рівномірного розподілу, помноживши її на (b - а), а потім змістити її на величину а.

Функція рівномірного розподілу широко застосовується для моделювання випадкових величин, для яких функція розподілу невідома, а відоме лише її середнє значення. У такому випадку припускають, що відомими є середнє значення випадкової величини та деяке розсіювання (+ -Δ) її значень відносно середнього. Це дає змогу стверджувати, що дана випадкова величина має рівномірний розподіл.

Прикладами реальних задач, в яких виникає необхідність моделювання рівномірно розподілених випадкових величин, можуть бути аналіз помилок округлення під час проведення числових розрахунків (точність задається як кількість десяткових знаків), час переміщення головок у магнітних накопичувачів (мінімальний та максимальний час), відхилення від розкладу руху транспортних засобів (наприклад, метро).

4.7.3 Експоненціальний розподіл

Експоненціальний закон розподілу набув широкого використання в теорії надійності складних систем. Функція щільності експоненціального розподілу випадкової величини має вигляд

Для її моделювання скористаємося методом оберненої функції. Маємо

(4.6)

З виразу (4.6) знаходимо значення х_i.

Можна показати, що випадкові величини (1 - r_i) мають такий самий розподіл, що і величини r_i. Тоді, замінивши 1 - r_i на r_i отримаємо

.

Випадкові величини з експоненціальним розподілом широко застосовуються в задачах моделювання та аналізу СМО, наприклад під час моделювання процесів виходу з ладу та ремонту обладнання, які виникають у складних системах, у разі визначення інтервалів часу між послідовними викликами абонентів у телефонній мережі або замовлень від незалежних клієнтів у будь-якій мережі обслуговування (швидка допомога, служби ремонту, виклик таксі і т. ін.)

4.7.4 Пуассонівський потік

Розглянемо моделювання пуассонівського потоку з інтенсивністю λ, основна властивість якого полягає в тому, що ймовірність надходження k вимог протягом інтервалу довжиною t становить

Для пуассонівського потоку інтервали часу між надходженням двох сусідніх вимог мають експоненціальний закон розподілу (див. розділ 2.1). Тому для його моделювання достатньо отримати ряд чисел з таким розподілом. Це можна реалі­зувати за допомогою методу оберненої функції, якщо ряд випадкових чисел r _i рівномірно розподілених у інтервалі [0, 1], перетворити згідно з функцією, оберненою до експоненціальної функції розподілу

де t_j — j-й проміжок часу між надходженнями двох сусідніх вимог; Τ = 1/λ — середнє значення проміжку часу між надходженнями двох сусідніх вимог; r_j — j-e число в послідовності випадкових чисел з рівномірним розподілом у інтервалі [0, 1].

4.5.5 Нормальний розподіл

Випадкова величина X має нормальний розподіл (розподіл Гаусса), якщо її щільність

де т — математичне сподівання, а σ — середньоквадратичне відхилення розподілу ймовірностей описується виразом

Функція розподілу нормально розподіленої величини X має вигляд

Для моделювання випадкової величини з нормальним законом розподілу безпосередньо скористатися методом оберненої функції не можна, оскільки немож­ливо аналітично виконати перетворення виду X = F^-1(r). Тому для моделювання слід скористатися методом згорток.

Метод згорток базується на центральній граничній теоремі — одному із найбільш видатних результатів теорії ймовірностей: за широких припущень відносно розподілів суми великої кількості взаємно незалежних малих випадкових величин має місце розподіл, який близький до нормального. Метод згорток передбачає зображення випадкової величини як суми незалежних однаково розподілених випадкових величин зі скінченними математичним сподіванням і дисперсією.

Найпростіший метод отримання значення випадкової величини, що має заданий нормальний розподіл, передбачає виконання таких кроків. Спочатку формують послідовність r_i (і = 1 ,n) незалежних, рівномірно розподілених у інтервалі [0, 1] величин і обчислюють суму 12 випадкових чисел, потім віднімають число 6. Величина п= 12 є хорошим наближенням до нормально розподіленої випадкової величини з нульовим математичним сподіванням m_z = 0 і одиничним середньоквадратичним відхиленням σ_ζ = 1. Нормальний розподіл з параметрами т_г - 0 та σ_ζ = 1 називається стандартним.

Недоліком розглянутих вище методів моделювання є те, що значення функції нормального розподілу, які лежать за межами т_х ± σ_х, суттєво відрізняються від точних значень. Щоб зменшити загальну похибку моделювання, треба використовувати більш точні методи отримання значень функції нормального розподілу. Ці методи базуються на такій властивості. Якщо X_t і Х₂ є незалежними нормально розподіленими випадковими величинами з нульовим математичним сподіванням і одиничним середньоквадратичним відхиленням, то величина кута між віссю абсцис і вершиною випадкового вектора з координатами {х\, х₂) має рівномірний розподіл і не залежить від довжини вектора () (рис. 4.7).

Х2

Рисунок 4.7 - Зображення вектора

для моделювання нормального розподілу

Квадрат довжини вектора в цьому випадку має розподіл χ² з двома ступеня­ми свободи і моделюється як окремий випадок показового розподілу з параметром λ = 1/2.

Існує два методи моделювання нормального розподілу, які використовують цю властивість:

Метод Бокса-Мюллера (Box-Muller). Генеруємо пару нормально розподілених чисел з т_х = 0 і σ_χ = 1 за допомогою двох випадкових чисел r_t і г₂:

Х₁ = - 2 In r₁ cos (2π r₁); χ₂ = - 2 In r₂ cos (2π r₂).

Таким чином, отримуємо два числа Х1 і х2 нормальним розподілом. Метод Марсальї-Брея (Marsaglia-Bray). Існує більш швидка модифікація цього методу. Генерують два випадкових числа г\ і г₂, вважаючи, що ν₁ = -1 + 2r_{1 ;}v₂ = -1 + 2r₂, обчислюють суму S = v₁ + v₂. Якщо S > 1, то повторюють процедуру, якщо S < 1, то одержують два нормально розподілених числа:

,

Щоб одержати за цим методом 100 пар нормально розподілених чисел, потрібно генерувати 127 пар випадкових чисел. Це простий та швидкий метод, у разі його застосування більша частина часу роботи алгоритму витрачається на обчислення логарифму.

Розподіл і потоки Ерланга

Випадкові величини з експоненціальним розподілом не завжди адекватно описують деякі реальні процеси та події, наприклад час обслуговування і моменти над­ходження вимог до СМО. Для більш точного моделювання таких процесів доцільніше використовувати гамма-розподілені випадкові величини або ті, що мають розподіл Ерланга. Розподіл Ерланга є результатом підсумовування взаємно незалежних і однаково розподілених експоненціальних випадкових величин і є окремим випадком гамма-розподілу.

Функція щільності розподілу Ерланга k-το порядку з інтенсивністю λ має та­кий вигляд

Математичне сподівання і дисперсія розподілу Ерланга визначаються як М[х] = 1 /k λ, D [x]= 1 /k λ².

Для моделювання розподілу Ерланга використовують метод згорток випадкових величин з експоненціальними функціями розподілу. Для цього треба лише обчислити суму k експоненціально розподілених випадко­вих величин. Зі збільшенням k розподіл Ерланга наближається до нормального.

4.8 Статистична обробка результатів моделювання

Основою для обчислення статистичної оцінки параметра системи є реалізація випадкової величини, яка формується під час прогонів імовірнісної імітаційної моделі. Статистична оцінка також є функцією від випадкових величин, які діста­ють у результаті прогонів моделі, тому і згадана оцінка є випадковою величиною, закон розподілу якої залежить від закону розподілу досліджуваної випадкової ве­личини та оцінюваного параметра. Чим більше реалізацій випадкової величини, тим точнішу статистичну оцінку параметра системи ми отримуємо. У таких умовах обробка результатів моделювання повинна провадитись лише з використанням методів і алгоритмів, які є оптимальними з погляду затрат часу та використання ресурсів комп'ютера. Під час вибору таких засобів необхідно враховувати, що всі статистичні оцінки мають бути ще й незміщеними, ефективними і спроможними.

Розглянемо деякі методи обчислення основних статистичних оцінок.

4.8.1 Оцінювання ймовірності

Оцінкою ймовірності ρ настання деякої події А є її частість:

де т — кількість випробувань, під час за яких випадкова подія спостерігалась; N — загальна кількість випробувань. Для її використання зазвичай на програмному рівні організовують два лічильники, один з яких призначено для підрахунку загальної кількості випробувань Ν, а другий — кількості успішних випробувань т.

4.8.2 Оцінювання розподілу випадкової величини

Для оцінювання функції розподілу випадкової величини, як звичайно, будується гістограма. Під час її побудови область можливих значень випадкової величини розбивають на η діапазонів і підраховують кількість попадання значень випадкової величини в конкретний інтервал - т_k (k = 1,..., n ). Оцінка ймовірності по­падання випадкової величини в k-й. інтервал має такий вигляд:

Цю величину називають відносною частістю. У процесі моделювання під час підрахунку значень т_k кожному інтервалу ставлять у відповідність окремий лічильник.

4.8.3 Оцінювання математичного сподівання

Для оцінювання математичного сподівання випадкової величини використовується формула:

де х_k — значення випадкової величини, що належить k-му інтервалу; т_k — кількість попадань значень випадкової величини в інтервал; N — загальна кількість випробувань. У більш простому випадку для оцінювання математичного сподівання випадкової величини можна використати звичайне середнє арифметичне:

Щоб запобігти непотрібному завантаженню пам'яті, суму доцільніше підраховувати шляхом поступового накопичення.

4.8.4 Оцінювання дисперсії

Для оцінювання дисперсії випадкової величини можна використати формулу:

де S² — оцінка дисперсії випадкової величини х.

4.8.5 Оцінювання кореляційного моменту

Для обчислення оцінки кореляційного моменту можна використовувати формулу

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

4.9 Визначення кількості реалізацій під час моделювання випадкових величин

Точність оцінок параметрів системи, які отримують під час обробки результатів моделювання, у першу чергу залежить від кількості випробувань N. Слід враховувати, що обсяг вибірки N завжди обмежений, тому вищезгадані оцінки матимуть різні похибки і дисперсії.

Якщо треба оцінити значення деякого параметра а за результатами моделювання x_i, то за його оцінку слід брати величину х_і, яка є функцією від усіх значень х_i. Статистична оцінка xтакож є випадковою величиною, тому вона буде відрізнятись від х_i., тобто

де ε - точність або похибка оцінки. Імовірність того, що ця нерівність виконується, позначимо через α:

У теорії ймовірностей ε — це довірчий інтервал для а, довжина якого фактично дорівнює 2ε, а α: — довірчий рівень, або надійність оцінки. Вираз (4.13) можна застосувати для визначення точності результатів статистичних випробувань.

4.9.1 Оцінювання ймовірності

Припустимо, що метою моделювання є оцінювання ймовірності настання деякої події А, яка визначає стан системи. У кожній з N реалізацій процесу настання події А є випадковою величиною ξ, що набуває значення x₁= 1 з імовірністю ρ і x ₂ = 0 з імовірністю 1- р. Тоді можна визначити математичне сподівання і дисперсію відповідно за формулами

Μ[ξ]=x₁ρ+x₂(1-ρ)=ρ, (4.14)

D [ξ] = (x₁ - Μ [ξ])²ρ + (x₂ - Μ [ξ])² (1- ρ)=ρ (1- ρ).

(4.15)

Як оцінку р використовують частість настання події А. Ця оцінка є незміщеною, спроможною та ефективною. За умови, що N задано, для отримання цієї оцінки достатньо накопичувати т:

де x_i — настання події А в реалізації і.

За формулами (4.14), (4.16) визначимо вибіркове математичне сподівання

Μ [m/N] = ρ і дисперсію D [m/Ν] =ρ(1-ρ)/(Ν -1).

Згідно з центральною граничною теоремою (у даному випадку її можна взяти у вигляді теореми Хінчина) випадкова величина m/N буде мати розподіл, близький до нормального (рис. 4.8). Тому для кожного рівня достовірності α з таблиць нормального розподілу можна

знайти таку величину t_α, при якій точність обчислюватиметься за формулою

Якщо α = 0,05, то t_a = 1,96, а якщо α = 0,003, то t_α= 3.

Підставимо у формулу (4.17) вираз дисперсії . Звідси

(4.18)

Рисунок 4.8 - Функція нормального розподілу для побудови довірчого інтервалу

З формули (4.18) видно, що при ρ = 1 або ρ = 0, кількість реалізацій, які необхідно провести для підтвердження того, що подія А настає (або ні), дорівнює одиниці. Але оскільки ймовірність ρ заздалегідь невідома, провадять випробування (Ν = 50... 100), оцінюють частість m /Ν і підставляють її значення у вираз (4.18) замість р, після чого визначають остаточну кількість реалізацій. Графік залежності числа реалізацій для α = 0,05 і різних значень р, якщо ε = 0,05, наведено на рис. 4.9.

Рисунок 4.9 - Залежність числа реалізацій від значень імовірності

4.9.2 Оцінювання середнього значення

Нехай випадкова величина має математичне сподівання а і дисперсію σ². У і -й реалізації вона набуває значення х_i. Як оцінку математичного сподівання а використаємо середнє арифметичне:

(4.19)

Згідно з центральною граничною теоремою при великих значеннях N середнє арифметичне (4.19) буде мати нормальний розподіл з математичним сподіванням а і дисперсією σ²/(Ν -1). Тоді

Звідси

N = t_α ² σ²/ε² + 1 (4.20)

Оскільки дисперсія σ² випадкової величини невідома, потрібно провести кілька десятків (50... 100) випробувань і знайти оцінку σ², а потім отримане значення підставити у формулу (4.20), щоб визначити необхідну кількість реалізацій N. У цьому випадку замість нормально розподіленої величини необхідно скориста­тись t-розподілом Стьюдента з N - 1 ступенями свободи для визначення t_α. Зауважимо, що за збільшення ступенів вільності t-розподіл наближається до нормального. З практичного погляду, якщо N більше 30, користуються нормальним розподілом.