SexOnR–04_ModelTypes. Типова побудова моделей у курсі (R-довідник II)

 

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

 

4.1 Концептуальні та реалізовані імітаційні моделі

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

  • Що саме слід спостерігати та досліджувати за допомогою моделі?

  • Як мають бути побудовані перерахунки розрахункових величин?

    • Що є входом для моделі?

    • З яким набором розрахункових величин працює модель?

    • За якими правилами здійснюються перерахунки протягом кожного циклу роботи моделі?

  • Як саме буде використано модель?

    • Які умови роботи моделі слід обирати?

    • Що є вихідними даними моделі?

    • Як будуть презентуватися вихідні дані?

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

“Ідеальна” концептуальна модель мала б бути незалежною від будь-якого засобу її реалізації, наприклад, R, Python, Excel або Calc чи, припустимо, перекладання вручну картонок з надписами зі стопочки у стопочку чи переписування записів олівцем на таких картках у залежності від кидання гральних кубиків (моделі, які ми будемо розглядати, цілком можна реалізувати і таким чином; єдиний недолік — це незручно). Втім, деякі з рішень, які слід приймати при побудові концептуальної моделі, звісно, залежать від особливостей засобу реалізації. Якщо б ми насправді розробляли модель, що реалізується перекладанням карток, це б відбилося на її задумі. Звісно, у подальшому викладенні ми будемо обговорвати моделі, “заточені” на реалізацію засобами R.

 

4.2 Типова структура моделі

Загалом, імітаційні моделі, які ми будемо використовувати, складаються з подібного набору функціональних блоків. Наведемо їх можливий перелік “на виріст”; з врахуванням моделей, які в цьому курсі якщо ми й будемо будувати, то лише у найкращому разі наприкінці.

 

  • I. ENTRANCE — ВХІД:

    • I.1. Initial script commands — початкові команди скрипту (необхідні бібліотеки, домашня директорія, читання збережених файлів, початкове число для генерації псевдовипадкових чисел, чистка середовища тощо);

    • I.2. Initial state of the system — початковий стан модельованої системи;

    • I.3. Parameters — параметри, за якими відбуваються перетворення;

    • I.4. Experimental conditions — умови експерименту з моделлю (кількість циклів, ітерацій тощо);

    • I.5. The changeable parameters combinations — комбінації змінюваних параметрів (в моделях ІІІ типу).

  • II. TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:

    • II.1. Users function setting — створення користувацьких функцій, якщо вони необхідні для подальших розрахунків;

    • II.2. Parameters combinations mechanism — механізм комбінації змінюваних параметрів (в моделях ІІІ типу);

    • II.3. Objects creation — створення об’єктів, що містять розрахункові величини, накопичують результати розрахунків для подальшого аналізу тощо; якщо змінювані параметри не відбиваються на створюваних об’єктах, може розташовуватися перед попереднім пунктом;

    • II.4. Calculation rules data — Дані з правилами перерахунків (наприклад, таблиці з результатами схрещувань, що будуть враховуватися у розрахунках тощо).

  • III. CALCULATIONS — РОЗРАХУНКИ:

    • III.1. Higher level cycles running — запуск циклів вищого рівня: циклів ітерацій у моделях II типу тощо; за необхідності — створення об'єктів, що потрібно створювати у кожній ітерації:

    • III.2. Initial composition creation — утворення початкового складу розрахункових величин залежно від вхідних даних;

    • III.3. Main work cycle — основний робочий цикл, що у типовому випадку достатньо складних моделей організовано як перерахунки певних чисельностей за фазами, що позначені грецькими літерами: αtN → βtN → γtN → δtN → … → ωtN → αt+1N…; ці перерахунки залежать від:

      • параметрів;

      • умов експерименту;

      • значень інших розрахункових величин;

      • псевдовипадкових, генерованих R за певним алгоритмом чисел;

    • III.4. Ending of higher level cycles — завершення циклів вищого рівня, за необхідності — з умовами дострокового припинення циклів.

  • IV. FINISHING — ЗАВЕРШЕННЯ:

    • IV.1. Creating objects for results integration — створення об’єктів, що інтегрують результати (за такої необхідності, наприклад, для візуалізації);

    • IV.2. Results saving — збереження об’єктів з результатами;

    • IV.3. Results viewing — перегляд підсумків — основних чи додаткових об’єктів або їх певних фрагментів;

    • IV.4. Results visualization — візуалізація підсумків; за необхідності — збереження результатів візуалізації.

 

4.3 Приклад найпростішої моделі: експоненційне зростання

4.3.1 Сутність і формули

Серед моделей, що описують популяційні процеси, найстарішою (з історію, що йде від XIII сторіччя, від Леонардо Фібоначчі) є модель експоненційного зростання; за необхідності, можна отримати додаткову інформацію про таке зростання у підручнику “Екологія: біологія взаємодії” Д. Шабанова та М. Кравченко у розділі “4.04. Експоненційне і логістичне зростання чисельності популяції”.

