• 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 — комбінації змінюваних параметрів (в моделях ІІІ типу).

• 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.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.1 Побудова концептуальної моделі та система позначень етапів у циклі

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

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

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

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

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

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

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

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

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

• Бета-популяція, ^βP, сукупність особин після загибелі особин максимального віку, ^β_{t+max_a}P = ^α_{t+max_a}P – ^γ_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.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 у коректному запису).

При повторних запусках того самого скрипту, або подібних один до одного скриптів з’являється певна небезпека. У робочому просторі програми, в 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”)