Скрипты для моделирования структуры сообществ

ДЛЯ ПЕРВОГО ЗАДАНИЯ


# задаю исходные параметры

z <- 0.5; S <- 100


# строю первую модель Мотомуры-Мэй по формуле Мотомуры
max_ob <- 900

# объявляю пустой вектор
n_1 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_1[i] <- max_ob * (z ^ (i - 1))
}

# строю вторую модель Мотомуры-Мэй по формуле Мотомуры
max_ob <- 900

# объявляю пустой вектор
n_2 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_2[i] <- max_ob * (z ^ (i - 1))
}

# строю третью модель Мотомуры-Мэй по формуле Мотомуры
max_ob <- 900

# объявляю пустой вектор
n_3 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_3[i] <- max_ob * (z ^ (i - 1))
}

# строю четвёртую модель Мотомуры-Мэй по формуле Мотомуры
max_ob <- 900

# объявляю пустой вектор
n_4 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_4[i] <- max_ob * (z ^ (i - 1))
}

# строю пятую модель Мотомуры-Мэй по формуле Мотомуры
max_ob <- 900

# объявляю пустой вектор
n_5 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_5[i] <- max_ob * (z ^ (i - 1))
}

# строю шестую модель Мотомуры-Мэй по формуле Мотомуры
max_ob <- 900

# объявляю пустой вектор
n_6 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_6[i] <- max_ob * (z ^ (i - 1))
}

# строю седьмую модель Мотомуры-Мэй по формуле Мотомуры
max_ob <- 900

# объявляю пустой вектор
n_7 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_7[i] <- max_ob * (z ^ (i - 1))
}

# строю график модели
plot(n_1, type = 'l', xlab = 'Ранг вида', ylab = 'Обилие вида', col = 'red')
lines(n_2, col = 'blue')
lines(n_3, col = 'green')
lines(n_4, col = 'yellow')
lines(n_5, col = 'orange')
lines(n_6, col = 'coral')
lines(n_7, col = 'pink')
title('Модель Мотомуры-Мэй, формула Мотомуры при разных значениях максимального обилия')


ДЛЯ ВТОРОГО ЗАДАНИЯ


# Построение первой модели Мотомуры-Мэй по формуле Мэй

# задаю исходные параметры
k <- 0.9
N <- 1000
S <- 100

# объявляю пустой вектор
n_1 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_1[i] <- (N * k * ((1 - k) ^ (i - 1))) / (1 - (1 - k) ^ S)
}

# объявляю пустой вектор
n_0 <- NULL

# рассчитываю индекс Шеннона
for(i in 1:S) {
  n_0[i] <- (n_1[i] / N) * log(n_1[i] / N)
}
H_1 <- -sum(n_0)

# рассчитываю выравненость индекса Шеннона
E_1 <- H_1 / log(S)

# Построение второй модели Мотомуры-Мэй по формуле Мэй

# задаю исходные параметры
k <- 0.9
N <- 1000
S <- 100

# объявляю пустой вектор
n_2 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_2[i] <- (N * k * ((1 - k) ^ (i - 1))) / (1 - (1 - k) ^ S)
}

# объявляю пустой вектор
n_0 <- NULL

# рассчитываю индекс Шеннона
for(i in 1:S) {
  n_0[i] <- (n_2[i] / N) * log(n_2[i] / N)
}
H_2 <- -sum(n_0)

# рассчитываю выравненость индекса Шеннона
E_2 <- H_2 / log(S)

# Построение третьей модели Мотомуры-Мэй по формуле Мэй

# задаю исходные параметры
k <- 0.9
N <- 1000
S <- 100

# объявляю пустой вектор
n_3 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_3[i] <- (N * k * ((1 - k) ^ (i - 1))) / (1 - (1 - k) ^ S)
}

# объявляю пустой вектор
n_0 <- NULL

# рассчитываю индекс Шеннона
for(i in 1:S) {
  n_0[i] <- (n_3[i] / N) * log(n_3[i] / N)
}
H_3 <- -sum(n_0)

# рассчитываю выравненость индекса Шеннона
E_3 <- H_3 / log(S)

