Суть методу Монте-Карло та моделювання випадкових величин

Введення

Метод Монте-Карло - це чисельний метод рішення математичних задач за допомогою моделювання випадкових величин.

Датою народження методу Монте-Карло прийнято вважати 1949 р., коли з'явилася стаття під назвою «Метод Монте-Карло» (Н. Метрополіс, С. Улам). Творцями цього методу вважають американських математиків Дж. Неймана і З. Улама. У нашій країні перші статті були опубліковані в 1955-56 рр.. (В. В. Чавчанідзе, Ю. А. Шрейдер, В. С. Владимиров)

Проте теоретична основа методу була відома давно. Крім того, деякі задачі статистики розраховувалися інколи з допомогою випадкових вибірок, тобто фактично методом Монте-Карло. Проте до появи ЕОМ цей метод не міг знайти скільки-небудь широкого застосування, так як моделювати випадкові величини вручну - дуже трудомістка робота. Таким чином, виникнення методу Монте-Карло як дуже універсального чисельного методу стала можлива лише завдяки появі ЕОМ.

Сама назва «Монте-Карло» походить від міста Монте-Карло в князівстві Монако, знаменитого своїм ігорним домом, а одним з найпростіших механічних приладів для отримання випадкових величин є рулетка.

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

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

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

1. Теоретична частина

1.1 Сутність методу Монте-Карло та моделювання випадкових величин

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

Рис. 1.1

Припустимо, що ця фігура розташована всередині одиничного квадрата.

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

Для того щоб вибирати точки випадково, необхідно перейти до поняття випадкова величина. Випадкова величина безперервна, якщо вона може приймати будь-яке значення з деякого інтервалу .

Неперервна випадкова величина визначається завданням інтервалу , Що містить можливі значення цієї величини, і функції , Яка називається щільністю ймовірностей випадкової величини (Щільністю розподілу ). Фізичний сенс наступний: нехай - Довільний інтервал, такий що , Тоді ймовірність того, що виявиться в інтервалі , Дорівнює інтегралу

(1.1)

Безліч значень може бути будь-яким інтервалом (можливий випадок ). Однак щільність повинна задовольняти двом умовам:

1. щільність позитивна:

; (1.2)

2. інтеграл від щільності по всьому інтервалу дорівнює 1:

(1.3)

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

(1.4)

Дисперсією випадкової величини називається число:

Нормальною випадковою величиною називається випадкова величина , Визначена на всій осі і має щільність

(1.5)

де - Числові параметри

Будь-які імовірності виду легко обчислюються за допомогою таблиці, в якій наведені значення функції

, Званої зазвичай інтегралом ймовірностей.

Згідно (1.1)

В інтегралі зробимо заміну змінної , Тоді отримаємо

,

де Звідси випливає, що Також

Нормальні випадкові величини дуже часто зустрічаються при дослідженні самих різних за своєю природою питань.

Вибравши , , Знайдемо . Отже,

(1.6)

Ймовірність настільки близька до 1, що іноді останню формулу інтерпретують так: при одному випробуванні практично неможливо отримати значення , Що відрізняється від більше ніж на .

Проводячи велику кількість дослідів, і отримуючи велику кількість випадкових величин можна скористатися центральною граничною теоремою теорії ймовірностей. Ця теорема вперше була сформульована П. Лапласом. Узагальненням цієї теореми займалися багато видатних математики, в тому числі П.Л. Чебишев, А.А. Марков, А.М. Ляпунов. Її доказ досить складно.

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

Позначимо

Суму всіх цих величин позначимо через

Використовуючи співвідношення

отримуємо

Розглянемо тепер нормальну випадкову величину з такими ж параметрами: .

У центральній граничній теоремі стверджується, що для будь-якого інтервалу при великих

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

Використовуючи ці дані з теорії ймовірностей можна перейти до опису загальної схеми методу Монте-Карло. Припустимо, що потрібно обчислити якусь невідому величину . Спробуємо придумати таку випадкову величину , Щоб . Нехай при цьому .

Розглянемо незалежних випадкових величин розподілу яких збігаються з розподілом . Якщо досить велике, то, згідно центральної граничної теореми, розподіл суми буде приблизно нормальним з параметрами . З (1.6) випливає, що .

Останнє співвідношення перепишемо у вигляді:

(1.7)

Це співвідношення дає і метод розрахунку , І оцінку похибки.

Справді, знайдемо значень випадкової величини . З (1.7) видно, що середньоарифметичне цих значень буде наближено одно . З великою ймовірністю похибка наближення не перевершує величини . Ця похибка прямує до нуля із зростанням . На практиці часто використовують не оцінку зверху , А на ймовірну помилку, яка приблизно дорівнює Саме такий зазвичай порядок фактичної похибки розрахунку, яка дорівнює

.

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