Побудуємо таку модель. Класичне рівняння приросту під час експоненційного зростання є таким:

dN/dt = r × N, де N — чисельність популяції, t — час, а r — параметр Мальтуса, що визначає здатність популяції до зростання. Імітаційне моделювання пов’язане з поділом модельованого процесу на окремі кроки. Наведеному диференційному рівнянню для розрахунку приросту відповідає таке різницеве рівняння:

t+1N = tN × (1 + r). Використаємо це рівняння.

 

4.3.2 Реалізація в R

Ми використаємо далеко не усі елементи структури моделі, порівняно з тими, що перелічені у пункті 4.2. Зрозуміло, що це є наслідком відносної простоти обговорюваної моделі.

# Експоненційне зростання
# I. ENTRANCE — ВХІД:
# Initial state of the system — початковий стан модельованої системи
al_1_N  <-  10 # Початкова чисельність модельної популяції

# Parameters — параметри
r <- 0.05 # Параметр Мальтуса, приріст чисельності на кожному кроці

# Experimental conditions — умови експерименту з моделлю
cycles <- 50 # I.3 Кількість циклів

# II TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:
# Objects creation — створення об’єктів
N <- rep(NA, cycles+1) # Створення вектору для розрахунків та збору результатів (у цьому випадку ці дві задачі вирішує тій самий вектор)

# III. CALCULATIONS — РОЗРАХУНКИ:
# Initial composition creation — утворення початкового складу
N[1] <- al_1_N # Перше значення визначається відповідно до початкового стану

# Main work cycle — основний робочий цикл
for (t in 1:cycles) { N[t+1] <- floor(N[t] * (1 + r) + runif(1)) } # Робочий цикл моделі повторює робочу формулу таку кількість разів, яка визначена умовами експерименту. Результат перерахунку залежить від значення розрахункової величини (попереднього значення в векторі), параметру (що визначає швидкість зростання популяції), умов експерименту (кількості кроків, яку відраховує лічильник у циклі) та випадкової величини(визначає ймовірнісне скорочення отриманої величини до цілих значень)

# IV. FINISHING — ЗАВЕРШЕННЯ:
# Results viewing — перегляд підсумків
N # Результати збережені в векторі N; на них можна просто подивитися
##  [1] 10 10 10 10 11 11 11 11 12 13 13 14 14 15 15 16 17 18 19 20 21 22 23 24 26
## [26] 27 29 30 31 33 35 37 39 41 43 45 48 50 52 55 57 60 63 67 71 75 79 83 87 91
## [51] 96
# Results visualization — візуалізація підсумків
plot(N) # Ця команда виводить результати на найпростішій графік

Рис. 4.3.2.1 Візуалізація моделі експоненційого зростання

 

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

 

4.3.3 Ймовірнісне округлення

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

Слід використовувати округлення. Яке? Стандартний спосіб округлювання не завжди підходить.

round(3.3)
## [1] 3
round(3.4999)
## [1] 3
round(3.5)
## [1] 4

Ви бачите, що існує певне межа, що чітко розподіляє округлювані значення. Чи має це біологічний сенс?

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

ProbabRound <- rep(NA, 1000)
for (i in 1:1000) {ProbabRound[i] <- floor(3.3 + runif(1))}
mean(ProbabRound)
## [1] 3.324

За допомогою ймовірнісного округлення 3.3 у трьох випадках з 10 округлиться до 4, а у сьома випадках з десяти — до 3. У прикладі вище ми розрахували середнє значення 1000 таких спроб. Воно відрізняється від одного перерахунку до іншого, але зазвичай є достатньо близьким до 3.3.

У переважній кількості випадків в наших моделях буде використовуватися саме ймовірнісне округлення.

 

4.3.4 Схема організації розглянутої моделі (експоненціального зростання)

Щоб краще зрозуміти логіку роботи описаної моделі, її можна представити у вигляді схеми.

 

Рис. 4.3.4.1 У графічній формі структуру описаної найпростішої моделі можна представити так

 

Спробуйте співставити текстовий опис типової структури моделі (пункт 4.2), приклад найпростішої моделі в R (пункт 4.3.2) і графічну схему!

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

 

4.4 Модель І типу: визначення динаміки процесу з певної кількості циклів, що складаються з етапів

4.4.1 Побудова концептуальної моделі та система позначень етапів у циклі

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

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

Логістичне рівняння у різницевий формі виглядає так:

t+1N = tN × (1 + r) × (K–tN/K), де N — чисельність популяції, t — час, r — параметр Мальтуса, що визначає здатність популяції до зростання, а K — параметр Ферхюльста, що задає обмеження чисельності популяції.

Це рівняння можна було б відбити і у моделі. яку ми побудували у пункті 4.3, достатньо було б додати до параметрів K та дещо ускладнити формулу робочого циклу. Але ми побудуємо більш складну модель. Давайте врахуємо в ній тривалість життя особин (позначимо її maxa, або max_a у тексті R-моделі).

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