# Построение четвёртой модели Мотомуры-Мэй по формуле Мэй

# задаю исходные параметры
k <- 0.9
N <- 1000
S <- 100

# объявляю пустой вектор
n_4 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_4[i] <- (N * k * ((1 - k) ^ (i - 1))) / (1 - (1 - k) ^ S)
}

# объявляю пустой вектор
n_0 <- NULL

# рассчитываю индекс Шеннона
for(i in 1:S) {
  n_0[i] <- (n_4[i] / N) * log(n_4[i] / N)
}
H_4 <- -sum(n_0)

# рассчитываю выравненость индекса Шеннона
E_4 <- H_4 / log(S)

# Построение пятой модели Мотомуры-Мэй по формуле Мэй

# задаю исходные параметры
k <- 0.9
N <- 1000
S <- 100

# объявляю пустой вектор
n_5 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_5[i] <- (N * k * ((1 - k) ^ (i - 1))) / (1 - (1 - k) ^ S)
}

# объявляю пустой вектор
n_0 <- NULL

# рассчитываю индекс Шеннона
for(i in 1:S) {
  n_0[i] <- (n_5[i] / N) * log(n_5[i] / N)
}
H_5 <- -sum(n_0)

# рассчитываю выравненость индекса Шеннона
E_5 <- H_5 / log(S)

# Построение шестой модели Мотомуры-Мэй по формуле Мэй

# задаю исходные параметры
k <- 0.9
N <- 1000
S <- 100

# объявляю пустой вектор
n_6 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_6[i] <- (N * k * ((1 - k) ^ (i - 1))) / (1 - (1 - k) ^ S)
}

# объявляю пустой вектор
n_0 <- NULL

# рассчитываю индекс Шеннона
for(i in 1:S) {
  n_0[i] <- (n_6[i] / N) * log(n_6[i] / N)
}
H_6 <- -sum(n_0)

# рассчитываю выравненость индекса Шеннона
E_6 <- H_6 / log(S)

# Построение седьмой модели Мотомуры-Мэй по формуле Мэй

# задаю исходные параметры
k <- 0.9
N <- 1000
S <- 100

# объявляю пустой вектор
n_7 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_7[i] <- (N * k * ((1 - k) ^ (i - 1))) / (1 - (1 - k) ^ S)
}

# объявляю пустой вектор
n_0 <- NULL

# рассчитываю индекс Шеннона
for(i in 1:S) {
  n_0[i] <- (n_7[i] / N) * log(n_7[i] / N)
}
H_7 <- -sum(n_0)

# рассчитываю выравненость индекса Шеннона
E_7 <- H_7 / log(S)

# собираю все значения индекса Шеннона и его выравненности в одну таблицу
H <- c(H_1, H_2, H_3, H_4, H_5, H_6, H_7)
tab_it <- as.data.frame(H)
tab_it$E <- c(E_1, E_2, E_3, E_4, E_5, E_6, E_7)
tab_it$id <- 1:nrow(tab_it)

# строю график распределения точек в пространстве H-E
plot(tab_it$H, tab_it$E, ylim = c(0, 1), col = "red", xlab = "Индекс Шеннона", ylab = "Выравненность индекса Шеннона")
text(tab_it$H, tab_it$E, labels = tab_it$id, cex = 0.75, pos = 3)
title('График распределения сообществ в пространстве индекса Шеннона и его выравненности')

ДЛЯ ТРЕТЬЕГО ЗАДАНИЯ


# задаю исходные параметры
z <- 0.5
N <- 1000
S <- 100
max_ob <- 100

# Построение модели Мотомуры-Мэй по формуле Мотомуры

# объявляю пустой вектор
n_1 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_1[i] <- max_ob * (z ^ (i - 1))
}

# объявляю пустой вектор
n_0 <- NULL

# рассчитываю индекс Шеннона
for(i in 1:S) {
  n_0[i] <- (n_1[i] / N) * log(n_1[i] / N)
}
H_1 <- -sum(n_0)

# рассчитываю выравненость индекса Шеннона
E_1 <- H_1 / log(S)

# Построение модели Мотомуры-Мэй по формуле Мэй

