Побудова лінійного графіку для трьох змінних (line chart)

Спочатку створимо матрицю з даними, які будемо візуалізувати (використання функції sample() пояснено далі).

mtr <- matrix(NA, 3, 20)
mtr[ , 1] <- sample(10:30, 3)
for (i in 2:20) {mtr[ , i] <- mtr[ , i-1] + sample(-2:2, 3)} 
mtr

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,]   29   27   28   28   26   24   24   24   23    22    23    23    21    22
## [2,]   10   12   14   13   13   15   17   15   17    19    19    20    19    21
## [3,]   11   10   10    8   10   11    9    8    9     9    11    10    12    11
##      [,15] [,16] [,17] [,18] [,19] [,20]
## [1,]    24    26    25    24    23    21
## [2,]    20    18    19    17    18    19
## [3,]     9     9     9    10     8     8

Тепер можна побудувати графік

plot(mtr[1, ], type="l", lty=2, col="red", # Перша лінія
     main = "Line chart", # Заголовок
     ylim = c(0, max(mtr)), # Масштаб осі ординат
     xlab="First axis", ylab="Second axis") # Підписи осей
lines(mtr[2, ], type="l", lty=3, col="blue") # Друга лінія
lines(mtr[3, ], type="l", lty=4, col="green") # Третя лінія
legend(2, 10, inset=.05, title="Позначення:", c("Перша лініяі","Друга лінія", "Третя лінія"), # Легенда
       lty=c(2, 3, 4), col=c("red", "blue", "green")) # Порівнювані елементи в легенді

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

Можна використовувати й функцію floor(), і функцію trunc(). До величини, яку треба ймовірнісно округлити, додається випадкова величина, що рівномірно розподілена між 0 та 1; сума округлюється до нижчого (у разі позитивних чисел) цілого.

У наведеному прикладі середнє зі 100 округлень величини 3.3 має бути достатньо близьким до цієї величини.

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

## [1] 3.296

Формування випадкової вибірки за допомогою функції sample ()

Функція sample(x, size, replace = FALSE, prob = NULL) може використовуватися як для генерації наборів випадкових цілих чисел, так і для вибору певної кількості об’єктів з їх сукупності. Два перші аргументи є обов’язковими: x — це діапазон, з якого здійснюється вибір, а size — розмір вибірки. Два наступних параметри є опціональними. Якщо вони не вказані, R використовує ті значення, що наведені у наведеному зразку повного запису аргументів. Так, параметр replace визначає, чи можна повторно обирати елемент, який вже був вибраний. У разі replace = FALSE після того, як обрано перший елемент, обране значення забирається з-поміж тих, що можуть бути обраними на наступному кроці. Якщо вказати replace=TRUE, команда обиратиме кожний наступний елемент з повної сукупності.

sample(1:10, 9)

## [1] 2 9 7 8 4 1 3 6 5

# sample(1:10, 11) # Тут мало б бути повідомлення про помилку
sample(1:10, 11, replace=T) # Тут усе нормально: елементи можна обирати повторно

##  [1] 6 2 8 1 9 4 3 1 6 9 7

Моделювання добору за допомогою функції sample () та “вектора долі”

У попередньому прикладі наведено зразок використання функція sample(). Її можна використати для моделювання добору. Припустимо, проти особин, що перелічені у векторі Al_Pop перелічені особини, проти яких діє добір s. Величина s лежить у межах від 0 (усі виживають) до 1 (усі гинуть). Особини, на користь яких працює добір, у цьому випадку виживають усі. Звісно, неконкурентна смертність, що однаково проріджує як тих, так і інших, може бути скільки завгодно потужною.

s <- sample(2:5, 1) / 10 # Випадково обирається значення добору
s

## [1] 0.4

Al_Pop <- sample(1:100, 12) # Сукупність, яка зазнає дії добору
Al_Pop

##  [1] 37 47 56 49 15 27 35 65 74 82 62  4

Fate <- c(rep(NA, s*10), rep(1, (10-s*10))) # Вектор "долі"
Fate # (один на весь скрипт; він не є випадковим)

##  [1] NA NA NA NA  1  1  1  1  1  1

Lot <- sample(Fate, length(Al_Pop), replace = T) # Вектор "жереба" 
Lot # (створюється окремо для кожного скорочення)

##  [1] NA NA NA  1  1 NA  1 NA  1  1  1 NA

Be_Pop <- Al_Pop * Lot # Скорочена сукупність; загиблі - NA
Be_Pop

##  [1] NA NA NA 49 15 NA 35 NA 74 82 62 NA