Нам знадобиться певна система позначень. Будемо позначати етапи в циклі за допомогою грецьких літер, а самі цикли — цифрами. Внесемо певну зміну у позначення порівняно з моделлю експоненційного зростання. В тій моделі ми характеризували модельну популяцію на циклі t як tN (маючи на увазі N — number). Якщо модельна популяція складається з ідентичних (з точки зору моделювання) особин, її можна характеризувати просто чисельністю. Якщо ми будемо враховувати вік особин, модельна популяція стає для нас неоднорідною; будемо позначати її tP (P — population).

Виділимо у циклі чотири етапи:

  • Альфа-популяція, αP, сукупність особин на початку циклу; на першому циклі залежить від початкового складу, заданого у вхідних параметрах: α1P = al_1_N, а надалі визначається омега-популяцією попереднього циклу: αt+1P = ωtP;

  • Бета-популяція, βP, сукупність особин після загибелі особин максимального віку, βt+max_aP = αt+max_aP – γtP;

  • Гама-(суб)популяція, γP, приплід: γtP = βtP × (1 + r) × (K –βtP) / K;

  • Омега-популяція, ωP, прикінцева чисельність у циклі: ωtP = βtP + γtP.

 

4.4.2 Система позначень етапів у циклі

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

Тепер, коли ми обговорили позначення етапів у циклі, ви зрозумієте, чому у моделі експоненційного зростання (пункт 4.3) ми позначили початкову чисельність як al_1_N. Крім іншого, слід зазначити, що імена об’єктів в R не можуть починатися з цифр, і саме тому у наших позначеннях ми будемо спочатку вказувати етап, а потім — цикл.

У деяких моделях ми будемо використовувати більшу кількість грецьких літер (наведемо їх перелік; у дужках будемо вказувати позначення з двох літер латиницею, які ми будемо використовувати в R-моделях; напівжирним виділені позначення, у яких можна заплутатися).

α — альфа (al); β — бета (be); γ — гамма (ga); δ — дельта (de); ε — епсилон (ep); ζ — дзета (ze); η — ета (et); θ — тета (th); ι — йота (io); κ — каппа (ka); λ — лямбда (la); μ — мю (mu); ν — ню (nu); ξ — ксі (xi); ο — омікрон (oc); π — пі (pi); ρ — ро (rh); ς — сигма (si); τ — тау (ta); υ — іпсилон (up); φ — фі (ph); χ — хі (ch); ψ — псі (ps); ω — омега (om).

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

“Я Альфа й Омега, Перший і Останній, Початок і Кінець”.
Об’явлення св. Івана Богослова, 22.13. Біблія або Книги Святого Письма Старого і Нового Заповіту. Переклад Івана Огієнка

 

4.4.3 R-модель логістичного зростання з врахуванням смертності у певному віці

# Логістичне зростання з врахуванням смертності у певному віці
# I. ENTRANCE — ВХІД:
# Initial state of the system — початковий стан модельованої системи
al_1_N  <-  10 # Початкова чисельність модельної популяції
# Parameters — параметри
r <- 0.3 # Параметр Мальтуса, приріст чисельності на кожному кроці
K <- 500 # Параметр Ферхюльста, що задає обмеження чисельності популяції
max_a <- 5 # Тривалість життя особини (у циклах)
# Experimental conditions — умови експерименту з моделлю
cycles <- 100 # Кількість циклів
# TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:
# Objects creation — створення об’єктів
al_P <- rep(NA, cycles+1) # 
ga_P <- rep(NA, cycles) #
# III. CALCULATIONS — РОЗРАХУНКИ:
# Initial composition creation — утворення початкового складу
al_P[1] <- al_1_N # Перше значення визначається відповідно до початкового стану
# Main work cycle — основний робочий цикл
for (t in 1:cycles) { # Початок робочого циклу
  Way <- (t!=max_a) + (t > max_a) +1 # Величина, що відбиває співвідношення t і max_a
  switch(Way, # Перемикач, що обирає один з трьох варіантів залежно від значення Way
         be_P <- al_P[t] - al_P[1], # У разі, якщо t == max_a
         be_P <- al_P[t], # У разі, якщо t < max_a
         be_P <- be_P <- al_P[t] - ga_P[t-max_a]) # У разі, якщо t > max_a
  ga_P[t] <- floor(be_P * r * (K - be_P) / K + runif(1)) # Приплід
  om_P <- be_P + ga_P[t] # Прикінцевий склад
  al_P[t+1] <- om_P # Перехід на наступний цикл роботи моделі
                    } # Кінець робочого циклу