Таблиці випадкових чисел використовують переважно при розрахунках вручну. Визначальну роль у цьому відіграють два факти: 1) при використанні ЕОМ легше і зручніше скористатися генератором випадкових чисел, одержуваних тут же, чим завантажувати з пам'яті значення таблиці, яка до того ж, буде займати там місце. 2) При підрахунку вручну немає необхідності використовувати ЕОМ, так як часто необхідно з'ясувати лише порядок шуканої величини.

Генератори випадкових чисел аналізують будь-який процес, доступний для них (шуми в електронних лампах, скачки напруги) і складають послідовність з 0 і 1, з яких складаються числа з певними розрядами, проте такий метод отримання випадкових величин має свої недоліки. По-перше, важко перевірити виробляються числа. Перевірки доводиться робити періодично, так як через якісь неполадки може виникнути так званий дрейф розподілу (нулі і одиниці в якомусь із розрядів стануть з'являтися не однаково часто). По-друге, зазвичай всі розрахунки на ЕОМ проводяться кілька разів, щоб виключити можливість збою. Але відтворити ті ж самі випадкові числа неможливо, якщо їх тільки не запам'ятовувати по ходу рахунку. А якщо запам'ятовувати, то знову з'являється випадок таблиць.

Таким чином, найефективнішим способом отримання випадкових чисел - це використання псевдовипадкових чисел.

Числа, отримані з будь-якої формулою і імітують значення випадкової величини , Називаються псевдовипадковими числами.

Перший алгоритм для отримання псевдовипадкових чисел був запропонований Дж. Нейманом. Він називається методом середини квадратів.

Нехай задано 4-значне число . Зведено його квадрат. Отримаємо 8-значне число . Виберемо 4 середні цифри цього числа і покладемо . Далі і т.д.

Але цей алгоритм не виправдав себе, оскільки виходить занадто багато малих значень. Тому були розроблені інші алгоритми. Найбільшого поширення набув алгоритм, званий методом порівнянь (Д. Лемер): визначається послідовність цілих чисел , В якій початкове число задано, а всі наступні числа обчислюються за однією і тією ж формулою

при (1.8)

За числах обчислюються псевдовипадкові числа

(1.9)

Формула (1.8) означає, що число дорівнює залишку, одержаному при розподілі на , Такий залишок називають найменшим позитивним вирахуванням за модулем Формули (1.8), (1.9) легко реалізувати на ЕОМ.

Переваги методу псевдовипадкових чисел досить очевидні. По-перше, на отримання кожного числа витрачається всього кілька простих операцій, так що швидкість генерування випадкових чисел має той же порядок, що і швидкість роботи ЕОМ. По-друге, програма займає не так багато місця в пам'яті. По-третє, будь-яке з чисел може бути легко відтворений. По-четверте, необхідно лише один раз перевірити «якість» такій послідовності, потім її можна багато разів безбоязно використовувати при розрахунку однотипних завдань.

Єдиний недолік методу - обмеженість кількості псевдовипадкових чисел, тому що якщо послідовність чисел обчислюється на ЕОМ за формулою виду

то ця послідовність обов'язково періодична. Втім, для найбільш поширених псевдовипадкових чисел період настільки великий, що перевершує будь-які практичні потреби. Переважна більшість розрахунків за методом Монте-Карло здійснюється з використанням псевдовипадкових чисел.

Значення будь випадкової величини можна отримати шляхом перетворення значень однієї будь-якої випадкової величини. Зазвичай роль такої випадкової величини грає випадкова величина , Рівномірно розподілена в . Процес знаходження значення якої-небудь випадкової величини шляхом перетворення одного або кількох значень називається розігруванням випадкової величини .

Припустимо, що необхідно отримувати значення випадкової величини , Розподіленої в інтервалі , З щільністю . Доведемо, що значення можна знаходити з рівняння

(1.10)

тобто вибравши чергове значення , Треба вирішити рівняння (1.10) і знайти чергове значення .

Для доказу розглянемо функцію

.

Із загальних властивостей щільності (1.2), (1.3) випливає, що

Значить, функція монотонно зростає від 0 до 1, і будь-яка пряма , Де , Перетинає графік в одній єдиній точці, абсциссу якої ми і приймаємо за . Таким чином, рівняння (1.10) завжди має одне й тільки одне рішення.

Виберемо тепер довільний інтервал , Що міститься усередині . Точкам цього інтервалу відповідають ординати кривої , Що задовольняють нерівності .

Тому, якщо належить інтервалу , То належить інтервалу , І навпаки. Значить

Так як рівномірно розподілена в , То

,

отже

,

а це і означає, що випадкова величина , Що є коренем рівняння (1.10), має щільність ймовірностей .

Може виявитися, що вирішити рівняння (1.10) відносно важко, наприклад, у випадках, коли інтеграл від не виражається через елементарні функції або коли щільність задана графічно. Припустимо, що випадкова величина визначена на кінцевому інтервалі і щільність її обмежена .

Розігрувати значення можна таким чином:

1) вибираються два значення і випадкової величини і будується випадкова точка з координатами

