Форум дисциплины

Частотные модели

 
Изображение пользователя Юрий Бобров
Частотные модели
от Юрий Бобров - Пятница, 24 сентября 2021, 14:10
 
# ввожу исходные данные по обилию видов на площадке
pp <- c(20,2,10,14,15,26,50,22,14,7,38,35,26,22,11,15,2)

# определяю уникальные обилия и сортирую их по возрастанию
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("bottomright", legend = c("Эмпирический график", "Теоретический график"), lty = c(1, 2))
title("Сравнение модели логарифмического распределения с теоретическим графиком")

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

# считаю Хи-квадрат
m1 <- as.table(cbind(ok_t1, ok_r))
chisq.test(m1)

# считаю критерий Колмогорова – Смирнова
ks.test(ok_t1, ok_r)

# логарифмирую числа видов в каждой октаве и убираю пустое значение
# (обратите внимание - у вас может его не быть, а может быть больше)
n_m <- log2(ok_r)
n_m <- n_m[-2]

# нахожу стандартное отклонение
s <- sqrt(sum((n_m - sum(n_m)/length(ok_r))^2) / length(ok_r))

# определяю величину а
a <- sqrt(2 * s^2)

# нахожу номер модальной октавы
mo_o <- match(max(ok_r), ok_r)    

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

# рассчитываю теоретические значения
for (i in 1:length(ok_r)) {
  ifelse(i == mo_o, ok_t2[i] <- ok_r[i], l <- 1)
    ifelse(i < mo_o, ok_t2[i] <- ok_r[mo_o] * exp(-(a^2) * (mo_o - i)^2),
           ok_t2[i] <- ok_r[mo_o] * exp(-(a^2) * (i - mo_o)^2))
}

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

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

# считаю Хи-квадрат
m2 <- as.table(cbind(ok_t2, ok_r))
chisq.test(m2)

# считаю критерий Колмогорова – Смирнова
ks.test(ok_t2, ok_r)