# IV. FINISHING — ЗАВЕРШЕННЯ:
# Results viewing — перегляд підсумків
al_P # Результати збережені в векторі al_P; на них можна просто подивитися
##   [1]  10  13  17  22  28  24  27  30  32  33  35  37  39  41  43  45  48  50
##  [19]  52  54  56  57  59  61  63  64  66  67  68  70  72  73  74  75  76  77
##  [37]  78  80  82  84  85  87  88  89  90  92  93  94  96  97  97  97  97  96
##  [55]  95  95  95  96  97  97  97  98  98  98  99 100 100 100 100 101 101 102
##  [73] 102 103 102 102 101 102 101 102 102 103 102 102 102 103 102 103 104 104
##  [91] 104 105 105 105 105 105 105 105 104 104 103
# Results visualization — візуалізація підсумків
plot(al_P) # Ця команда виводить результати на найпростішій графік

Рис. 4.4.3.1 Динаміка моделі I типу

 

Спочатку обговоримо технічні деталі, пов’язані з тим, як побудована модель. У ряді αP ⟶ βP ⟶ γP ⟶ ωP нема сенсу зберігати усі значення. Усі альфа-чисельності популяції необхідно зберігати, тому що на їх підставі буде побудовано графік. Гама-чисельності слід запам’ятовувати, тому що до їх попередніх значень доводиться повертатися, щоб імітувати смертність старих особин. А ось бета- і омега-чисельності можна записувати не для усієї послідовності, а перераховувати їх у кожному циклі. Саме тому на етапі II.1 створюються два вектори, а не чотири.

Ще одна складність полягає в тому, що чисельність особин, які лишаються після загибелі старих представників, може розраховуватися трьома різними способами залежно від того, як співвідносяться t (номер циклу) та max_a (максимальна тривалість життя). Щоб обрати формулу, використано функцію switch(). Побудована вона так: switch(об’єкт згідно зі значенням якого здійснюється вибір, вибір при значенні 1, вибір при значенні 2, …).

Об’єкт, який керує цим перемикачем, має назву Way (це вектор одиничної довжини). Щоб його розрахувати, використана сума, яка включає логічні вирази. Так, вираз (t!=max_a) (тобто t не дорівнює max_a) у арифметичних виразах може приймати значення 1 (TRUE) або 0 (FALSE). Роздивіться, як використана функція switch(). Розгляньте три випадки: t = max_a, t < max_a та t > max_a. Розрахуйте, яке значення буде мати Way і яка формула буде використана для розрахунку be_P.

Якщо ви зрозуміли, як працює ця R-модель, можна замислитися над цікавими питаннями, що стосуються її властивостей.

Що буде, якщо тривалість життя особин зменшиться — хоча б до 4 років замість 5? Чому?

Ємність середовища у наведеному прикладі — 500 особин. Чому ж зростання чисельності модельної популяції зупиняється на рівні, близькому до 100?

Що є причиною коливань, які можна побачити на графіку?

 

4.4.4 Схема організації розглянутої моделі І типу (логістичного зростання)

Як і в попередньому випадку, спробуємо показати, як організована наша модель, за допомогою схеми.

Головна новизна цієї схеми, у порівнянні з попередньою, у тому, що вона показує організацію перерахунків як повторення певних циклів, що складаються з етапів.

Рис. 4.4.4.1 Схема, що показує структуру описаної моделі І типу

 

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

4.5.1 Призначення моделі ІІ типу

Як ви знаєте, моделі можуть бути детермінованими та недетермінованими (такими, що включають стохастичні, випадкові елементи). У тому разі, коли ми можемо чітко і впевнено встановити усі чинники, що впливають на очікуваний результат, ми можемо використовувати детерміновані моделі. Більш поширеними є випадки, коли і повний набір чинників нам невідомий, і точне визначення сили цих чинників для нас є неможливим…

Рис. 4.5.1.1 Фрагмент презентації до іншого курсу: два типи прогнозів

 

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

Припустимо, за допомогою моделі І типу ми встановили якусь динаміку. І що? Наскільки ми можемо довіряти результату? Чи є він закономірним наслідком дії чинників, які ми врахували, чи відбиває якісь малоймовірний випадок? Щоб встановити це, моделювання краще повторити, і не один раз. Повтори з тими самими параметрами характерні для моделей ІІ типу, а зі змінними — для моделей ІІІ типу.

Звісно, щоб узнати, чи завжди ми отримуємо однаковий результат, достатньо і моделі І типу. Ми можемо просто раз за разом запускати її, та занотовувати те, що вийшло. Але навіщо робити це “вручну”, коли це можна зробити за допомогою невеликої перебудови самої моделі?

Завдяки переходу від моделі І типу до ІІ типу ми перейдемо від одиничного спостереження результату модельного експерименту до можливості проводити аналіз розподілу ймовірностей результатів.

 

4.5.2 Перебудова моделі І типу у модель ІІ типу

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

Перелічені зміни зроблені в наведеному скрипті. На даній html-сторінці ми виділимо внесені зміни цегляним кольором (звісно, при роботі з R в RStudio це неможливо)

# Логістичне зростання з врахуванням смертності у певному віці (Модель ІІ типу)
# I. ENTRANCE — ВХІД:
# Initial state of the system — початковий стан модельованої системи
al_1_N  <-  5 # Початкова чисельність модельної популяції