Ga_Pop <- na.omit(Be_Pop) # Остаточно скорочена сукупність
Ga_Pop

## [1] 49 15 35 74 82 62

Перемикач у разі множинного вибору: switch()

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

Randoms = sample(10:99, 5)
Randoms

## [1] 27 32 63 93 56

Rando <- sort(unique(Randoms))
Rando

## [1] 27 32 56 63 93

x <- sample(Rando, 1)
x

## [1] 56

Way <- which(x == Rando)
switch(Way, 
       print(Rando[1]),
       {print(Rando[2]); print("Second")},
       {print(Rando[3]); (a <- x^2)},
       print(Rando[4]),
       print(Rando[5]))

## [1] 56

## [1] 3136

Перебір значень початкових параметрів

Застосовується у моделях ІІІ типу

vector_First <- c(100, 200, 300, 400)
vector_Second <- c(0.1, 0.2, 0.3, 0.4, 0.5)
Outcome <- matrix(NA, nrow = length(vector_First), ncol = length(vector_Second), dimnames = list(vector_First, vector_Second))
for (a in 1:length(vector_First)){
   First <- vector_First[a] 
   for (b in 1:length(vector_Second)){ 
      Second <- vector_Second[b] 
      Outcome[a, b] <- First + Second * 100 # Тут розташовані основні розрахунки результату
   }}
Outcome

##     0.1 0.2 0.3 0.4 0.5
## 100 110 120 130 140 150
## 200 210 220 230 240 250
## 300 310 320 330 340 350
## 400 410 420 430 440 450

Тривимірна візуалізація результатів перебору параметрів

У попередньому прикладі згенеровано матрицю Outcome. Візуалізувати її можна таким чином

# install.packages("plot3D")
library(plot3D)
persp3D(z = Outcome, x = vector_First, y = vector_Second, theta = 330, phi = 30,
        zlab = "\nРезультат", xlab = "First", ylab = "Second",
        expand = 0.6, facets = FALSE, scale = TRUE,  ticktype = "detailed",nticks = 5, cex.axis=0.9, 
        clab = "Кольорова\nшкала", colkey = list(side = 2, length = 0.5))

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

Можна використовувати і функцію floor(), і функцію trunc(). До величини, яку треба ймовірнісно округлити, додається випадкова величина, що рівномірно розподілена між 0 та 1; сума округлюється до нижчого (у разі позитивних чисел) цілого.

У наведеному прикладі середнє з 100 округлень величини 3.3 має бути достатньо близьким до цієї величини.

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

## [1] 3.309

Формування випадкової вибірки за допомогою функції sample

Функція sample(x, size, replace = FALSE, prob = NULL) може використовуватися як для генерації наборів випадкових цілих чисел, так і для вибору певної кількості об’єктів з їх сукупності. Два перші аргументи є обов’язковими: x — це діапазон, з якого здійснюється вибір, а size — розмір вибірки. Два наступних параметри є опціональними. Якщо вони не вказані, R використовує ті значення, що наведені у наведеному зразку повного запису аргументів. Так, параметр replace визначає, чи можна повторно обирати елемент, який вже був вибраний. У разі replace = FALSE після того, як обрано перший елемент, обране значення забирається з числа тих, які можуть бути обраними. Якщо вказати replace=TRUE, команда обиратиме кожний наступний елемент з повної сукупності.

sample(1:10, 9)

## [1]  6  1  4  7  8  3 10  9  5

# sample(1:10, 11) # Тут мало б бути повідомлення про помилку
sample(1:10, 11, replace=T) # Тут усе нормально: елементи можна обирати повторно

##  [1]  6 10  8  7  1  3  9  1  3 10  5

Моделювання добору за допомогою функції sample та “вектора долі”

У попередньому прикладі наведено зразок використання функція sample(). Її можна використати для моделювання добору. Припустимо, проти особин, що перелічені у векторі Al_Pop перелічені особини, проти яких діє добір s. Величина s лежить у межах від 0 (усі виживають) до 1 (усі гинуть). Особини, на користь яких працює добір, у цьому випадку виживають усі. Звісно, неконкурентна смертність, що однаково проріджує як тих, так і інших, може бути скільки завгодно потужною.

s <- sample(2:5, 1) / 10 # Випадково обирається значення добору
s

## [1] 0.2

Al_Pop <- sample(1:100, 12) # Сукупність, яка зазнає дії добору
Al_Pop

##  [1] 10 49 46 74 89 65 67 80  7 31 14 25

Fate <- c(rep(NA, s*10), rep(1, (10-s*10))) # Вектор "долі"
Fate # (один на весь скрипт; він не є випадковим)