# задаю исходные параметры
k <- z

# объявляю пустой вектор
n_2 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_2[i] <- (N * k * ((1 - k) ^ (i - 1))) / (1 - (1 - k) ^ S)
}

# объявляю пустой вектор
n_0 <- NULL

# рассчитываю индекс Шеннона
for(i in 1:S) {
  n_0[i] <- (n_2[i] / N) * log(n_2[i] / N)
}
H_2 <- -sum(n_0)

# рассчитываю выравненость индекса Шеннона
E_2 <- H_2 / log(S)


# Построение гиперболической модели Левича

# задаю исходные параметры
b <- z
K <- max_ob

# объявляю пустой вектор
n_3 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_3[i] <- K * (i ^ -b)
}

# объявляю пустой вектор
n_0 <- NULL

# рассчитываю индекс Шеннона
for(i in 1:S) {
  n_0[i] <- (n_3[i] / N) * log(n_3[i] / N)
}
H_3 <- -sum(n_0)

# рассчитываю выравненость индекса Шеннона
E_3 <- H_3 / log(S)

# Построение дзета-модели Левича

# задаю исходные параметры
x <- z
y <- z
N_1 <- max_ob

# объявляю пустой вектор
n_4 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_4[i] <- N_1 * ((x ^ (i-1)) / (i ^ y))
}

# объявляю пустой вектор
n_0 <- NULL

# рассчитываю индекс Шеннона
for(i in 1:S) {
  n_0[i] <- (n_4[i] / N) * log(n_4[i] / N)
}
H_4 <- -sum(n_0)

# рассчитываю выравненость индекса Шеннона
E_4 <- H_4 / log(S)


# Построение модели "разломанного стержня" Макартура(1)

# объявляю пустой вектор
n_5 <- NULL
n <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  for(l in i:S){
    n[l] <- (l ^ (-1)) / S
  }
  n <- n[!is.na(n)]
  n_5[i] <- (S ^ 0.5) * sum(n)
  n <- NULL
}

# объявляю пустой вектор
n_0 <- NULL

# рассчитываю индекс Шеннона
for(i in 1:S) {
  n_0[i] <- (n_5[i] / N) * log(n_5[i] / N)
}
H_5 <- -sum(n_0)

# рассчитываю выравненость индекса Шеннона
E_5 <- H_5 / log(S)


# Построение модели "разломанного стержня" Макартура(2)

# объявляю пустой вектор
n_6 <- NULL
n <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  for(j in i:S){
    n[j] <- 1 / j
  }
  n <- n[!is.na(n)]
  n_6[i] <- (N / S) * sum(n)
  n <- NULL
}

# объявляю пустой вектор
n_0 <- NULL

# рассчитываю индекс Шеннона
for(i in 1:S) {
  n_0[i] <- (n_6[i] / N) * log(n_6[i] / N)
}
H_6 <- -sum(n_0)

# рассчитываю выравненость индекса Шеннона
E_6 <- H_6 / log(S)


# Построение модели "разломанного стержня" Макартура(3)

# объявляю пустой вектор
n_7 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_7[i] <- (((i + 1) ^ 0.5) - (i ^ 0.5)) / (S ^ 0.5)
}

# объявляю пустой вектор
n_0 <- NULL

# рассчитываю индекс Шеннона
for(i in 1:S) {
  n_0[i] <- (n_7[i] / N) * log(n_7[i] / N)
}
H_7 <- -sum(n_0)

# рассчитываю выравненость индекса Шеннона
E_7 <- H_7 / log(S)

# строю график модели
plot(n_1, type = 'l', xlab = 'Ранг вида', ylab = 'Обилие вида', col = 'red')
lines(n_2, col = 'blue')
lines(n_3, col = 'green')
lines(n_4, col = 'yellow')
lines(n_5, col = 'orange')
lines(n_6, col = 'coral')
lines(n_7, col = 'black')
title('Модели сообщества вида "ранг - обилие"')