# Parameters — параметри
r <- 0.3 # Параметр Мальтуса, приріст чисельності на кожному кроці
K <- 500 # Параметр Ферхюльста, що задає обмеження чисельності популяції
max_a <- 5 # Тривалість життя особини (у циклах)

# Experimental conditions — умови експерименту з моделлю
cycles <- 100 # Кількість циклів
iterat <- 10 # Кількість ітерацій

# II. TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:
# Objects creation — створення об’єктів
al_P <- rep(NA, cycles+1) 
ga_P <- rep(NA, cycles)
Rezult <- matrix(NA, nrow = cycles+1, ncol = iterat) # У разі використання ітерацій потрібна матриця, а не вектор, як у моделі І типу

# III. CALCULATIONS — РОЗРАХУНКИ:
# Higher level cycles running — запуск циклів вищого рівня
for (i in 1:iterat) { # Початок ітерації

# Initial composition creation — утворення початкового складу
al_P[1] <- al_1_N # Перше значення визначається відповідно до початкового стану

# Main work cycle — основний робочий цикл
    Rezult[1, i] <- al_1_N
    for (t in 1:cycles) { # Початок робочого циклу
        Way <- (t!=max_a) + (t > max_a) +1 # Величина, що відбиває співвідношення t і max_a
        switch(Way, # Перемикач, що обирає один з трьох варіантів залежно від значення Way
           be_P <- al_P[t] - al_P[1], # У разі, якщо t == max_a
           be_P <- al_P[t], # У разі, якщо t < max_a
           be_P <- be_P <- al_P[t] - ga_P[t-max_a]) # У разі, якщо t > max_a
        ga_P[t] <- floor(be_P * r * (K - be_P) / K + runif(1)) # Приплід
        om_P <- be_P + ga_P[t] # Прикінцевий склад
        al_P[t+1] <- om_P # Перехід на наступний цикл роботи моделі
        Rezult[t+1, i] <- al_P[t+1]
                         } # Кінець робочого циклу

# Ending of higher level cycles — завершення циклів вищого рівня   
                     } # Кінець циклу ітерацій

# IV. FINISHING — ЗАВЕРШЕННЯ:
# Results viewing — перегляд підсумків
head(Rezult) # Підсумки збережені в цій матриці; оскільки вона велика, можна подивитися на її початок...
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]    5    5    5    5    5    5    5    5    5     5
## [2,]    6    7    6    6    7    6    7    6    6     7
## [3,]    8    9    8    7    9    8    9    7    8     9
## [4,]   11   11   11    9   11   11   11    9   11    12
## [5,]   14   14   14   11   14   14   14   11   15    15
## [6,]   12   12   12    8   12   11   12    8   13    13
tail(Rezult) # ...та кінець
##        [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [96,]  105  105  105  105  105  104  105  105  103   103
##  [97,]  105  105  105  105  105  103  105  105  104   104
##  [98,]  105  105  105  105  105  103  105  105  104   104
##  [99,]  105  105  105  105  105  103  105  105  105   104
## [100,]  105  105  105  105  105  102  105  105  105   105
## [101,]  105  105  105  105  105  102  105  105  105   105
# Results visualization — візуалізація підсумків
plot(Rezult[, 1], type="l", lty=1, col="red", ylim=c(0, 110),
     main = "Динаміка чисельності модельної  популяції\n у 10 різних ітераціях",
     xlab="Цикли імітації", ylab="Чисельність модельної популяції")
lines(Rezult[, 2], type="l", lty=1, col="blue")
lines(Rezult[, 3], type="l", lty=1, col="brown")
lines(Rezult[, 4], type="l", lty=1, col="chartreuse")
lines(Rezult[, 5], type="l", lty=1, col="coral")
lines(Rezult[, 6], type="l", lty=1, col="darkviolet")
lines(Rezult[, 7], type="l", lty=1, col="gold")
lines(Rezult[, 8], type="l", lty=1, col="deepskyblue")
lines(Rezult[, 9], type="l", lty=1, col="hotpink")
lines(Rezult[, 10], type="l", lty=1, col="lightsalmon")

Рис. 4.5.2.1 Результат роботи моделі II типу. Ми бачимо 10 окремих імітацій

 

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

 

4.5.3 Схема організації розглянутої моделі ІІ типу

Як і в попередніх випадках, наведемо схему, де показано, як працює обговорена модель. Ми бачимо “стопку” з ітерацій, кожна з яких відповідає моделі І типу. На кожну з них впливає комплекс вхідних значень, параметрів та умов. На підставі співставлення виходів з кожної ітерації здійснюється аналіз розподілу ймовірностей результатів моделювання за даних параметрів.

Рис. 4.5.3.1 Схема, що показує структуру описаної моделі ІІ типу

 

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

 

4.6 Модель ІІІ типу: встановлення впливу параметрів на результат моделювання