##  [1] NA NA  1  1  1  1  1  1  1  1

Lot <- sample(Fate, length(Al_Pop), replace = T) # Вектор "жереба" 
Lot # (створюється окремо для кожного скорочення)

##  [1] NA  1  1  1  1  1  1  1  1  1  1 NA

Be_Pop <- Al_Pop * Lot # Скорочена сукупність; загиблі - NA
Be_Pop

##  [1] NA 49 46 74 89 65 67 80  7 31 14 NA

Ga_Pop <- na.omit(Be_Pop) # Остаточно скорочена сукупність
Ga_Pop

##  [1] 49 46 74 89 65 67 80  7 31 14
## attr(,"na.action")
## [1]  1 12
## attr(,"class")
## [1] "omit"

Перемикач у разі множинного вибору: switch()

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

Randoms = sample(10:99, 5)
Randoms

## [1] 85 47 43 93 26

Rando <- sort(unique(Randoms))
Rando

## [1] 26 43 47 85 93

x <- sample(Rando, 1)
x

## [1] 93

Way <- which(x == Rando)
switch(Way, 
       print(Rando[1]),
       {print(Rando[2]); print("Second")},
       {print(Rando[3]); (a <- x^2)},
       print(Rando[4]),
       print(Rando[5]))

## [1] 93

Перебір значень початкових параметрів

Застосовується у моделях ІІІ типу

vector_First <- c(100, 200, 300, 400)
vector_Second <- c(0.1, 0.2, 0.3, 0.4, 0.5)
Outcome <- matrix(NA, nrow = length(vector_First), ncol = length(vector_Second), dimnames = list(vector_First, vector_Second))
for (a in 1:length(vector_First)){
   First <- vector_First[a] 
   for (b in 1:length(vector_Second)){ 
      Second <- vector_Second[b] 
      Outcome[a, b] <- First + Second * 100 # Тут розташовані основні розрахунки результату
   }}
Outcome

##     0.1 0.2 0.3 0.4 0.5
## 100 110 120 130 140 150
## 200 210 220 230 240 250
## 300 310 320 330 340 350
## 400 410 420 430 440 450

Тривимірна візуалізація результатів перебору параметрів

У попередньому прикладі згенеровано матрицю Outcome. Візуалізувати її можна таким чином

# install.packages("plot3D")
library(plot3D)
persp3D(z = Outcome, x = vector_First, y = vector_Second, theta = 330, phi = 30, # Побудова графіку
        zlab = "\nРезультат", xlab = "First", ylab = "Second",
        expand = 0.6, facets = FALSE, scale = TRUE,  ticktype = "detailed",nticks = 5, cex.axis=0.9, 
        clab = "Кольорова\nшкала", colkey = list(side = 2, length = 0.5))

Власна функція конкурентного скорочення чисельності (імітація на рівні груп)

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

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

Виходом є нові, скорочені чисельності груп.

Наведемо два варіанти. Перший (CompetitiveReduction0) носить детермінований характер, та таз за разом буде призводити до того ж самого результату.

CompetitiveReduction0 <- function(Capacity, Amounts, Competitiveness){
  RelativCompetit <- Competitiveness/max(Competitiveness); Quoted <- RelativCompetit * Amounts
  Way <- ifelse(sum(Amounts) <= Capacity, 1, ifelse(sum(Quoted) >= Capacity, 2, 3))
  Reduced <- switch(Way, # Defines a formula for calculations based on value Way
                    Amounts, # If Way==1
                    round(Quoted*Capacity/sum(Quoted)), # If Way==2, & on next row — if Way==3:
                    round(Amounts-(Amounts-Quoted)*(sum(Amounts)-Capacity)/(sum(Amounts)-sum(Quoted))) )}

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

CompetitiveReduction <- function(Capacity, Amounts, Competitiveness){
  RelativCompetit <- Competitiveness/max(Competitiveness); Quoted <- RelativCompetit * Amounts
  Way <- ifelse(sum(Amounts) <= Capacity, 1, ifelse(sum(Quoted) >= Capacity, 2, 3))
  Reduced <- switch(Way, # Defines a formula for calculations based on value Way
                    Amounts, # If Way==1
                    floor(Quoted*Capacity/sum(Quoted) + runif(1)), # If Way==2, & on next row — if Way==3:
                    floor(Amounts-(Amounts-Quoted)*(sum(Amounts)-Capacity)/(sum(Amounts)-sum(Quoted)) + runif(1)) )}