# собираю все значения индекса Шеннона и его выравненности в одну таблицу
H <- c(H_1, H_2, H_3, H_4, H_5, H_6, H_7)
tab_it <- as.data.frame(H)
tab_it$E <- c(E_1, E_2, E_3, E_4, E_5, E_6, E_7)
tab_it$id <- 1:nrow(tab_it)

# строю график распределения точек в пространстве H-E
plot(tab_it$H, tab_it$E, ylim = c(0, 1), col = "red", xlab = "Индекс Шеннона", ylab = "Выравненность индекса Шеннона")
text(tab_it$H, tab_it$E, labels = tab_it$id, cex = 0.75, pos = 3)

ДЛЯ ЧЕТВЁРТОГО ЗАДАНИЯ


# ввожу исходные данные и сортирую по убыванию
soob <- c(0.5, 15, 0.5, 2, 12, 1, 7, 4, 1, 20, 0.5, 1, 0.5, 7, 0.5)
soob <- sort(soob, decreasing = T)

# задаю исходные параметры
z <- 0.5
N <- sum(soob)
S <- length(soob)
max_ob <- max(soob)

# Построение модели Мотомуры-Мэй по формуле Мотомуры

# объявляю пустой вектор
n_1 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_1[i] <- max_ob * (z ^ (i - 1))
}

# Построение модели Мотомуры-Мэй по формуле Мэй

# задаю исходные параметры
k <- z

# объявляю пустой вектор
n_2 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_2[i] <- (N * k * ((1 - k) ^ (i - 1))) / (1 - (1 - k) ^ S)
}

View(itog)
# Построение гиперболической модели Левича

# задаю исходные параметры
b <- z
K <- max_ob

# объявляю пустой вектор
n_3 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_3[i] <- K * (i ^ -b)
}

# Построение дзета-модели Левича

# задаю исходные параметры
x <- z
y <- z
N_1 <- max_ob

# объявляю пустой вектор
n_4 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_4[i] <- N_1 * ((x ^ (i-1)) / (i ^ y))
}

# Построение модели "разломанного стержня" Макартура(1)

# объявляю пустой вектор
n_5 <- NULL
n <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  for(l in i:S){
    n[l] <- (l ^ (-1)) / S
  }
  n <- n[!is.na(n)]
  n_5[i] <- (S ^ 0.5) * sum(n)
  n <- NULL
}

# Построение модели "разломанного стержня" Макартура(2)

# объявляю пустой вектор
n_6 <- NULL
n <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  for(j in i:S){
    n[j] <- 1 / j
  }
  n <- n[!is.na(n)]
  n_6[i] <- (N / S) * sum(n)
  n <- NULL
}

# Построение модели "разломанного стержня" Макартура(3)

# объявляю пустой вектор
n_7 <- NULL

# рассчитываю саму модель
for(i in 1:S) {
  n_7[i] <- (((i + 1) ^ 0.5) - (i ^ 0.5)) / (S ^ 0.5)
}

# рассчитываю корреляции Спирмена
S_1 <- cor(n_1, soob, method = "spearman")
S_2 <- cor(n_2, soob, method = "spearman")
S_3 <- cor(n_3, soob, method = "spearman")
S_4 <- cor(n_4, soob, method = "spearman")
S_5 <- cor(n_5, soob, method = "spearman")
S_6 <- cor(n_6, soob, method = "spearman")
S_7 <- cor(n_7, soob, method = "spearman")

# рассчитываю корреляции Пирсона
P_1 <- cor(log(n_1), log(soob))
P_2 <- cor(log(n_2), log(soob))
P_3 <- cor(log(n_3), log(soob))
P_4 <- cor(log(n_4), log(soob))
P_5 <- cor(log(n_5), log(soob))
P_6 <- cor(log(n_6), log(soob))
P_7 <- cor(log(n_7), log(soob))

# собираю корреляции в единую таблицу и вывожу её
CS <- c(S_1, S_2, S_3, S_4, S_5, S_6, S_7)
itog <- as.data.frame(CS)
CP <- c(P_1, P_2, P_3, P_4, P_5, P_6, P_7)
itog$CP <- CP
itog$id <- 1:nrow(itog)
itog