2) якщо точка лежить під кривою , То вважаємо , Якщо ж точка лежить над кривою , То пара відкидається і вибирається нове значення.

1.2 Обчислення інтегралів

Розглянемо функцію , Задану на інтервалі , Потрібно приблизно обчислити інтеграл

(2.1)

Цей інтеграл може бути невласним, але абсолютно збіжним.

Виберемо довільну щільність розподілу , Визначену на інтервалі . Поряд з випадковою величиною , Визначеної в інтервалі з щільністю , Необхідно визначити випадкову величину

Згідно співвідношенню одержимо

Розглянемо тепер однакових незалежних випадкових величин і застосуємо до їх суми центральну граничну теорему. Формула (1.7) в цьому випадку запишеться так:

Останнє співвідношення означає, що якщо вибирати значень , То при досить великому

(2.2)

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

Для розрахунку інтеграла (2.1) можна використовувати будь-яку випадкову величину . Певну в інтервалі з щільністю . У будь-якому випадку . Однак дисперсія , А з нею і оцінка похибки формули (2.2) залежать від того, яка величина використовується, так як

(2.3)

Доведемо, що цей вираз буде мінімальним тоді, коли пропорційна .

Для цього скористаємося нерівністю

, В яких покладемо , . Отримаємо нерівність

(2.4)

З (2.3), (2.4) випливає, що

(2.5)

Залишається довести, що нижня межа дисперсії (2.5) реалізується при виборі щільності . Так як

.

Отже,

,

і права частина (2.3) звертається в праву частину (2.5)

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

Звичайно, вибирати дуже складні не можна, оскільки процедури розігрування стане дуже трудомісткою. Оцінку (2.2) з щільністю , Подібної , Називають істотною вибіркою.

Також якщо стоїть завдання обчислити інтеграл (2.1), перетворимо його до виду

(2.6)

Якщо тепер позначити (2.7)

Те інтеграл приймає вид

(2.8)

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

В окремому випадку, якщо і кінцеві або їх можна вважати кінцевими приблизно, як доцільно вибрати рівномірний закон розподілу.

Як відомо, щільність ймовірності рівномірного закону розподілу в інтервалі дорівнює:

(2.9)

Підставимо в інтеграл (2.6) значення з формули (2.9) і отримаємо:

(2.10)

і розглянемо процедуру обчислення:

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

(2.11)

функції на інтервалі

Таким чином, величина інтеграла (2.10) може бути представлена ​​у вигляді такої формули

(2.12)

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

1.3 Обчислення кратних інтегралів

Зазвичай при обчисленні кратних інтегралів методом Монте-Карло використовують один з двох способів.

Перший спосіб.

Нехай потрібно обчислити кратний інтеграл

(3.1)

по області G, що лежить в мірному одиничному кубі

Виберемо рівномірно розподілених на відрізку послідовностей випадкових чисел

Тоді точки можна розглядати як випадкові, рівномірно розподілені в мірному одиничному кубі.

Нехай із загального числа випадкових точок точок потрапили в область G, інші опинилися поза G. Тоді при достатньо великому має місце наближена формула:

(3.2)

де під розуміється мірний обсяг області інтегрування. Якщо обчислення обсягу скрутно, то можна прийняти , І для наближеного обчислення інтеграла отримаємо:

(3.3)

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

Другий спосіб.

Якщо функція , То інтеграл (3.1) можна розглядати як об'єм тіла в мірному просторі, тобто

(3.5)

де область інтегрування визначається умовами

Якщо в області , То ввівши нову змінну , Одержимо

де область лежить в одиничному мірному кубі

Візьмемо рівномірно розподілених на відрізку випадкових послідовностей

Складемо відповідну послідовність випадкових точок

Нехай із загального числа випадкових точок точок належать обсягу , Тоді має місце наближена формула

(3.6)

2. Практична частина

2.1 Приклад 1

Обчислимо наближено інтеграл

Точне значення його відомо:

Використовуємо для обчислення дві різні випадкові величини , З постійною щільністю (Тобто рівномірна розподілена в інтервалі ) І з лінійною щільністю . Лінійна щільність більше відповідає рекомендації про пропорційність і . Тому слід очікувати, що другий спосіб обчислення дасть кращий результат.

1) Нехай , Формула для розігрування має вигляд . А формула (2.2) прийме вигляд .

Нехай . Як значення використовуємо трійки чисел з табл. 1 (див. додаток), помножені на 0.001. Проміжні результати зведені в табл. 2.1. Результат розрахунку

Таблиця 2.1