Таке скорочення чисельності задовольняє двом важливим умовам: 1) сумарна чисельність груп до скорочення знижується до величини, яка відповідає забезпеченості ресурсами; 2) чисельність кожної групи до конкурентного скорочення Amounts[i] знижується до такої чисельності після конкурентного скорочення Reduced[i], що в кожній групі частка особин, які збереглися після конкурентного скорочення, пропорційна конкурентоспроможності представників цієї групи: Competitiveness[i]/max(Competitiveness).

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

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

set.seed(123456)
K <- 20
rws <- c("A", "B", "C", "D", "E", "F")
cls <- c("Codes", "Compet", "s", "Al", "Be")
ModelPop <- matrix(NA, nrow = length(rws), ncol = length(cls), dimnames = list(rws, cls))
ModelPop[, "Codes"] <- 1:6
ModelPop[, "Compet"] <- c(1, 0.9, 0.1, 0.8, 0.5, 0.4)
ModelPop[, "Al"] <- c(12, 20, 15, 6, 5, 2)
ModelPop[, "Be"] <- CompetitiveReduction(K, ModelPop[, "Al"], ModelPop[, "Compet"])
ModelPop

##   Codes Compet  s Al Be
## A     1    1.0 NA 12  6
## B     2    0.9 NA 20  9
## C     3    0.1 NA 15  1
## D     4    0.8 NA  6  3
## E     5    0.5 NA  5  2
## F     6    0.4 NA  2  1

sum(ModelPop[, "Be"])

## [1] 22

Від альфа-складу популяції ми перейшли до бета-складу. Чому бета-чисельність не чітко відповідає обмеженню ємності середовища (22 особини замість 20)? Це наслідок випадковості при реалізації алгоритму скорочення чисельності.

Скорочення чисельності популяції внаслідок добору (імітація на рівні особин)

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

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

ModelPop[, "s"] <- 1- ModelPop[, "Compet"] 

Створимо вектор Al_Pop, що складається з кодів кожної особини окремо.

Al_Pop <- c(rep(ModelPop[1, "Codes"], ModelPop[1, "Al"]), 
            rep(ModelPop[2, "Codes"], ModelPop[2, "Al"]), 
            rep(ModelPop[3, "Codes"], ModelPop[3, "Al"]), 
            rep(ModelPop[4, "Codes"], ModelPop[4, "Al"]), 
            rep(ModelPop[5, "Codes"], ModelPop[5, "Al"]), 
            rep(ModelPop[6, "Codes"], ModelPop[6, "Al"]))
Al_Pop

##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3
## [39] 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 5 5 5 5 5 6 6

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

Selection <- rep(NA, length(Al_Pop))
Randoms <- runif(length(Al_Pop), 0, 1)
for (n in 1:length(Al_Pop)) {Selection[n] <- ModelPop[which(ModelPop[, "Codes"]==Al_Pop[n]), "s"]}
Be_Pop <- Al_Pop * (floor(Randoms - Selection) + 1)
Be_Pop <- Be_Pop[-which(Be_Pop==0)] 
Be_Pop

##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 4 4 4 4 4 4
## [39] 5 5 5 5

table(Be_Pop)

## Be_Pop
##  1  2  3  4  5 
## 12 19  1  6  4

length(Be_Pop)

## [1] 42

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

Неконкурентне імітаційне скорочення чисельності (імітація на рівні особин)

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

if(length(Be_Pop) > K) {
  Randoms <- runif(length(Be_Pop), 0, 1)  
  Ga_Pop <- Be_Pop * (floor(Randoms - (length(Be_Pop)-K)/length(Be_Pop)) + 1)
  Ga_Pop <- Ga_Pop[-which(Ga_Pop==0)]   }
Ga_Pop

##  [1] 1 1 1 1 2 2 2 2 2 2 2 2 4 4 5

table(Ga_Pop)

## Ga_Pop
## 1 2 4 5 
## 4 8 2 1

length(Ga_Pop)

## [1] 15

Лишилося менше особин, ніж мало залишитися. Випадковість…

Cons <- 10
z <- 0
for (f in 1:Cons) {
  print(paste(f, "from", Cons)) 
  y <- 0
  for (s in 1:Cons) {
    x <- 0
    for (t in 1:Cons) {
      w <- 0
      for (f in 1:Cons) {
        w <- w + f}
      x <- x + w}
    y <- y + x}
  z <- z + y}

## [1] "1 from 10"
## [1] "2 from 10"
## [1] "3 from 10"
## [1] "4 from 10"
## [1] "5 from 10"
## [1] "6 from 10"
## [1] "7 from 10"
## [1] "8 from 10"
## [1] "9 from 10"
## [1] "10 from 10"

z

## [1] 55000