# строю график модели
plot(x = log(soob), y = log(1:length(soob)), type = 'l', xlab = 'Логарифм ранга вида',
     ylab = 'Логарифм обилия вида', col = 'black')
lines(log(n_1), y = log(1:length(soob)), col = 'red', lty = 1)
lines(log(n_2), y = log(1:length(soob)), col = 'red', lty = 2)
lines(log(n_3), y = log(1:length(soob)), col = 'green', lty = 1)
lines(log(n_4), y = log(1:length(soob)), col = 'green', lty = 2)
lines(log(n_5), y = log(1:length(soob)), col = 'blue', lty = 1)
lines(log(n_6), y = log(1:length(soob)), col = 'blue', lty = 2)
lines(log(n_7), y = log(1:length(soob)), col = 'blue', lty = 5)
title('Логарифмированные модели сообщества вида "ранг - обилие"
      (чёрным цветом показан график реального сообщества)')


ДЛЯ ЧАСТОТНОЙ (ЛОГАРИФМИЧЕСКОЙ) МОДЕЛИ


# ввожу исходные данные по обилию видов на площадке
pp <- c(0.5, 15, 0.5, 2, 12, 1, 7, 4, 1, 20, 0.5, 1, 0.5, 7, 0.5)

# определяю уникальные обилия и сортирую их по возрастанию
un_p <- unique(pp)
un_p <- sort(un_p)

# задаю пустой вектор для частот и равную нулю переменную для счёта
ch_p <- NULL
l <- 0

# подсчитываю частоты для каждого значения обилия
for (i in 1:length(un_p)) {
  for (j in 1:length(pp)) {
    if(un_p[i] == pp[j]) {
      l <- l + 1
    }
    else {
      l <- l
    }
  }
  ch_p[i] <- l
  l <- 0
}

# рассчитываю общее число видов и общее число особей
S <- length(pp)
N <- sum(pp)

# задаю исходное значение для х
x <- 0.001

# подбираю подходящее значение х
while(S / N < (1 - x) / x * (-log(1-x))) {
  x <- x + 0.001
}

# определяю а
a <- (N * (1 - x)) / x

# определяю n
n <- max(pp)

# задаю пустой вектор для теоретических значений обилий
pp_t <- NULL

# рассчитываю теоретические значения обилий
for(i in 1:n) {
  pp_t[i] <- (a * x^i) / i
}

# задаю вектор с границами октав
gr <- c(0, 2, 4, 8, 16, 32, 64)

# задаю пустой вектор для теоретических октав и равный нулю счётчик
ok_t1 <- NULL
l <- 0

# суммирую теоретические значения по октавам
for(i in 2:length(gr)) {
  for(j in 1:length(pp_t)) {
    if(gr[i] >= j & j > gr[i-1]) {
      l = l + pp_t[j]
    }
    else {
      l <- l
    }
  }
  ok_t1[i] <- l
  l <- 0
}

# удаляю NA значение
ok_t1 <- ok_t1[!is.na(ok_t1)]

# задаю пустой вектор для эмпирических октав и равный нулю счётчик
ok_r <- NULL
l <- 0

# суммирую эмпирические значения по октавам
for(i in 2:length(gr)) {
  for(j in 1:length(un_p)) {
    if(gr[i] >= un_p[j] & un_p[j] > gr[i-1]) {
      l = l + ch_p[j]
    }
    else {
      l <- l
    }
  }
  ok_r[i] <- l
  l <- 0
}

# удаляю NA значение
ok_r <- ok_r[!is.na(ok_r)]

# задаю нумерацию шкалы для оси х
id <- c(1:length(ok_t1))

# рисую графики
plot(id, ok_r, type = "l", xlab = "Номер октавы", ylab = "Число видов в октаве", lty = 1)
lines(id, ok_t1, lty = 2)
legend("topright", legend = c("Эмпирический график", "Теоретический график"), lty = c(1, 2))
title("Сравнение модели логарифмического распределения с теоретическим графиком")

# рассчитываю корреляции Спирмена и Кендалла
cor(ok_t1, ok_r, method = "spearman")
cor(ok_t1, ok_r, method = "kendall")

Последнее изменение: Понедельник, 18 октября 2021, 06:45