Модель Шеллінга з різними значеннями порогів

У першій та другій частині йшлося про модель Шеллінга, де всі агенти мали однакове значення порогу:
https://rpubs.com/ruslana/709723 (Частина 1)
https://rpubs.com/ruslana/711645 (Частина 2)
Зараз поговоримо про модель Шеллінга, де агенти мають різні значення порогу. Значення порогу для кожного агента генерується випадковим чином в межах заданого інтервалу. Межі інтервалу встановлюються за допомогою параметрів rmin (мінімальне значення порогу) та rmax (максимальне значення порогу). Очевидно, що rmin не може бути більшим за rmax. Якщо rmin = rmax, всі агенти матимуть однаковий поріг. За замовчуванням встановлено мінімальне значення порогу 0.

# Маємо двовимірну решітку 32 Х 32
# rmin - мінімальне значення порогу, 0<=rmin<=1
# rmax - максимальне значення порогу, 0<=rmax<=1
# chs - початкове розташування агентів
# *якщо chs = F, агенти розміщені випадковим чином
# *якщо chs = T, агенти розміщені у шаховому порядку
# g1 - кількість агентів у 1-ій групі, ціле число > 0 
# *якщо chs = T, параметр g1 не має сенсу 
# n - кількість порожніх клітинок, ціле число > 0
# *n + g1 має бути меншим за 1024
# vz - візуалізація
# *якщо vz = "none", візуалізація не виводиться
# *якщо vz = "static", виводимо лише початкове і кінцеве розміщення
# *якщо vz = "dynamic", спостерігаємо рух агентів на кожному кроці
# res - потрібно виводити результати симуляцій чи ні
# *якщо res = F, результати не виводяться
# *якщо res = T, результати виводяться
Schelling_model_unq <- function(rmin = 0, rmax, chs = F, g1 = 462, n = 100, 
                                vz = c("none", "static", "dynamic"), res = F) {
  # Перший крок (початкове розміщення)
  count <- 1
  if (g1 <= 0 | n <= 0 | g1 + n >= 1024 | rmin < 0 | rmin > 1 | rmax < 0 | rmax > 1 | rmin > rmax | 
      is.numeric(rmin) == F | is.numeric(rmax) == F | is.na(rmin) == T | is.na(rmax) == T | is.logical(chs) == F | 
      vz != "none" & vz != "static" & vz != "dynamic" | is.logical(res) == F | is.na(res) == T) {
  } else {
    # Початкове розміщення агентів
    if (chs == F) {
      group <- c(rep(0, n), rep(1, g1), rep(2, 1024 - g1 - n))
      grid <- matrix(sample(group, 1024), 32, 32)
    } else {
      grid <- matrix(rep(c(rep(c(1, 2), 16), rep(c(2, 1), 16)), 16), 32, 32)
      k <- sample(c(floor(n/2), ceiling(n/2)), 1)
      grid[sample(which(grid == 1), k)] <- 0
      grid[sample(which(grid == 2), n - k)] <- 0
    }
  }
  # Візуалізація початкового розміщення агентів
  # Gray - вільні клітинки
  # Red - представники 1-ої групи
  # Dark green - агенти,  представники 2-ої групи
  if (vz == "static") {
    par(mfrow = c(1, 2))
    image(t(apply(grid, 2, rev)), col = c("gray", "red", "dark green"), axes = F, asp = 1)
  }
  if (vz == "dynamic") {
    par(mfrow = c(1, 1))
    image(t(apply(grid, 2, rev)), col = c("gray", "red", "dark green"), axes = F, asp = 1)
  }
  # Крім візуалізації ще будуть результати симуляцій
  df <- NULL
  # Значення порогів
  r <- round(sapply(1:1024, function(x) ifelse(grid[x] == 0, NA, runif(1, min = rmin, max = rmax))), 2)
  names(r) <- grid
  repeat {
    # Встановлюємо сусідів для кожної з клітинок
    w <- cbind(rep(c(1:dim(grid)[1]), dim(grid)[2]), sort(rep(c(1:dim(grid)[1]), dim(grid)[2])))
    d <- as.matrix(dist(w, method = "maximum"))
    a <- apply(d, 1, function(i) grid[i == 1])
    names(a) <- grid
    # Розраховуємо загальну кількість сусідів для кожного агента
    t_ns <- sapply(1:length(grid), function(x) sum(table(a[[x]], exclude = 0)))
    names(t_ns) <- grid
    t_ns <- ifelse(names(t_ns) == "0", NA, t_ns)
    # Розраховуємо кількість сусідів агента, що належать до тієї ж групи
    g_s <- as.data.frame(sapply(1:length(grid), function(x) table(factor(a[[x]], levels = c("0", "1", "2")))))
    colnames(g_s) <- grid
    g_ns <- sapply(1:1024, function(x) g_s[match(colnames(g_s)[x], rownames(g_s)), x])
    names(g_ns) <- grid
    g_ns <- ifelse(names(g_ns) == "0", NA, g_ns)
    # Розраховуємо частку сусідів "своєї" групи
    p <- ifelse(t_ns == 0, 0, g_ns/t_ns)
    names(p) <- grid
    # Визначимо середню частку сусідів "своєї" групи
    p_sg <- round(mean(p, na.rm = T), 3)
    # Визначаємо агентів, задоволених і незадоволених своїм розташуванням
    satisfied <- p
    satisfied[satisfied < r] <- 0
    satisfied[satisfied >= r] <- 1
    # Розрахуємо частку агентів, задоволених своїм розташуванням
    p_st <- round(mean(satisfied, na.rm = T), 3)
    # Формування результатів симуляцій
    df <- rbind(df, c(p_sg, p_st))
    # Незадоволені агенти переходять на вільну клітинку
    if (length(which(is.na(satisfied) == T)) <= length(which(satisfied == 0 & is.na(satisfied) == F))) {
      ls1 <- sample(which(satisfied == 0 & is.na(satisfied) == F), length(which(is.na(satisfied) == T)))
      names(satisfied)[which(is.na(satisfied) == T)] <- names(satisfied)[ls1]
      r[is.na(r) == T] <- r[ls1]
      r[ls1] <- NA
      names(satisfied)[ls1] <- "0"
      names(r) <- names(satisfied)
    } else {
      ls2 <- sample(which(is.na(satisfied) == T), length(which(satisfied == 0 & is.na(satisfied) == F)))
      r[ls2] <- r[which(satisfied == 0 & is.na(satisfied) == F)]
      names(satisfied)[ls2] <- names(satisfied)[which(satisfied == 0 & is.na(satisfied) == F)] 
      r[which(satisfied == 0 & is.na(satisfied) == F)] <- NA
      names(satisfied)[which(satisfied == 0 & is.na(satisfied) == F)] <- "0"
      names(r) <- names(satisfied)
    }
    # Наступний крок
    count <- count + 1
    # Нове розміщення
    grid <- matrix(as.numeric(names(satisfied)), 32, 32)
    # Візуалізація розміщення агентів на кожному кроці
    if (vz == "dynamic") {
      image(t(apply(grid, 2, rev)), col = c("gray", "red", "dark green"), axes = F, asp = 1)
    }
    # Умови для завершення руху агентів
    if (mean(satisfied, na.rm = T) == 1 | count == 300) break
  }
  # Візуалізація кінцевого розміщення агентів
  if (vz == "static") {
    image(t(apply(grid, 2, rev)), col = c("gray", "red", "dark green"), axes = F, asp = 1)
  }
  # Результати симуляцій
  if (res == T) {
    df <- as.data.frame(df)
    colnames(df) <- c("Середня частка своїх", "Частка задоволених")
    df
  }
}