4.6.1 Призначення моделі ІІІ типу

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

  • чи може бути спостережувана поведінка (динаміка, ознаки) природних систем наслідком дії певного механізму?

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

  • у якому разі певний набір чинників може викликати ті феномени, що ми їх спостерігаємо?

  • як зміни певного чинника впливають на поведінку досліджуваної системи, у разі, якщо наші припущення про механізми її функціювання є правильними?

Для відповіді на перелічені питання використання моделі ІІ типу є недостатнім. Припустимо, ми визначили розподіл ймовірностей вихідних станів модельної системи за визначеного набору вхідних параметрів. Але чи можемо ми без перевірки бути впевнені, що за інших значень система буде вести себе подібним чином? Врахуйте, що моделі ми використовуємо саме у таких випадках, коли поведінка досліджуваної системи не є тривіальною, коли для її дослідження треба врахувати взаємодію між певним набором чинників.

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

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

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

 

4.6.2 Перебудова моделі ІІ типу у модель ІІІ типу

В моделі ІІІ типу нам доведеться “натягнути” на ті два цикли, що вже є у моделі ІІ типу, ще два: для кожного зі змінюваних параметрів.

Загалом, у порівнянні з моделлю ІІ типу необхідними є перелічені зміни (ми будемо змінювати параметр Мальтуса та максимальний вік особин):

  • для візуалізації результатів буде використана додаткова бібліотека для тривимірної графіки plot3D; у такому разі, на початку роботи моделі цю бібліотеку слід активувати, а якщо вона не інстальована у системі — попередньо інсталювати (у наведеному прикладі команда інсталювання вимкнена “ґраткою”);

  • команди, що задають постійні значення параметрів, що будуть змінюватися, слід прибрати або вимкнути “ґраткою”;

  • слід задати вектори (vector_r та vector_max_a), що будуть містити значення змінюваних параметрів, які будуть комбінуватися;

  • оскільки матриця Rezult, у яку збиралися підсумки ітерацій, працює лише з певним сполученням змінюваних параметрів, слід створити об’єкт, де буде збиратися інтегральна величина, яка буде характеризувати результат кожної групи ітерацій; назвемо таку матрицю Outcome;

  • у якості інтегральної величини можна використовувати середнє значення чисельності модельної популяції наприкінці кожної ітерації, коли коливання чисельності стають випадковими; команда Outcome[a, b] <- mean(Rezult[(cycles-20):(cycles+1), ]) забезпечить запис середнього значення 21 останньо стану чисельності модельної популяції в кожній ітерації;

  • візуалізація результатів з використанням функції persp3D дозволить побачити загальний результат.

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

# Логістичне зростання з врахуванням смертності у певному віці (Модель ІІІ типу)
# I. ENTRANCE — ВХІД:
# Initial script commands — початкові команди скрипту
# install.packages("plot3D")
library(plot3D)

# Initial state of the system — початковий стан модельованої системи
al_1_N  <-  5 # Початкова чисельність модельної популяції

# Parameters — параметри
# r <- 0.3 # Параметр Мальтуса, приріст чисельності на кожному кроці; ЗМІНЮВАНИЙ ПАРАМЕТР 1
K <- 500 # Параметр Ферхюльста, що задає обмеження чисельності популяції
# max_a <- 5 # Тривалість життя особини (у циклах); ЗМІНЮВАНИЙ ПАРАМЕТР 2

# Experimental conditions — умови експерименту з моделлю
cycles <- 150 # Кількість циклів
iterat <- 20 # Кількість ітерацій

# The changeable parameters combinations — комбінації змінюваних параметрів
# Створення векторів зі значеннями змінних параметрів, які будуть перебиратися у ході роботи:
vector_r <- c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6) # Варіанти значень параметру Мальтуса
vector_max_a <- c(2, 3, 4, 5, 6, 7) # Варіанти максимальної тривалості життя

# II TRANSFORMATIONS SYSTEM CREATION — СТВОРЕННЯ СИСТЕМИ ПЕРЕТВОРЕНЬ:
# Parameters combinations mechanism — механізм комбінації змінюваних параметрів
Outcome <- matrix(NA, nrow = length(vector_r), ncol = length(vector_max_a), dimnames = list(vector_r, vector_max_a))

# Цикли, що перебирають досліджувані сполучення параметрів:
for (a in 1:length(vector_r)){ # Цикл, що змінює значення параметру Мальтуса
    r <- vector_r[a] # Змінюване значення чисельності популяції
    for (b in 1:length(vector_max_a)){ # Цикл, що змінює максимальну тривалість життя
        max_a <- vector_max_a[b] # Змінюване значення сили добору

# Objects creation — створення об’єктів
al_P <- rep(NA, cycles+1) # 
ga_P <- rep(NA, cycles) #
Rezult <- matrix(NA, nrow = cycles+1, ncol = iterat) # У разі використання ітерацій, потрібна матриця, а не вектор, як у моделі І типу

# III. CALCULATIONS — РОЗРАХУНКИ:
# Higher level cycles running — запуск циклів вищого рівня
for (i in 1:iterat) { # Початок ітерації

# Initial composition creation — утворення початкового складу
  al_P[1] <- al_1_N # Перше значення визначається відповідно до початкового стану
  Rezult[1, i] <- al_1_N

# Main work cycle — основний робочий цикл  
  for (t in 1:cycles) { # Початок робочого циклу
      Way <- (t!=max_a) + (t > max_a) +1 # Величина, що відбиває співвідношення t і max_a
      switch(Way, # Перемикач, що обирає один з трьох варіантів залежно від значення Way
         be_P <- al_P[t] - al_P[1], # У разі, якщо t == max_a
         be_P <- al_P[t], # У разі, якщо t < max_a
         be_P <- be_P <- al_P[t] - ga_P[t-max_a]) # У разі, якщо t > max_a
      ga_P[t] <- floor(be_P * r * (K - be_P) / K + runif(1)) # Приплід
      om_P <- be_P + ga_P[t] # Прикінцевий склад
      al_P[t+1] <- om_P # Перехід на наступний цикл роботи моделі
      Rezult[t+1, i] <- al_P[t+1]
                         } # Кінець робочого циклу

# Ending of higher level cycles — завершення циклів вищого рівня  
                     } # Кінець циклу ітерацій

Outcome[a, b] <- mean(Rezult[(cycles-20):(cycles+1), ]) # Запис певного значення у матрицю зі сполученнями параметрів 
                                      } # Кінець циклу, що змінює максимальну тривалість життя
                              } # Кінець циклу, що змінює значення параметру Мальтуса

# IV. FINISHING — ЗАВЕРШЕННЯ:
# Results saving — збереження об’єктів з результатами
# save(Outcome, file = "Outcome.RData")

# Results visualization — візуалізація підсумків
persp3D(z = Outcome, x = vector_r, y = vector_max_a, theta = 330, phi = 30, # Побудова графіку
       xlab = "r", ylab = "max_a",
        expand = 0.8, facets = FALSE, scale = TRUE,  ticktype = "detailed",nticks = 6, cex.axis=0.9, 
        clab = "mean N", colkey = list(side = 2, length = 0.5))

Рис. 4.6.2.1 Результат роботи моделі III типу. Ми бачимо, як чисельність модельної популяції (mean N) залежить від тривалості життя особин (max_a) та від їхнього репродуктивного потенціала, тобто кількості потомків (r)

 

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

 

4.6.3 Схема організації розглянутої моделі ІІІ типу

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

На схемі показана таблиця, “комірки” якої відповідають різним сполученням досліджуваних параметрів. У кожній такій “комірці” — фактично, модель ІІ типу. На виході кожної такої моделі — вже не розподіл ймовірностей, а якась одна величина. Значення такої величини можуть стати, наприклад, окремими ділянками, що формують тривимірну поверхню. Ця поверхня відбиватиме залежність результату моделювання від значень параметрів. Звісно, інформативність, користь для розуміння процесів, у такої моделі набагато вища, ніж у попередніх.

Рис. 4.5.3.1 Схема, що показує структуру описаної моделі ІІІ типу

 

Звісно, не є обов’язковою вимогою, щоб у моделі ІІІ типу досліджувалися зміни саме двох параметрів. Три параметри — забагато; незрозуміло, як візуалізувати результат. Один переметр — цілком можливо, і звичайний графік (на вісі абсцис — значення досліджуваного параметру, на вісі ординат — інтегруча величина, що характеризує результат моделювання) є цілком задовільним варіантом подання результатів такого дослідження.

 

4.7 Деякі особливостей R-моделей

4.7.1 Зовнішній вигляд R-скрипту

Автор даного тексту має признатися, що має дуже погану пам’ять на конкретні деталі, пов’язані з моделями і загалом з R-скриптами. Під час їх написання навіть найпростіші команди R доводиться шукати у різноманітних пам’ятках. Дуже швидко втрачається пам’ять про конкретні прийоми, що були застосовані для вирішення тієї або іншої задачі. Коли через місяць-два доводиться подивитися на власний скрипт R, він сприймається як незнайомий…

Під час написання скрипту доводиться пробувати різні варіанти послідовностей команд. Автор зазвичай лишає їх після основного тексту скрипту. Коли шукаєш можливості, пишеш якомога швидше без усіляких пояснень, з використанням назв об’єктів у одну букву. Знайдено рішення або ні — коли дивишся на ці спроби через деякий час, зрозуміти нічого неможливо. Те, що зберігав “на майбутнє” без пояснень, виявляється непотребом.