Переходимо до симуляцій і розпочнемо з випадкового розміщення. Розглянемо такі варіанти максимального значення порогу: 0.3, 0.5, 0.7.

Schelling_model_unq(rmax = 0.3, vz = "static")

Schelling_model_unq(rmax = 0.5, vz = "static")

Schelling_model_unq(rmax = 0.7, vz = "static")

Для шахового розміщення розглянемо ті ж значення rmax: 0.3, 0.5, 0.7.

Schelling_model_unq(rmax = 0.3, chs = T, vz = "static")

Schelling_model_unq(rmax = 0.5, chs = T, vz = "static")

Schelling_model_unq(rmax = 0.7, chs = T, vz = "static")

Для випадкового розміщення тримали такі результати. За значення порогу 0.3 візуально різниці між початковим та кінцевим розміщенням агентів майже не видно. А для значень порогів 0.5 і 0.7 вже можна спостерігати сегрегацію. Для шахового розміщення за значення порогу 0.3 відмінностей між початковим і кінцевим розміщенням теж практично не видно. За значення порогу 0.5 сегрегація наявна, однак візуально є досить невираженою. За значення 0.7 чітко видно утворення "кварталів", де проживають представники переважно однієї групи. В кого буде натхення, той може поекспериментувати з іншими значеннями rmin та rmax.

P. S. Це оновлена версія моделі Шеллінга, де агенти мають різні значення порогів. Попередня версія коду моделі містила помилку:(