Як зменшити руйнівний ефект від такого швидкого забування? Приділяти значну увагу зовнішньому вигляду скрипту. Прикладом того, як це робить автор (коли в нього є час акуратно поставитися до результатів власної роботи) є попередній фрагмент R-коду. На які його особливості варто звернути увагу?

  • Скрипт містить багато коментарів. Деякі коментарі розміщуються після команд R у одному з ними рядку, деякі — перед робочим рядком (і закінчуються двокрапкою).

  • Скрипт має заголовок.

  • Об’єкти мають виразні імена (Rezult, Fate, Lot) або утворюються за чіткими правилами (Al_Pop, Be_Pop, Ga_Pop, Om_Pop, al_1_N_a, al_1_N_A, Be_N_a, Be_N_A), які відповідають позначенням при описі моделі.

  • Рядки, що належать до певного циклу або певної умови, розташовуються з відступом, що показує залежність від рядків, де починаються ці функції (for, if тощо).

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

  • У тексті скрипту невипадково розставлені пробіли, що полегшує його сприйняття. Пробіли стоять перед та після символів присвоєння <-, знаків = та ґраток #, після ком та ком з точкою. Знаки +, - тощо розміщуються між пробілів, якщо поєднують відносно складні назви (наприклад, (length(which(Al_Pop==1)) + length(which(Al_Pop==0))), та без пробілів — у разі коротких назв (Rezult[cycles+1, i]).

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

 

4.7.2 Початок R-документу

Початок R-документу може не містити якихось додаткових елементів (попередній наведений фрагмент R-коду є цілком працездатним!), а може починатися з кількох корисних команд.

У багатьох випадках R моделі використовують можливості додаткових пакетів. Ці пакети слід один раз встановити на тому комп’ютері, де працює R, та кожного разу активувати при запуску моделі. Для цього використовуються такі команди.

# install.packages("plot3D")
library(plot3D)

Пакет, якій встановлює перша команда та активує друга, знадобиться нам при візуалізації результатів створеної нами моделі. Після того, як цей пакет встановлений, першу команду можна закрити “ґраткою”.

Коли на комп’ютері користувача за допомогою R розв’язується кілька різних задач, результати роботи починають плутатися. Часто зручно розділити їх по окремих теках (директоріях). Зрозуміло, для кожного користувача шлях до робочої директорії є своїм.

setwd("~/!_Courses/Sex_on_R") # Робоча директорія (лише на комп'ютері Д.Ш.!!!)

Спосіб записувати шлях до робочої директорії відрізняється у різних операційних системах (у попередньому прикладі показано варіант, характерний для Linux; у Windows адреса виглядає іншим чином). Щоб встановити, як треба записувати адресу робочої директорії, можна встановити її “вручну” (перейти у меню Files, обрати потрібну теку та виконати команду More / Set As Working Directory). У консолі RStudio, як це показано на рисунку, з’явиться відповідна команда (RStudio “перекладе” вибір, зроблений за допомогою курсору, на мову команд R у коректному запису).

Рис. 4.7.2.1 Після вибору опції Set As Working Directory у консолі RStudio з’явилася відповідна команда

 

При повторних запусках того самого скрипту, або подібних один до одного скриптів з’являється певна небезпека. У робочому просторі програми, в Environment, створюються певні об’єкти. У якомусь наступному запуску скрипту або його фрагменту можна не помітити, що у розрахунках використовується раніше створений об’єкт. Це є одним з можливих джерел помилок. Щоб цього не відбувалося, на початку скрипту можна вставити команду, що очистить Environment — “простір”, у якому виконуються команди. Звісно, можна не використовувати наведену команду, а очистити Environment “вручну”.

rm(list = ls()) # Очищення раніше збережених об'єктів в Environment

Ця лекція написана в R Markdown. Якщо ви працюєте з документом Markdown (він має розширення .Rmd), при його запуску усі R-розрахунки відбуваються заново. Якщо модель є недетермінованою, тобто використовує випадкові числа, ви побачите не ті результати, що бачить автор.

Є й більш важливі випадки, коли бажано зробити результати робот R передбачуваними. Припустимо, ви створили модель. При її запусках вона у деяких випадках видає помилки — а у деяких випадках не видає. Помилки пов’язані з певним сполученням розрахункових величин у моделі… Щоб покроково дослідити цю ситуацію, слід передивитися проміжні результати розрахунків, щоб зрозуміти, де щось пішло не так, як треба. Але з кожним запуском R буде вибирати інші варіанти псевдовипадкових чисел, що використовуються у розрахунках. Як бути?

“Заморозити” початкову точку розрахунків псевдовипадкових чисел. Це можна зробити за допомогою команди set.seed(); у дужки слід ввести якесь число.

# set.seed(1234567)

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

Ще одним важливим доповненням до основного “тіла” програми можуть бути команди збереження та читання результатів тривалих розрахунків. Припустимо, у наших розрахунках довелося тривалий час обчислювати матрицю Rezult. Щоб використати її у подальшому раціонально зберегти її у вигляді окремого файлу. Якщо не вказувати шлях до певної директорії, цей файл буде збережено у робочій директорії, з якою працює скрипт. Після того, як ми створили матрицю Rezult, ми могли б її зберегти.

# save(Rezult, file = "Rezult.RData")

Збережений файл, що містить матрицю Rezult, можна прочитати за допомогою команди load(“Rezult.RData”)