Прогноз численности русских к 2050 году
Migration Tendencies Analysis Center
2025-11-05
Настройка окружения
Этот блок конфигурирует knitr и подключает используемые пакеты, создавая единое рабочее окружение для расчёта демографической модели.
Он не требует входных данных, но задаёт глобальные опции вывода и делает доступными библиотеки dplyr, readr, ggplot2 и другие, на которые опираются все последующие шаги
Горизонт моделирования и индексы
Здесь определяется временной горизонт прогнозирования, возрастная сетка и обозначения полов, формируя основу для всех многомерных массивов модели.
Ключевые выходные структуры: - years, ages,
sexes — последовательности лет, возрастов и кодов пола; -
idxF, idxM, colS — индексы и
утилита для обращения к соответствующим срезам матриц.
## =========================================================
## МОДЕЛЬ С ПОЛО-ВОЗРАСТНОЙ СТРУКТУРОЙ (1-год. возраста, 0..100; пол: F/M)
## Шаг 2: национальный уровень, M1 / M+ / русские / прочие коренные, конвергенция TFR
## =========================================================
## ---------------------------
## 0) Горизонт, сетки и утилиты
## ---------------------------
years <- 2021:2050; ny <- length(years) # модельный горизонт и число моделируемых лет
ages <- 0:100; na <- length(ages) # возрастная сетка (0..100) и её длина
sexes <- c("F","M"); ns <- length(sexes) # обозначения полов и их количество
idxF <- 1; idxM <- 2 # индексы для обращения к женским/мужским столбцам в матрицах
colS <- function(s) ifelse(s == "F", idxF, idxM) # утилита: возвращает индекс пола по буквенному кодуСвод параметров модели
Блок формирует список model_inputs с основными
параметрами населения, миграции, рождаемости и сценарными
коэффициентами, которые используются в расчётах.
На выходе получается иерархическая структура с численностями базового года, характеристиками мигрантов, настройками второго поколения и траекториями сценариев.
model_inputs <- list(
population = list(
total_2021 = 147590600,
male_2021 = 68725249,
female_2021 = 78865351
),
mplus = list(
total_nominal_2021 = 12631000,
include_base = FALSE
),
m1 = list(
target_2021 = 6895947,
male_share_2021 = 0.55,
age_min = 18,
age_max = 49,
mu_male = 27,
sigma_male = 5,
mu_female = 26,
sigma_female = 6
),
sg2 = list(
total_2021 = 1.05e6,
cohort_shares = c(`2001-2005` = 0.05, `2006-2010` = 0.16, `2011-2015` = 0.38, `2016-2020` = 0.41),
srb = 105,
share_both_migrants = 0.32,
share_mixed_mother = 0.33,
share_mixed_father = 0.35,
origin_shares = c(
Kazakhstan = 0.29,
CentralAsia_exclKZ = 0.36,
Ukraine = 0.15,
Caucasus = 0.10,
Other = 0.10
)
),
scenario = list(
TFR_nat = seq(1.50, 1.65, length.out = ny),
mort_impr = seq(1.00, 0.85, length.out = ny),
nm_scenario = toupper(getOption("nm_scenario", Sys.getenv("NM_SCENARIO", "MMM"))),
nm_2021 = getOption("nm_2021", default = 150000),
nm_mid_peak = 250000,
nm_high_peak = 430000,
nm_low_peak = 90000,
nm_plateau_year = 2035,
nm_drift_start = 2041,
nm_drift_high = 0.005,
nm_drift_low = -0.005,
rho = rep(0.11, ny),
kappa_nat = rep(0.05, ny),
lambda = seq(1, 0, length.out = ny),
alpha_assim = rep(0.40, ny),
beta_M1_core = rep(0.50, ny),
m_M1_start = 1.15,
m_M1_target = 1.00,
m_Mplus_start = 1.08,
m_Mplus_target = 1.00,
m_RUS_start = 0.95,
m_RUS_target = 1.00,
m_CORE_start = 0.95,
m_CORE_target = 1.00
)
)Вспомогательные функции расчёта
В этом разделе определяются функции нормировки, интерполяции, проверки размерностей и другие утилиты, облегчающие обработку данных и выполнение демографических операций.
Каждая функция не принимает внешних данных напрямую, но предоставляется глобально для повторного использования в дальнейшем коде модели.
## ---------------------------
## ВСПОМОГАТЕЛЬНЫЕ ФУНКЦИИ
## ---------------------------
norm1 <- function(x) if (sum(x) == 0) x else x / sum(x)
# См. таблицу Function catalog.
blend <- function(start, target, lam) {
if (!length(start) || !length(target)) {
stop("Параметры start и target не должны быть пустыми")
}
target + (start - target) * lam
}
# См. таблицу Function catalog.
nm_path <- function(years,
scenario = c("MMM", "HHH", "LLL"),
nm_start,
mid_peak,
high_peak,
low_peak,
plateau_year,
drift_start,
drift_high,
drift_low,
clamp0 = TRUE) {
scenario <- match.arg(toupper(scenario), c("MMM", "HHH", "LLL"))
if (!length(years)) {
stop("Список лет для расчёта траектории миграции пуст")
}
y0 <- if (2021 %in% years) 2021 else min(years)
target <- switch(
scenario,
"MMM" = mid_peak,
"HHH" = high_peak,
"LLL" = low_peak
)
ramp_years <- years[years <= plateau_year]
post_years <- years[years > plateau_year]
ramp <- approx(
x = c(y0, plateau_year),
y = c(nm_start, target),
xout = ramp_years,
rule = 2
)$y
post <- rep(target, length(post_years))
if (length(post) && any(post_years >= drift_start)) {
idx <- which(post_years >= drift_start)
t <- seq_along(idx) - 1
rate <- switch(
scenario,
"MMM" = 0,
"HHH" = drift_high,
"LLL" = drift_low
)
post[idx] <- target * (1 + rate)^t
}
nm <- c(ramp, post)
if (clamp0) nm <- pmax(0, nm)
as.numeric(nm)
}
# См. таблицу Function catalog.
scenario_defaults <- function(years, overrides = list(), assumptions = model_inputs$scenario) {
ny <- length(years)
scenario <- modifyList(assumptions, overrides)
# применяем пользовательские переопределения (если заданы)
if (is.null(overrides$net_mig)) {
scenario$net_mig <- nm_path(
years = years,
scenario = scenario$nm_scenario,
nm_start = scenario$nm_2021,
mid_peak = scenario$nm_mid_peak,
high_peak = scenario$nm_high_peak,
low_peak = scenario$nm_low_peak,
plateau_year = scenario$nm_plateau_year,
drift_start = scenario$nm_drift_start,
drift_high = scenario$nm_drift_high,
drift_low = scenario$nm_drift_low
)
}
if (length(scenario$net_mig) != ny) {
stop(sprintf(
"Вектор net_mig имеет длину %d, ожидалось %d",
length(scenario$net_mig),
ny
))
}
scenario$net_mig_paths <- list(
LLL = nm_path(
years = years,
scenario = "LLL",
nm_start = scenario$nm_2021,
mid_peak = scenario$nm_mid_peak,
high_peak = scenario$nm_high_peak,
low_peak = scenario$nm_low_peak,
plateau_year = scenario$nm_plateau_year,
drift_start = scenario$nm_drift_start,
drift_high = scenario$nm_drift_high,
drift_low = scenario$nm_drift_low
),
MMM = nm_path(
years = years,
scenario = "MMM",
nm_start = scenario$nm_2021,
mid_peak = scenario$nm_mid_peak,
high_peak = scenario$nm_high_peak,
low_peak = scenario$nm_low_peak,
plateau_year = scenario$nm_plateau_year,
drift_start = scenario$nm_drift_start,
drift_high = scenario$nm_drift_high,
drift_low = scenario$nm_drift_low
),
HHH = nm_path(
years = years,
scenario = "HHH",
nm_start = scenario$nm_2021,
mid_peak = scenario$nm_mid_peak,
high_peak = scenario$nm_high_peak,
low_peak = scenario$nm_low_peak,
plateau_year = scenario$nm_plateau_year,
drift_start = scenario$nm_drift_start,
drift_high = scenario$nm_drift_high,
drift_low = scenario$nm_drift_low
)
)
if (is.null(scenario$m_M1)) {
scenario$m_M1 <- blend(scenario$m_M1_start, scenario$m_M1_target, scenario$lambda)
}
if (is.null(scenario$m_Mplus)) {
scenario$m_Mplus <- blend(scenario$m_Mplus_start, scenario$m_Mplus_target, scenario$lambda)
}
if (is.null(scenario$m_RUS)) {
scenario$m_RUS <- blend(scenario$m_RUS_start, scenario$m_RUS_target, scenario$lambda)
}
if (is.null(scenario$m_CORE)) {
scenario$m_CORE <- blend(scenario$m_CORE_start, scenario$m_CORE_target, scenario$lambda)
}
scenario
}
# См. таблицу Function catalog.
age_up <- function(A, Sur) {
A <- as.matrix(A)
Sur <- as.matrix(Sur)
if (all(dim(Sur) == rev(dim(A)))) {
Sur <- t(Sur)
}
if (is.null(dim(Sur)) || (any(dim(Sur) == 1) && !all(dim(Sur) == dim(A)))) {
Sur <- matrix(Sur, nrow = nrow(A), ncol = ncol(A), byrow = FALSE)
}
if (!all(dim(Sur) == dim(A))) {
stop(sprintf(
"Размерность Sur (%s) не совпадает с размерностью A (%s)",
paste(dim(Sur), collapse = "×"),
paste(dim(A), collapse = "×")
))
}
B <- A * Sur
res <- matrix(0, nrow = nrow(B), ncol = ncol(B))
if (nrow(B) > 1) {
res[2:nrow(B), ] <- B[1:(nrow(B)-1), ]
}
res[nrow(B), ] <- res[nrow(B), ] + B[nrow(B), ]
res[1, ] <- 0
res
}
# См. таблицу Function catalog.
check_tensor_shape <- function(arr, expected_dim, name) {
actual <- dim(arr)
if (!identical(actual, expected_dim)) {
stop(sprintf(
"%s имеет размерность %s, ожидалась %s",
name,
paste(actual, collapse = "×"),
paste(expected_dim, collapse = "×")
))
}
invisible(arr)
}
# См. таблицу Function catalog.
make_trunc_gauss <- function(ages, mu, sigma, min_age, max_age) {
weights <- exp(-0.5 * ((ages - mu) / sigma)^2)
weights[ages < min_age | ages > max_age] <- 0
total <- sum(weights)
if (total <= 0) {
stop("Профиль мигрантов обнулился: скорректируйте параметры mu/sigma или возрастные границы")
}
weights / total
}
# См. таблицу Function catalog.
allocate_migrants <- function(pop_matrix, profile_matrix, target, tol = 1e-6, maxit = 1000) {
profile_matrix <- pmax(profile_matrix, 0)
profile_sum <- sum(profile_matrix)
if (profile_sum <= 0) {
stop("Профиль распределения мигрантов имеет нулевую сумму")
}
profile_matrix <- profile_matrix / profile_sum
alloc <- profile_matrix * target
for (iter in seq_len(maxit)) {
overfull <- alloc > pop_matrix
if (any(overfull)) {
alloc[overfull] <- pop_matrix[overfull]
}
gap <- target - sum(alloc)
if (abs(gap) < tol) break
capacity <- pmax(pop_matrix - alloc, 0)
capacity_weights <- capacity * profile_matrix
cap_sum <- sum(capacity_weights)
if (cap_sum <= 0) {
alloc <- pmin(alloc, pop_matrix)
alloc <- alloc * (target / sum(alloc))
next
}
alloc <- alloc + gap * (capacity_weights / cap_sum)
alloc <- pmax(alloc, 0)
}
list(M_as = alloc, N_as = pmax(pop_matrix - alloc, 0))
}
# См. таблицу Function catalog.
build_age_weights_from_cohorts <- function(ages, cohort_shares, year_ref) {
if (!length(cohort_shares)) {
stop("Не заданы доли когорт для построения возрастного профиля")
}
weights <- numeric(length(ages))
names(weights) <- ages
norm_shares <- cohort_shares / sum(cohort_shares)
for (label in names(norm_shares)) {
bounds <- strsplit(gsub("\\s", "", label), "-", fixed = TRUE)[[1]]
if (grepl("+", label, fixed = TRUE)) {
start_year <- suppressWarnings(as.integer(gsub("[^0-9]", "", label)))
end_year <- year_ref
} else if (length(bounds) == 2) {
start_year <- suppressWarnings(as.integer(bounds[1]))
end_year <- suppressWarnings(as.integer(bounds[2]))
} else {
start_year <- suppressWarnings(as.integer(bounds[1]))
end_year <- start_year
}
if (is.na(start_year) || is.na(end_year)) next
cohort_years <- seq.int(start_year, end_year)
share_per_year <- norm_shares[[label]] / length(cohort_years)
for (birth_year in cohort_years) {
age <- year_ref - birth_year
if (as.character(age) %in% names(weights)) {
weights[as.character(age)] <- weights[as.character(age)] + share_per_year
}
}
}
if (sum(weights) <= 0) {
stop("Возрастной профиль второго поколения обнулился: проверьте настройки когорт")
}
weights / sum(weights)
}
# См. таблицу Function catalog.
distribute_sg2_details <- function(base_matrix, union_shares, origin_shares) {
if (sum(base_matrix) == 0) {
empty_union <- array(0, dim = c(nrow(base_matrix), ncol(base_matrix), 0))
empty_origin <- array(0, dim = c(nrow(base_matrix), ncol(base_matrix), 0, 0))
return(list(union = empty_union, origin = empty_origin))
}
norm_union <- norm1(pmax(union_shares, 0))
norm_origin <- norm1(pmax(origin_shares, 0))
union_arr <- array(0,
dim = c(nrow(base_matrix), ncol(base_matrix), length(norm_union)),
dimnames = list(age = rownames(base_matrix), sex = colnames(base_matrix), union = names(norm_union))
)
for (i in seq_along(norm_union)) {
union_arr[,, i] <- base_matrix * norm_union[[i]]
}
origin_arr <- array(0,
dim = c(nrow(base_matrix), ncol(base_matrix), length(norm_union), length(norm_origin)),
dimnames = list(
age = rownames(base_matrix),
sex = colnames(base_matrix),
union = names(norm_union),
origin = names(norm_origin)
)
)
for (j in seq_along(norm_origin)) {
origin_arr[,,, j] <- union_arr * norm_origin[[j]]
}
list(union = union_arr, origin = origin_arr)
}
# См. таблицу Function catalog.
expand_age_span <- function(label, max_age) {
lbl <- trimws(label)
lbl <- gsub("\u2013", "-", lbl, fixed = TRUE) # en dash → hyphen
lbl <- gsub("\u2014", "-", lbl, fixed = TRUE) # em dash → hyphen
lbl <- gsub("[[:space:]]*-[[:space:]]*", "-", lbl) # убираем пробелы вокруг дефиса
if (grepl("\\+", lbl)) {
start <- suppressWarnings(as.integer(gsub("[^0-9]", "", lbl)))
end <- max_age
} else if (grepl("-", lbl, fixed = TRUE)) {
parts <- strsplit(lbl, "-", fixed = TRUE)[[1]]
start <- suppressWarnings(as.integer(parts[1]))
end <- suppressWarnings(as.integer(parts[length(parts)]))
} else {
start <- suppressWarnings(as.integer(lbl))
end <- start
}
if (is.na(start) || is.na(end)) {
return(integer(0))
}
start <- max(0, start)
end <- min(max_age, end)
if (end < start) {
return(integer(0))
}
seq.int(start, end)
}
# См. таблицу Function catalog.
clean_numeric_value <- function(x) {
if (is.na(x) || trimws(x) %in% c("", "-")) return(0)
val <- gsub("[\u00A0[:space:]]", "", as.character(x))
val <- gsub(",", "", val, fixed = TRUE)
suppressWarnings(as.numeric(val))
}
# См. таблицу Function catalog.
sanitize_names <- function(nms) {
clean <- tolower(trimws(as.character(nms)))
clean <- gsub("[^[:alnum:]_]+", "_", clean)
empty <- is.na(clean) | clean == ""
if (any(empty)) {
clean[empty] <- paste0("col", which(empty))
}
make.unique(clean, sep = "_")
}
# См. таблицу Function catalog.
load_age_sex_distribution <- function(path, label, max_age) {
if (!file.exists(path)) {
return(NULL)
}
readers <- list(
function() readr::read_csv(
path,
col_types = readr::cols(.default = readr::col_character()),
trim_ws = TRUE
),
function() readr::read_tsv(
path,
col_types = readr::cols(.default = readr::col_character()),
trim_ws = TRUE
),
function() readr::read_delim(
path,
delim = ";",
col_types = readr::cols(.default = readr::col_character()),
trim_ws = TRUE
)
)
raw <- NULL
for (reader in readers) {
raw <- tryCatch(reader(), error = function(e) NULL)
if (!is.null(raw)) break
}
if (is.null(raw)) {
warning(sprintf("Не удалось прочитать %s", path))
return(NULL)
}
names(raw) <- sanitize_names(names(raw))
age_col <- names(raw)[1]
female_col <- names(raw)[grepl("female|famale|жен", names(raw))][1]
male_col <- names(raw)[grepl("male|men|муж", names(raw))][1]
if (is.na(female_col) || is.na(male_col)) {
warning(sprintf(
"%s должен содержать столбцы с численностью женщин и мужчин",
basename(path)
))
return(NULL)
}
female_vec <- numeric(max_age + 1)
male_vec <- numeric(max_age + 1)
for (i in seq_len(nrow(raw))) {
label_i <- trimws(raw[[age_col]][i])
if (is.na(label_i) || label_i == "") next
ages_span <- expand_age_span(label_i, max_age)
if (!length(ages_span)) next
f_val <- clean_numeric_value(raw[[female_col]][i])
m_val <- clean_numeric_value(raw[[male_col]][i])
if (length(ages_span) > 0 && (f_val > 0 || m_val > 0)) {
female_vec[ages_span + 1] <- female_vec[ages_span + 1] + f_val / length(ages_span)
male_vec[ages_span + 1] <- male_vec[ages_span + 1] + m_val / length(ages_span)
}
}
total <- sum(female_vec) + sum(male_vec)
if (total <= 0) {
warning(sprintf("В %s нет ненулевых значений для %s", basename(path), label))
return(NULL)
}
list(female = female_vec, male = male_vec, total = total)
}
# См. таблицу Function catalog.
load_tfr_status <- function(path) {
if (!file.exists(path)) {
return(NULL)
}
possible_delims <- c(",", ";", "|", "\t")
tfr_raw <- NULL
for (delim in possible_delims) {
tfr_raw <- tryCatch(
readr::read_delim(
path,
delim = delim,
trim_ws = TRUE,
show_col_types = FALSE,
locale = readr::locale(decimal_mark = ".", grouping_mark = ",")
),
error = function(e) NULL
)
if (!is.null(tfr_raw)) break
}
if (is.null(tfr_raw)) {
warning(sprintf("Не удалось прочитать %s — данные по статусам не будут использованы", path))
return(NULL)
}
names(tfr_raw) <- sanitize_names(names(tfr_raw))
status_col <- names(tfr_raw)[grepl("статус", names(tfr_raw), ignore.case = TRUE)]
tfr_col <- names(tfr_raw)[grepl("tfr", names(tfr_raw), ignore.case = TRUE)]
ethnos_col <- names(tfr_raw)[grepl("этнос", names(tfr_raw), ignore.case = TRUE)]
if (length(status_col) == 0) {
warning(sprintf("В файле %s отсутствует столбец со статусом — используем значения по умолчанию", path))
return(list(raw = tfr_raw, status_by_ethnos = tibble(), status_summary = tibble(), lookup = c()))
}
tfr_tbl <- tfr_raw %>%
mutate(
status = tolower(trimws(as.character(.data[[status_col[1]]]))),
ethnos = if (length(ethnos_col) > 0) trimws(as.character(.data[[ethnos_col[1]]])) else NA_character_,
tfr_value = if (length(tfr_col) > 0) readr::parse_number(as.character(.data[[tfr_col[1]]])) else NA_real_
) %>%
filter(!is.na(status) & status != "")
status_summary <- tfr_tbl %>%
filter(!is.na(tfr_value)) %>%
group_by(status) %>%
summarise(tfr = mean(tfr_value, na.rm = TRUE), .groups = "drop")
status_map <- c(
russian = "русские",
other_core = "другие коренные народы",
migrant = "народы-мигранты"
)
tfr_lookup <- setNames(rep(NA_real_, length(status_map)), names(status_map))
for (nm in names(status_map)) {
matched <- status_summary %>% filter(status == status_map[[nm]])
if (nrow(matched) > 0) {
tfr_lookup[[nm]] <- matched$tfr[1]
}
}
list(
raw = tfr_tbl,
status_by_ethnos = tfr_tbl %>% select(ethnos, status) %>% filter(!is.na(ethnos) & ethnos != ""),
status_summary = status_summary,
lookup = tfr_lookup
)
}
# См. таблицу Function catalog.
load_inputs <- function(inputs) {
if (!length(inputs)) {
return(list(errors = list()))
}
format_message <- function(template, default_message, values) {
safe_call <- function(fn, args, fallback) {
tryCatch(do.call(fn, args), error = function(e) fallback)
}
if (is.function(template)) {
return(safe_call(template, values, default_message))
}
if (is.character(template)) {
return(safe_call(sprintf, c(list(template), values), default_message))
}
default_message
}
notify <- function(handler, msg) {
if (is.function(handler)) {
handler(msg)
} else {
warning(msg)
}
}
results <- list()
errors <- list()
for (name in names(inputs)) {
spec <- inputs[[name]]
if (is.null(spec$loader) || !is.function(spec$loader)) {
stop(sprintf("Для источника \"%s\" требуется функция-загрузчик", name))
}
path <- spec$path
label <- if (!is.null(spec$label)) spec$label else name
missing_message <- spec$missing_message
missing_handler <- if (!is.null(spec$notify_missing)) spec$notify_missing else warning
error_message <- spec$error_message
error_handler <- if (!is.null(spec$notify_error)) spec$notify_error else warning
if (is.null(path) || is.na(path) || !nzchar(path)) {
default_msg <- sprintf("Файл для \"%s\" не указан", label)
msg <- format_message(
missing_message,
default_msg,
list(label)
)
notify(missing_handler, msg)
errors[[name]] <- msg
results[[name]] <- NULL
next
}
if (!file.exists(path)) {
default_msg <- sprintf("Файл %s не найден", path)
msg <- format_message(
missing_message,
default_msg,
list(path, label)
)
notify(missing_handler, msg)
errors[[name]] <- msg
results[[name]] <- NULL
next
}
res <- tryCatch(spec$loader(path), error = function(e) e)
if (inherits(res, "error")) {
default_msg <- sprintf("Ошибка при обработке файла %s: %s", path, conditionMessage(res))
msg <- format_message(
error_message,
default_msg,
list(path = path, label = label, message = conditionMessage(res))
)
notify(error_handler, msg)
errors[[name]] <- msg
results[[name]] <- NULL
} else {
results[[name]] <- res
}
}
results$errors <- errors
results
}
# См. таблицу Function catalog.
compute_mortality_improvement <- function(target_years, le_table) {
required_cols <- c("year", "female", "male")
if (!all(required_cols %in% names(le_table))) {
stop(sprintf(
"Таблица продолжительности жизни должна содержать столбцы: %s",
paste(required_cols, collapse = ", ")
))
}
if (!length(target_years)) {
stop("Список лет для расчёта улучшения смертности пуст")
}
female_vals <- le_table$female
male_vals <- le_table$male
years_src <- le_table$year
if (any(!is.finite(female_vals)) || any(!is.finite(male_vals))) {
stop("В таблице продолжительности жизни присутствуют некорректные значения")
}
interp_f <- approx(years_src, female_vals, xout = target_years, rule = 2, ties = mean)$y
interp_m <- approx(years_src, male_vals, xout = target_years, rule = 2, ties = mean)$y
base_year <- min(target_years)
base_f <- interp_f[target_years == base_year][1]
base_m <- interp_m[target_years == base_year][1]
if (!is.finite(base_f) || base_f <= 0 || !is.finite(base_m) || base_m <= 0) {
stop("Базовые значения ожидаемой продолжительности жизни должны быть положительными")
}
female_factor <- base_f / interp_f
male_factor <- base_m / interp_m
female_factor[target_years == base_year] <- 1
male_factor[target_years == base_year] <- 1
combined <- pmax(0, (female_factor + male_factor) / 2)
combined[target_years == base_year] <- 1
list(
combined = combined,
female = female_factor,
male = male_factor,
life_expectancy = tibble(
year = target_years,
female = interp_f,
male = interp_m
)
)
}
# См. таблицу Function catalog.
make_survival_matrix <- function(q_t) {
q_mat <- as.matrix(q_t)
matrix(
pmax(0, 1 - as.vector(q_mat)),
nrow = nrow(q_mat),
ncol = ncol(q_mat),
byrow = FALSE
)
}
# См. таблицу Function catalog.
advance_with_mortality <- function(states, Sur) {
lapply(states, age_up, Sur = Sur)
}
# См. таблицу Function catalog.
apply_migration <- function(M1_matrix, net_mig_t, mig_prof) {
check_tensor_shape(mig_prof, dim(M1_matrix), "mig_prof")
M1_matrix + net_mig_t * mig_prof
}
# См. таблицу Function catalog.
compute_births <- function(states, asfr_shape, TFR_nat_t, multipliers, rho_t,
fert_age_mask, srb, idxF, idxM) {
group_names <- names(states)
asfr_nat_t <- asfr_shape * TFR_nat_t
B_raw <- setNames(numeric(length(group_names)), group_names)
for (grp in group_names) {
mult <- if (grp %in% names(multipliers)) multipliers[[grp]] else NULL
if (is.null(mult) || is.na(mult)) mult <- 1
asfr_grp <- asfr_nat_t * mult
B_raw[[grp]] <- sum(states[[grp]][fert_age_mask, idxF] * asfr_grp[fert_age_mask])
}
total_fertile <- Reduce(
`+`,
lapply(states, function(mat) mat[fert_age_mask, idxF])
)
B_nat_target <- sum(total_fertile * asfr_nat_t[fert_age_mask])
total_raw <- sum(B_raw)
sf <- if (total_raw > 0) B_nat_target / total_raw else 1
B_scaled <- B_raw * sf
numM <- sum(states$M1[fert_age_mask, idxM] + states$Mpl[fert_age_mask, idxM])
denM <- sum(
Reduce(`+`, lapply(states, function(mat) mat[fert_age_mask, idxM]))
)
share_M_male <- ifelse(denM == 0, 0, numM / denM)
core_groups <- setdiff(group_names, c("M1", "Mpl"))
core_mixed_share <- (1 - rho_t) * share_M_male
if (!length(core_groups)) {
core_mixed_share <- 0
}
core_births_total <- sum(B_scaled[core_groups])
core_to_mplus <- core_mixed_share * core_births_total
mplus_newborns <- (B_scaled[["M1"]] + B_scaled[["Mpl"]] + core_to_mplus)
states$Mpl[1, idxM] <- states$Mpl[1, idxM] + mplus_newborns * srb["m"]
states$Mpl[1, idxF] <- states$Mpl[1, idxF] + mplus_newborns * srb["f"]
for (grp in core_groups) {
newborns <- B_scaled[[grp]] * (1 - core_mixed_share)
states[[grp]][1, idxM] <- states[[grp]][1, idxM] + newborns * srb["m"]
states[[grp]][1, idxF] <- states[[grp]][1, idxF] + newborns * srb["f"]
}
list(
states = states,
births = list(
components = B_scaled,
total = sum(B_scaled)
)
)
}
# См. таблицу Function catalog.
update_identity_shares <- function(states, alpha, beta, fallback_share_rus_core) {
totals <- vapply(states, sum, numeric(1))
pop_total <- sum(totals)
if (pop_total <= 0) {
return(c(RUS = NA_real_, IND = NA_real_, MIGID = NA_real_))
}
rus_groups <- intersect(names(states), c("RUS", "RUS_UND"))
other_core_groups <- setdiff(names(states), c("M1", "Mpl", rus_groups))
core_base_total <- sum(totals[c(rus_groups, other_core_groups)])
fallback <- ifelse(is.na(fallback_share_rus_core), 0.5, fallback_share_rus_core)
share_rus_core <- if (core_base_total > 0) {
sum(totals[rus_groups]) / core_base_total
} else {
fallback
}
core_from_Mplus <- alpha * (totals[["Mpl"]] / pop_total)
core_from_M1 <- beta * (totals[["M1"]] / pop_total)
core_total_id <- core_from_Mplus + core_from_M1 + (core_base_total / pop_total)
id_RUS_share <- core_total_id * share_rus_core
id_IND_share <- core_total_id * (1 - share_rus_core)
id_MIGID_share <- 1 - (id_RUS_share + id_IND_share)
c(RUS = id_RUS_share, IND = id_IND_share, MIGID = id_MIGID_share)
}
# См. таблицу Function catalog.Каталог определённых функций
Раздел собирает справочный tibble function_catalog,
который описывает назначение вспомогательных функций и упрощает
навигацию по коду.
На выходе формируется таблица с именами и текстовыми пояснениями, пригодная для быстрой проверки доступных инструментов.
function_catalog <- tribble(
~name, ~description,
"norm1", "Нормирует вектор/матрицу к сумме 1, сохраняя форму аргумента",
"blend", "Линейное сглаживание между стартовым и целевым значением по фактору lambda (0..1)",
"nm_path", "Строит траекторию чистой миграции по сценариям Юмагузина (LLL/MMM/HHH)",
"scenario_defaults", "Возвращает структуру сценария с возможностью переопределения параметров",
"age_up", "Переводит матрицу численности на следующий возраст, применяя выживаемость",
"check_tensor_shape", "Проверяет, что тензор имеет ожидаемую размерность",
"make_trunc_gauss", "Строит усечённое нормальное распределение по возрастам",
"allocate_migrants", "Распределяет целевой поток мигрантов с учётом профиля и вместимости",
"build_age_weights_from_cohorts", "Формирует однолеточный возрастной профиль из долей когорт",
"distribute_sg2_details", "Раскладывает массив SG2 по типу союза и происхождению",
"expand_age_span", "Преобразует текстовый возрастной интервал в набор однолетних возрастов",
"clean_numeric_value", "Парсит числовые значения, очищая пробелы и разделители",
"sanitize_names", "Очищает и унифицирует имена столбцов",
"load_age_sex_distribution", "Загружает половозрастное распределение и разворачивает его по возрастам",
"load_tfr_status", "Загружает таблицу TFR по национальностям и возвращает статусы",
"compute_mortality_improvement", "Строит траекторию снижения смертности из прогноза продолжительности жизни",
"make_survival_matrix", "Преобразует вероятности смерти в матрицу выживания",
"advance_with_mortality", "Применяет старение ко всем демографическим состояниям",
"apply_migration", "Добавляет чистую миграцию в массив M1 согласно профилю",
"compute_births", "Рассчитывает рождаемость и распределяет новорождённых по группам",
"update_identity_shares", "Вычисляет доли самоидентификации на конец года"
)
function_catalog## # A tibble: 21 × 2
## name description
## <chr> <chr>
## 1 norm1 Нормирует вектор/матрицу к сумме 1, сохраняя …
## 2 blend Линейное сглаживание между стартовым и целевы…
## 3 nm_path Строит траекторию чистой миграции по сценария…
## 4 scenario_defaults Возвращает структуру сценария с возможностью …
## 5 age_up Переводит матрицу численности на следующий во…
## 6 check_tensor_shape Проверяет, что тензор имеет ожидаемую размерн…
## 7 make_trunc_gauss Строит усечённое нормальное распределение по …
## 8 allocate_migrants Распределяет целевой поток мигрантов с учётом…
## 9 build_age_weights_from_cohorts Формирует однолеточный возрастной профиль из …
## 10 distribute_sg2_details Раскладывает массив SG2 по типу союза и проис…
## # ℹ 11 more rowsКонфигурация источников данных
Этот блок определяет базовую папку ввода и формирует словарь путей
input_files для всех требуемых CSV- и TSV-файлов.
Он не загружает данные напрямую, но подготавливает структуру путей и сообщений об ошибках, которая используется загрузчиками ниже.
if (!exists("input_dir", inherits = FALSE)) {
input_dir <- tryCatch({
current <- knitr::current_input()
if (is.null(current)) getwd() else dirname(normalizePath(current))
}, error = function(e) getwd())
}
input_files <- list(
rus_profile = list(
path = file.path(input_dir, "age_gender_nationality.csv"),
label = "распределения русских",
loader = function(path) load_age_sex_distribution(path, "Русские", max(ages))
),
other_core_profile = list(
path = file.path(input_dir, "age_gender_nationality_other.csv"),
label = "распределения прочих коренных народов",
loader = function(path) load_age_sex_distribution(path, "Прочие коренные", max(ages))
),
migpast_profile = list(
path = file.path(input_dir, "age_gender_nationality_migpast.csv"),
label = "распределения мигрантских народов",
loader = function(path) load_age_sex_distribution(path, "Мигрантские народы", max(ages))
),
unknown_profile = list(
path = file.path(input_dir, "age_gender_nationality_unknown.csv"),
label = "распределения населения с неопределённой национальностью",
loader = function(path) load_age_sex_distribution(path, "Россияне-неопределённые", max(ages))
),
tfr_status = list(
path = file.path(input_dir, "TFR_nationality.csv"),
label = "таблицы TFR по национальностям",
loader = load_tfr_status,
missing_message = "Файл %s не найден — используются сценарные мультипликаторы TFR по умолчанию",
notify_missing = message,
error_message = function(path, label, message) message
),
age_gender = list(
path = file.path(input_dir, "age_gender_2021_final.csv"),
label = "пирамиды населения 2021 года",
loader = function(path) {
attempts <- list(
function() readr::read_csv(path, show_col_types = FALSE,
locale = readr::locale(grouping_mark = ",")),
function() readr::read_csv2(path, show_col_types = FALSE,
locale = readr::locale(decimal_mark = ",", grouping_mark = " ")),
function() readr::read_tsv(path, show_col_types = FALSE)
)
age_gender_raw <- NULL
for (reader in attempts) {
age_gender_raw <- try(reader(), silent = TRUE)
if (!inherits(age_gender_raw, "try-error")) {
break
}
}
if (inherits(age_gender_raw, "try-error") || is.null(age_gender_raw)) {
stop(sprintf("Не удалось прочитать файл %s в поддерживаемых форматах CSV/TSV", path))
}
names(age_gender_raw) <- sanitize_names(names(age_gender_raw))
female_col <- dplyr::case_when(
"female" %in% names(age_gender_raw) ~ "female",
"famale" %in% names(age_gender_raw) ~ "famale",
TRUE ~ NA_character_
)
male_col <- dplyr::case_when(
"male" %in% names(age_gender_raw) ~ "male",
"men" %in% names(age_gender_raw) ~ "men",
TRUE ~ NA_character_
)
if (is.na(female_col) || is.na(male_col) || !("age" %in% names(age_gender_raw))) {
stop(
sprintf(
"Файл %s не содержит обязательные столбцы 'age', 'female'/'famale', 'male'",
basename(path)
)
)
}
age_gender <- age_gender_raw %>%
transmute(
age = as.integer(age),
female = readr::parse_number(as.character(.data[[female_col]])),
male = readr::parse_number(as.character(.data[[male_col]]))
) %>%
arrange(age)
age_grid <- tibble(age = ages)
age_gender <- age_grid %>%
left_join(age_gender, by = "age")
if (anyNA(age_gender$female) || anyNA(age_gender$male)) {
stop(
sprintf(
"В %s отсутствуют данные для некоторых возрастов 0..100",
basename(path)
)
)
}
age_gender
},
missing_message = "Файл %s не найден. Используем сглаженные распределения как приближение.",
notify_missing = warning,
error_message = function(path, label, message) message
),
life_expectancy = list(
path = file.path(input_dir, "life_expectancy_prediction.csv"),
label = "прогноза продолжительности жизни",
loader = function(path) {
possible_delims <- c(",", ";", "\t", "|")
lifeexp_raw <- NULL
for (delim in possible_delims) {
lifeexp_raw <- tryCatch(
readr::read_delim(
path,
delim = delim,
trim_ws = TRUE,
show_col_types = FALSE
),
error = function(e) NULL
)
if (!is.null(lifeexp_raw)) break
}
if (is.null(lifeexp_raw)) {
stop(sprintf("Не удалось прочитать файл %s — используется сценарный ряд mort_impr по умолчанию", path))
}
names(lifeexp_raw) <- sanitize_names(names(lifeexp_raw))
required <- c("years", "female", "male")
if (!all(required %in% names(lifeexp_raw))) {
stop(sprintf(
"В файле %s отсутствуют необходимые столбцы (%s) — используется mort_impr из сценария",
basename(path),
paste(required, collapse = ", ")
))
}
lifeexp_tbl <- lifeexp_raw %>%
mutate(
year = suppressWarnings(as.numeric(.data[["years"]])),
female = readr::parse_number(as.character(.data[["female"]])),
male = readr::parse_number(as.character(.data[["male"]]))
) %>%
filter(!is.na(year) & is.finite(female) & is.finite(male)) %>%
arrange(year) %>%
select(year, female, male)
if (nrow(lifeexp_tbl) == 0) {
stop(sprintf("В таблице %s нет валидных строк — используется mort_impr из сценария", path))
}
lifeexp_tbl
},
missing_message = "Файл %s не найден — используется mort_impr по умолчанию",
notify_missing = message,
error_message = function(path, label, message) message
)
)Стартовая база населения 2021 года
Раздел читает входные таблицы, извлекает профиль населения и инициализирует матрицы для различных групп населения на стартовый год.
На выходе формируются массивы rus_2021_as,
core_2021_as, M1_2021_as,
Mpl_2021_as и другие заготовки, а также фиксируются
потенциальные ошибки загрузки.
## ---------------------------
## 1) БАЗА 2021
## ---------------------------
pop_total_2021 <- model_inputs$population$total_2021 # численность населения РФ на стартовый год
male_2021 <- model_inputs$population$male_2021 # численность мужчин в 2021
female_2021 <- model_inputs$population$female_2021 # численность женщин в 2021
## Происхождение (M+, русские, прочие коренные) на 2021 (возрастное распределение загрузим из файлов)
Mplus_2021 <- model_inputs$mplus$total_nominal_2021 # потомки мигрантов (номинальная оценка)
include_Mplus_base <- model_inputs$mplus$include_base # при TRUE в расчёт базы войдут бывшие мигранты
Mplus_2021_effective <- if (include_Mplus_base) Mplus_2021 else 0
inputs <- load_inputs(input_files)
input_errors <- inputs$errors
tfr_nat_path <- input_files$tfr_status$path
rus_profile <- inputs$rus_profile
other_core_profile <- inputs$other_core_profile
migpast_profile <- inputs$migpast_profile
unknown_profile <- inputs$unknown_profile
tfr_status_data <- inputs$tfr_status
rus_2021_as <- array(0, dim = c(na, ns))
core_2021_as <- array(0, dim = c(na, ns))
rus_unknown_2021_as <- array(0, dim = c(na, ns))
M1_2021_as <- array(0, dim = c(na, ns))
share_RUS_2021 <- NA_real_Профили и балансировка стартового распределения
Здесь создаются возрастно-половые профили для мигрантов и коренных групп, распределяется вместимость и выравнивается база населения относительно контрольных итогов 2021 года.
Ключевые результаты включают профили mig_profile_target,
sg2_profile, детали sg2_details, а также
скорректированные массивы русских, коренных народов и неопределённых
групп.
## ---------------------------
## 2) ТЕСТОВЫЕ ПРОФИЛИ ПО ВОЗРАСТУ (для разложения базы)
## (здесь — реалистичные заглушки; подставишь реальную пирамиду ВПН-2021)
## ---------------------------
## Возрастные «веса» (F/M) — подгружаем фактическую пирамиду 2021 года из CSV
age_gender <- inputs$age_gender
if (!is.null(age_gender)) {
wF <- norm1(age_gender$female)
wM <- norm1(age_gender$male)
} else {
wF <- norm1(0.6 * dnorm(ages, mean = 34, sd = 18) + 0.4 * dnorm(ages, mean = 66, sd = 14))
wM <- norm1(0.6 * dnorm(ages, mean = 33, sd = 19) + 0.4 * dnorm(ages, mean = 64, sd = 15))
}
## Пирамида 2021 (масштабируем веса к муж/жен тоталам)
pop2021_weights <- array(0, dim = c(na, ns)) # стартовая возрастно-половая матрица населения (веса)
pop2021_weights[, idxF] <- female_2021 * wF
pop2021_weights[, idxM] <- male_2021 * wM
rownames(pop2021_weights) <- ages
colnames(pop2021_weights) <- sexes
P_as <- pop2021_weights
## Мигранты первого поколения: трудовой профиль и раскладка по переписи
M1_target_2021 <- model_inputs$m1$target_2021
M1_male_share_2021 <- model_inputs$m1$male_share_2021
M1_age_min <- model_inputs$m1$age_min
M1_age_max <- model_inputs$m1$age_max
M1_mu_m <- model_inputs$m1$mu_male
M1_sigma_m <- model_inputs$m1$sigma_male
M1_mu_f <- model_inputs$m1$mu_female
M1_sigma_f <- model_inputs$m1$sigma_female
age_profile_m <- make_trunc_gauss(ages, mu = M1_mu_m, sigma = M1_sigma_m, min_age = M1_age_min, max_age = M1_age_max)
age_profile_f <- make_trunc_gauss(ages, mu = M1_mu_f, sigma = M1_sigma_f, min_age = M1_age_min, max_age = M1_age_max)
mig_profile_target <- cbind(
(1 - M1_male_share_2021) * age_profile_f,
M1_male_share_2021 * age_profile_m
)
mig_profile_target <- mig_profile_target / sum(mig_profile_target)
mig_allocation <- allocate_migrants(P_as, mig_profile_target, M1_target_2021)
M1_2021_as <- mig_allocation$M_as
non_migrant_capacity <- mig_allocation$N_as
dimnames(M1_2021_as) <- dimnames(P_as)
dimnames(non_migrant_capacity) <- dimnames(P_as)
## Потомки мигрантов (второе поколение) по оценкам Смирновой
SG2_cfg <- model_inputs$sg2
sg2_age_weights <- build_age_weights_from_cohorts(ages, SG2_cfg$cohort_shares, year_ref = 2021)
sg2_male_share <- SG2_cfg$srb / (SG2_cfg$srb + 100)
sg2_female_share <- 1 - sg2_male_share
sg2_profile <- cbind(
sg2_age_weights * sg2_female_share,
sg2_age_weights * sg2_male_share
)
colnames(sg2_profile) <- sexes
rownames(sg2_profile) <- ages
sg2_alloc <- allocate_migrants(non_migrant_capacity, sg2_profile, SG2_cfg$total_2021)
Mpl_2021_as <- sg2_alloc$M_as
non_migrant_capacity <- sg2_alloc$N_as
dimnames(Mpl_2021_as) <- dimnames(P_as)
dimnames(non_migrant_capacity) <- dimnames(P_as)
sg2_details <- distribute_sg2_details(Mpl_2021_as, c(
both_migrants = SG2_cfg$share_both_migrants,
mixed_mother = SG2_cfg$share_mixed_mother,
mixed_father = SG2_cfg$share_mixed_father
), SG2_cfg$origin_shares)
zero_like <- function(mat) matrix(0, nrow = nrow(mat), ncol = ncol(mat))
profile_to_matrix <- function(profile) {
if (is.null(profile)) return(NULL)
cbind(profile$female, profile$male)
}
allocate_group_from_profile <- function(capacity_matrix, profile_matrix, label) {
result <- zero_like(capacity_matrix)
if (is.null(profile_matrix)) {
warning(sprintf("Не удалось загрузить распределение для %s — значения будут определены по остаточной вместимости", label))
return(list(alloc = result, capacity = capacity_matrix))
}
desired <- pmax(profile_matrix, 0)
desired_total <- sum(desired)
if (desired_total <= 0) {
warning(sprintf("В профиле %s отсутствуют ненулевые значения", label))
return(list(alloc = result, capacity = capacity_matrix))
}
cap_total <- sum(capacity_matrix)
if (cap_total <= 0) {
return(list(alloc = result, capacity = capacity_matrix))
}
if (desired_total > cap_total + 1e-6) {
desired <- desired * (cap_total / desired_total)
warning(sprintf("Численность %s превышает доступную вместимость и будет пропорционально уменьшена", label))
}
alloc <- pmin(desired, capacity_matrix)
list(alloc = alloc, capacity = pmax(capacity_matrix - alloc, 0))
}
capacity_matrix <- non_migrant_capacity
rus_res <- allocate_group_from_profile(capacity_matrix, profile_to_matrix(rus_profile), "русских")
rus_2021_as <- rus_res$alloc
capacity_matrix <- rus_res$capacity
if (is.null(dimnames(rus_2021_as))) {
dimnames(rus_2021_as) <- dimnames(P_as)
}
core_res <- allocate_group_from_profile(capacity_matrix, profile_to_matrix(other_core_profile), "прочих коренных народов")
core_2021_as <- core_res$alloc
capacity_matrix <- core_res$capacity
if (is.null(dimnames(core_2021_as))) {
dimnames(core_2021_as) <- dimnames(P_as)
}
unknown_res <- allocate_group_from_profile(capacity_matrix, profile_to_matrix(unknown_profile), "населения с неопределённой национальностью")## Warning in allocate_group_from_profile(capacity_matrix,
## profile_to_matrix(unknown_profile), : Численность населения с неопределённой
## национальностью превышает доступную вместимость и будет пропорционально
## уменьшенаrus_unknown_2021_as <- unknown_res$alloc
capacity_matrix <- unknown_res$capacity
if (is.null(dimnames(rus_unknown_2021_as))) {
dimnames(rus_unknown_2021_as) <- dimnames(P_as)
}
if (sum(capacity_matrix) > 0) {
rus_unknown_2021_as <- rus_unknown_2021_as + capacity_matrix
capacity_matrix[] <- 0
}
if (!exists("Mpl_2021_as", inherits = FALSE)) {
Mpl_2021_as <- array(0, dim = c(na, ns))
}
if (Mplus_2021_effective > 0) {
base_f <- if (!is.null(migpast_profile)) norm1(migpast_profile$female) else wF
base_m <- if (!is.null(migpast_profile)) norm1(migpast_profile$male) else wM
Mpl_2021_as[, idxF] <- Mpl_2021_as[, idxF] + Mplus_2021_effective * base_f
Mpl_2021_as[, idxM] <- Mpl_2021_as[, idxM] + Mplus_2021_effective * base_m
}
if (is.null(dimnames(Mpl_2021_as))) {
dimnames(Mpl_2021_as) <- dimnames(P_as)
}
## Остаток населения → «Россияне-неопределённые»
known_female <- sum(rus_2021_as[, idxF]) + sum(core_2021_as[, idxF]) +
sum(M1_2021_as[, idxF]) + sum(Mpl_2021_as[, idxF]) + sum(rus_unknown_2021_as[, idxF])
known_male <- sum(rus_2021_as[, idxM]) + sum(core_2021_as[, idxM]) +
sum(M1_2021_as[, idxM]) + sum(Mpl_2021_as[, idxM]) + sum(rus_unknown_2021_as[, idxM])
residual_female <- max(female_2021 - known_female, 0)
residual_male <- max(male_2021 - known_male, 0)
if (female_2021 + male_2021 - (known_female + known_male) < -1e-3) {
warning("Сумма численностей по национальным группам превышает общую численность населения. Проверьте входные данные.")
}
if (residual_female > 0) {
if (sum(rus_unknown_2021_as[, idxF]) > 0) {
scale_f <- (sum(rus_unknown_2021_as[, idxF]) + residual_female) /
sum(rus_unknown_2021_as[, idxF])
rus_unknown_2021_as[, idxF] <- rus_unknown_2021_as[, idxF] * scale_f
} else {
rus_unknown_2021_as[, idxF] <- residual_female * wF
}
}
if (residual_male > 0) {
if (sum(rus_unknown_2021_as[, idxM]) > 0) {
scale_m <- (sum(rus_unknown_2021_as[, idxM]) + residual_male) /
sum(rus_unknown_2021_as[, idxM])
rus_unknown_2021_as[, idxM] <- rus_unknown_2021_as[, idxM] * scale_m
} else {
rus_unknown_2021_as[, idxM] <- residual_male * wM
}
}
## При превышении суммарного населения корректируем «неопределённых» пропорционально
adj_female <- female_2021 - (sum(rus_2021_as[, idxF]) + sum(core_2021_as[, idxF]) +
sum(M1_2021_as[, idxF]) + sum(Mpl_2021_as[, idxF]) + sum(rus_unknown_2021_as[, idxF]))
adj_male <- male_2021 - (sum(rus_2021_as[, idxM]) + sum(core_2021_as[, idxM]) +
sum(M1_2021_as[, idxM]) + sum(Mpl_2021_as[, idxM]) + sum(rus_unknown_2021_as[, idxM]))
if (abs(adj_female) > 1) {
if (sum(rus_unknown_2021_as[, idxF]) > 0) {
rus_unknown_2021_as[, idxF] <- pmax(
rus_unknown_2021_as[, idxF] + adj_female * norm1(rus_unknown_2021_as[, idxF]),
0
)
} else if (adj_female > 0) {
rus_unknown_2021_as[, idxF] <- adj_female * wF
} else {
warning("Отрицательное отклонение женской численности не может быть распределено между 'Россиянами-неопределёнными'")
}
}
if (abs(adj_male) > 1) {
if (sum(rus_unknown_2021_as[, idxM]) > 0) {
rus_unknown_2021_as[, idxM] <- pmax(
rus_unknown_2021_as[, idxM] + adj_male * norm1(rus_unknown_2021_as[, idxM]),
0
)
} else if (adj_male > 0) {
rus_unknown_2021_as[, idxM] <- adj_male * wM
} else {
warning("Отрицательное отклонение мужской численности не может быть распределено между 'Россиянами-неопределёнными'")
}
}
rus_unknown_initial_sum <- sum(rus_unknown_2021_as)
rus_2021_as <- rus_2021_as + rus_unknown_2021_as
rus_unknown_2021_as <- matrix(0,
nrow = nrow(rus_2021_as),
ncol = ncol(rus_2021_as),
dimnames = dimnames(rus_2021_as)
)
## Итоговые суммы и контрольные соотношения
rus_sum <- sum(rus_2021_as)
core_sum <- sum(core_2021_as)
mig_sum <- sum(M1_2021_as)
rus_unknown_sum <- sum(rus_unknown_2021_as)
pop2021 <- M1_2021_as + Mpl_2021_as + rus_2021_as + core_2021_as + rus_unknown_2021_as
total_diff <- abs(sum(pop2021) - pop_total_2021)
stopifnot(total_diff < 1e-6 * max(1, pop_total_2021))
M1_share_sex <- c(F = 0.5, M = 0.5)
if (mig_sum > 0) {
M1_share_sex <- c(F = sum(M1_2021_as[, idxF]) / mig_sum,
M = sum(M1_2021_as[, idxM]) / mig_sum)
}
core_total <- rus_sum + core_sum + rus_unknown_sum
if (core_total > 0) {
share_RUS_2021 <- (rus_sum + rus_unknown_sum) / core_total
share_RUS_2021 <- min(max(share_RUS_2021, 0), 1)
}
msg_fmt <- scales::comma_format(accuracy = 1, big.mark = " ")
message(
sprintf(
paste0(
"База 2021: русские (включая неопределённую национальность) = %s, ",
"другие коренные = %s, мигранты = %s, бывшие мигранты = %s (номинально %s)"
),
msg_fmt(rus_sum),
msg_fmt(core_sum),
msg_fmt(mig_sum),
msg_fmt(sum(Mpl_2021_as)),
msg_fmt(Mplus_2021)
)
)## База 2021: русские (включая неопределённую национальность) = 125 608 418, другие коренные = 14 036 235, мигранты = 6 895 947, бывшие мигранты = 1 050 000 (номинально 12 631 000)if (rus_unknown_initial_sum > 0) {
message(
sprintf(
"В состав русских добавлено %s человек с неопределённой национальностью.",
msg_fmt(rus_unknown_initial_sum)
)
)
}## В состав русских добавлено 15 793 834 человек с неопределённой национальностью.if (pop_total_2021 > 0) {
pct_rus <- 100 * rus_sum / pop_total_2021
pct_core <- 100 * core_sum / pop_total_2021
pct_mig <- 100 * mig_sum / pop_total_2021
message(
sprintf(
"Доли 2021: русские = %.2f%%, другие коренные = %.2f%%, мигранты = %.2f%%",
pct_rus,
pct_core,
pct_mig
)
)
}## Доли 2021: русские = 85.11%, другие коренные = 9.51%, мигранты = 4.67%if (!is.null(sg2_details$union) && length(dim(sg2_details$union)) == 3) {
union_totals <- apply(sg2_details$union, 3, sum)
if (sum(union_totals) > 0) {
union_pct <- round(100 * union_totals / sum(union_totals), 1)
message(
sprintf(
"Структура союзов SG2: %s",
paste(sprintf("%s = %.1f%%", names(union_pct), union_pct), collapse = ", ")
)
)
}
}## Структура союзов SG2: both_migrants = 32.0%, mixed_mother = 33.0%, mixed_father = 35.0%if (!is.null(sg2_details$origin) && length(dim(sg2_details$origin)) == 4) {
origin_totals <- apply(sg2_details$origin, 4, sum)
if (sum(origin_totals) > 0) {
origin_pct <- round(100 * origin_totals / sum(origin_totals), 1)
message(
sprintf(
"Происхождение SG2: %s",
paste(sprintf("%s = %.1f%%", names(origin_pct), origin_pct), collapse = ", ")
)
)
}
}## Происхождение SG2: Kazakhstan = 29.0%, CentralAsia_exclKZ = 36.0%, Ukraine = 15.0%, Caucasus = 10.0%, Other = 10.0%Описательная статистика 2021
В этом разделе собираем ключевые срезы по возрастно-половым матрицам базового года, чтобы убедиться в согласованности исходных данных.
Далее приведены сводные таблицы и иллюстрации, которые помогают проверить баланс групп населения и выявить возможные расхождения с контрольной численностью 2021 года.
group_labels <- c(
rus = "Русские (включая неопределённую национальность)",
core = "Прочие коренные",
M1 = "Мигранты первого поколения",
Mpl = "Потомки мигрантов"
)
matrices_named <- list(
rus = get0("rus_2021_as", ifnotfound = NULL, inherits = FALSE),
core = get0("core_2021_as", ifnotfound = NULL, inherits = FALSE),
M1 = get0("M1_2021_as", ifnotfound = NULL, inherits = FALSE),
Mpl = get0("Mpl_2021_as", ifnotfound = NULL, inherits = FALSE)
)
matrices_tbls <- lapply(names(matrices_named), function(name) {
mat <- matrices_named[[name]]
if (is.null(mat) || length(dim(mat)) != 2) {
return(NULL)
}
df <- as.data.frame.table(mat, responseName = "pop", stringsAsFactors = FALSE)
if (!all(c("Var1", "Var2", "pop") %in% names(df))) {
return(NULL)
}
tibble(
age = as.integer(as.character(df$Var1)),
sex = as.character(df$Var2),
group = name,
pop = as.numeric(df$pop)
)
})
matrices_tbls <- Filter(Negate(is.null), matrices_tbls)
base_as_tbl <- bind_rows(matrices_tbls)
pop_total_2021 <- get0("pop_total_2021", ifnotfound = NA_real_, inherits = FALSE)
if (nrow(base_as_tbl) == 0) {
base_as_tbl <- tibble(age = integer(), sex = character(), group = character(), pop = numeric())
base_group_tbl <- tibble(group = character(), pop = numeric(), share_pct = numeric())
base_group_sex_tbl <- tibble(group = character(), sex = character(), pop = numeric(), share_pct = numeric())
pyramid_plot <- NULL
group_share_plot <- NULL
summary_points <- character()
} else {
base_as_tbl <- base_as_tbl %>%
mutate(
sex = factor(sex, levels = sexes, labels = c("Женщины", "Мужчины")),
group = factor(group, levels = names(matrices_named)),
group_label = group_labels[as.character(group)]
) %>%
arrange(group, sex, age)
pop_total_2021 <- sum(base_as_tbl$pop)
base_group_tbl <- base_as_tbl %>%
group_by(group, group_label) %>%
summarise(
pop = sum(pop),
share_pct = 100 * pop / pop_total_2021,
.groups = "drop"
)
base_group_sex_tbl <- base_as_tbl %>%
group_by(group, group_label, sex) %>%
summarise(
pop = sum(pop),
share_pct = 100 * pop / pop_total_2021,
.groups = "drop"
)
base_group_tbl <- base_group_tbl %>% filter(pop > 0)
positive_groups <- as.character(base_group_tbl$group)
base_group_sex_tbl <- base_group_sex_tbl %>% filter(group %in% positive_groups)
base_as_tbl <- base_as_tbl %>% filter(group %in% positive_groups)
pyramid_plot <- base_as_tbl %>%
mutate(
pop_signed = if_else(sex == "Мужчины", -pop, pop)
) %>%
ggplot(aes(x = age, y = pop_signed, fill = sex)) +
geom_col(width = 1) +
facet_wrap(~ group_label, ncol = 2, scales = "free_y") +
scale_y_continuous(
labels = function(x) scales::comma(abs(x), big.mark = " "),
name = "Численность"
) +
scale_fill_manual(values = c("Женщины" = "#f28e2b", "Мужчины" = "#4e79a7")) +
labs(x = "Возраст", fill = "Пол") +
theme_minimal()
group_share_plot <- base_group_tbl %>%
mutate(share = share_pct / 100) %>%
ggplot(aes(x = reorder(group_label, -share), y = share, fill = group_label)) +
geom_col(width = 0.7, show.legend = FALSE) +
scale_y_continuous(
labels = scales::percent_format(scale = 1),
expand = expansion(mult = c(0, 0.05))
) +
labs(x = "Группа", y = "Доля в населении 2021") +
theme_minimal()
diff_total <- if (!is.null(model_inputs$population$total_2021)) pop_total_2021 - model_inputs$population$total_2021 else NA_real_
diff_label <- dplyr::case_when(
is.na(diff_total) ~ "нет контрольной оценки для сравнения",
abs(diff_total) < 5 ~ "совпадает с контрольной численностью",
diff_total > 0 ~ paste0("на ", scales::comma(diff_total, accuracy = 1, big.mark = " "), " больше контрольной численности"),
TRUE ~ paste0("на ", scales::comma(abs(diff_total), accuracy = 1, big.mark = " "), " меньше контрольной численности")
)
dominant_group <- base_group_tbl %>% arrange(desc(pop)) %>% slice(1)
summary_points <- c(
sprintf(
"%s занимают %.1f%% стартовой базы.",
dominant_group$group_label,
dominant_group$share_pct
),
sprintf(
"Совокупная доля миграционных групп (M1 + Mpl) составляет %.1f%% населения.",
sum(base_group_tbl$share_pct[base_group_tbl$group %in% c("M1", "Mpl")], na.rm = TRUE)
),
sprintf("Сводная численность базы %s.", diff_label)
)
}
if (nrow(base_group_tbl) > 0) {
knitr::kable(
base_group_tbl %>%
mutate(
pop = scales::comma(pop, big.mark = " "),
share_pct = sprintf("%.1f", share_pct)
) %>%
select(`Группа` = group_label, `Численность` = pop, `Доля, %` = share_pct),
caption = "Итоги по группам населения базы 2021 года"
)
knitr::kable(
base_group_sex_tbl %>%
mutate(
pop = scales::comma(pop, big.mark = " "),
share_pct = sprintf("%.1f", share_pct)
) %>%
select(`Группа` = group_label, `Пол` = sex, `Численность` = pop, `Доля в населении, %` = share_pct),
caption = "Распределение по полу внутри групп"
)
} else {
message("Отсутствуют матрицы для построения описательной статистики.")
}| Группа | Пол | Численность | Доля в населении, % |
|---|---|---|---|
| Русские (включая неопределённую национальность) | Женщины | 67 263 331 | 45.6 |
| Русские (включая неопределённую национальность) | Мужчины | 58 345 087 | 39.5 |
| Прочие коренные | Женщины | 7 986 649 | 5.4 |
| Прочие коренные | Мужчины | 6 049 586 | 4.1 |
| Мигранты первого поколения | Женщины | 3 103 176 | 2.1 |
| Мигранты первого поколения | Мужчины | 3 792 771 | 2.6 |
| Потомки мигрантов | Женщины | 512 195 | 0.3 |
| Потомки мигрантов | Мужчины | 537 805 | 0.4 |
if (!is.null(pyramid_plot)) pyramid_plot
if (!is.null(group_share_plot)) group_share_plot
Краткие выводы по базе 2021
- Русские (включая неопределённую национальность) занимают 85.1% стартовой базы.
- Совокупная доля миграционных групп (M1 + Mpl) составляет 5.4% населения.
- Сводная численность базы совпадает с контрольной численностью.
Далее начинается моделирование сценариев.
Сценарные ряды рождаемости и миграции
Этот раздел формирует сценарий scenario, объединяя
значения по умолчанию и, при наличии данных, калибруя мультипликаторы
TFR и миграции.
На выходе создаётся обновлённый список с траекториями
net_mig, m_M1, m_Mplus,
m_RUS и m_CORE, а также вспомогательные
метаданные по входным таблицам.
## ---------------------------
## 3) СЦЕНАРНЫЕ РЯДЫ (MMM) — СОБРАНЫ В ОДНОЙ СТРУКТУРЕ
## ---------------------------
scenario <- scenario_defaults(years) # базовый набор сценарных рядов MMM
if (is.na(share_RUS_2021)) {
share_RUS_2021 <- 0.5
}
if (!is.null(tfr_status_data)) {
tfr_lookup <- tfr_status_data$lookup
status_summary <- tfr_status_data$status_summary
if (length(tfr_lookup) == 0) {
warning(sprintf("В %s отсутствуют агрегированные значения TFR по статусам — используются значения сценария по умолчанию", tfr_nat_path))
} else if (any(scenario$TFR_nat <= 0)) {
warning("Национальный TFR содержит неположительные значения — данные по статусам не применены")
} else {
if (!is.na(tfr_lookup[["migrant"]])) {
m1_path <- tfr_lookup[["migrant"]] / scenario$TFR_nat
scenario$m_M1_start <- m1_path[1]
scenario$m_M1 <- blend(
start = m1_path,
target = rep(scenario$m_M1_target, ny),
lam = scenario$lambda
)
}
if (!is.na(tfr_lookup[["russian"]])) {
rus_path <- tfr_lookup[["russian"]] / scenario$TFR_nat
scenario$m_RUS_start <- rus_path[1]
scenario$m_RUS <- blend(
start = rus_path,
target = rep(scenario$m_RUS_target, ny),
lam = scenario$lambda
)
}
if (!is.na(tfr_lookup[["other_core"]])) {
core_path <- tfr_lookup[["other_core"]] / scenario$TFR_nat
scenario$m_CORE_start <- core_path[1]
scenario$m_CORE <- blend(
start = core_path,
target = rep(scenario$m_CORE_target, ny),
lam = scenario$lambda
)
}
if (!is.na(tfr_lookup[["migrant"]]) || !is.na(tfr_lookup[["other_core"]])) {
base_core <- if (!is.na(tfr_lookup[["other_core"]])) tfr_lookup[["other_core"]] else tfr_lookup[["migrant"]]
mix_value <- dplyr::case_when(
!is.na(tfr_lookup[["migrant"]]) && !is.na(tfr_lookup[["other_core"]]) ~
(tfr_lookup[["migrant"]] + tfr_lookup[["other_core"]]) / 2,
!is.na(tfr_lookup[["migrant"]]) ~ tfr_lookup[["migrant"]],
TRUE ~ base_core
)
mpl_path <- mix_value / scenario$TFR_nat
scenario$m_Mplus_start <- mpl_path[1]
scenario$m_Mplus <- blend(
start = mpl_path,
target = rep(scenario$m_Mplus_target, ny),
lam = scenario$lambda
)
}
scenario$tfr_inputs <- list(
source = tfr_nat_path,
status_summary = status_summary,
lookup = tfr_lookup
)
}
} else {
message(sprintf("Файл %s не найден — используются сценарные мультипликаторы TFR по умолчанию", tfr_nat_path))
}Корректировка смертности на основе ожидаемой жизни
Блок загружает прогноз продолжительности жизни и преобразует его в коэффициенты улучшения смертности для сценария.
На выходе обновляются поля scenario$mort_impr и
scenario$mortality_profile, обеспечивая связь с исходной
таблицей.
## ---------------------------
## 3A) ПРОГНОЗ ПРОДОЛЖИТЕЛЬНОСТИ ЖИЗНИ → УЛУЧШЕНИЕ СМЕРТНОСТИ
## ---------------------------
lifeexp_path <- input_files$life_expectancy$path
lifeexp_tbl <- inputs$life_expectancy
if (!is.null(lifeexp_tbl)) {
lifeexp_profile <- compute_mortality_improvement(years, lifeexp_tbl)
scenario$mort_impr <- lifeexp_profile$combined
scenario$mortality_profile <- list(
source = lifeexp_path,
life_expectancy = lifeexp_profile$life_expectancy,
factors = tibble(
year = years,
mort_impr = lifeexp_profile$combined,
female_factor = lifeexp_profile$female,
male_factor = lifeexp_profile$male
)
)
}Подготовка данных по фертильности
Раздел загружает таблицу возрастной фертильности, проверяет наличие
ключевых столбцов и строит форму asfr_shape для
однолеточных возрастов.
Результатом становятся очищенные данные по паритетам рождения и набор
сценарных мультипликаторов m_M1, m_Mplus,
m_RUS, m_CORE для дальнейшего цикла.
## ======================
## БЛОК 4. ФЕРТИЛЬНОСТЬ (обновлённый)
## ======================
# (1) ФОРМА ASFR ИЗ ТВОЕЙ ГОТОВОЙ ТАБЛИЦЫ
# Формат таблицы: age_band, women_parity_0, ..., women_parity_5plus, women_total
if (!exists("input_dir", inherits = FALSE)) {
input_dir <- tryCatch({
current <- knitr::current_input()
if (is.null(current)) getwd() else dirname(normalizePath(current))
}, error = function(e) getwd())
}
path_ready <- file.path(input_dir, "Birth_women_final.csv") # путь к CSV-файлу с фактическими ASFR
levels_5y <- sprintf("%d-%d", seq(15, 45, 5), seq(19, 49, 5)) # ожидаемые интервалы 5-леток
need_cols <- c("age_band","women_parity_0","women_parity_1","women_parity_2",
"women_parity_3","women_parity_4","women_parity_5plus","women_total") # обязательные столбцы
if (file.exists(path_ready)) {
read_birth_data <- function(path) {
attempts <- list(
function() readr::read_csv(path, show_col_types = FALSE,
locale = readr::locale(grouping_mark = ",")),
function() readr::read_csv2(path, show_col_types = FALSE,
locale = readr::locale(decimal_mark = ",", grouping_mark = " ")),
function() readr::read_tsv(path, show_col_types = FALSE)
)
for (reader in attempts) {
result <- try(reader(), silent = TRUE)
if (!inherits(result, "try-error")) {
return(result)
}
}
stop(sprintf("Не удалось прочитать файл %s в поддерживаемых форматах CSV/TSV", path))
}
agg_raw <- read_birth_data(path_ready) # исходные данные рождаемости из файла
missing_cols <- setdiff(need_cols, names(agg_raw))
if (length(missing_cols) > 0) {
stop(sprintf(
"В файле %s отсутствуют необходимые колонки: %s",
path_ready,
paste(missing_cols, collapse = ", ")
))
}
agg <- agg_raw %>%
mutate(age_band = factor(age_band, levels = levels_5y)) %>%
arrange(age_band)
} else {
warning(sprintf(
"Файл %s не найден. Используем демонстрационную таблицу возрастной фертильности.",
path_ready
))
agg <- tribble(
~age_band, ~women_parity_0, ~women_parity_1, ~women_parity_2,
~women_parity_3, ~women_parity_4, ~women_parity_5plus, ~women_total,
"15-19", 950, 45, 5, 0, 0, 0, 1000,
"20-24", 520, 300, 130, 35, 10, 5, 1000,
"25-29", 180, 280, 340, 150, 35, 15, 1000,
"30-34", 90, 190, 320, 240, 110, 50, 1000,
"35-39", 70, 140, 260, 260, 170, 100, 1000,
"40-44", 60, 110, 210, 240, 210, 170, 1000,
"45-49", 55, 95, 180, 210, 210, 250, 1000
) %>%
mutate(age_band = factor(age_band, levels = levels_5y)) %>%
arrange(age_band)
}
if (any(agg$women_total <= 0)) {
stop("Колонка women_total содержит нулевые или отрицательные значения. Проверь источник данных.")
}
# Конструируем форму asfr_shape (сумма по 15–49 = 1)
idx1549 <- which(ages >= 15 & ages <= 49) # индексы фертильного возраста
asfr_shape <- rep(0, length(ages)) # заготовка под распределение по однолеткам
share0 <- agg$women_parity_0 / pmax(agg$women_total, 1) # доля бездетных женщин по возрастным группам
drop0 <- pmax(0, share0 - dplyr::lead(share0, default = 0)) # снижение бездетности ~ первые рождения
w5 <- if (sum(drop0) > 0) drop0 / sum(drop0) else rep(0, length(drop0)) # веса по 5-леткам
for (i in seq_along(w5)) {
lo <- 15 + (i-1)*5; hi <- lo + 4
idx <- which(ages >= lo & ages <= hi)
asfr_shape[idx] <- asfr_shape[idx] + w5[i] / length(idx)
}
if (sum(asfr_shape[idx1549]) > 0) asfr_shape[idx1549] <- asfr_shape[idx1549] / sum(asfr_shape[idx1549])
message(sprintf("Проверка: sum(asfr_shape[15:49]) = %.3f (должно быть 1.000)",
sum(asfr_shape[idx1549])))## Проверка: sum(asfr_shape[15:49]) = 1.000 (должно быть 1.000)# (2) УРОВЕНЬ РОЖДАЕМОСТИ И ИДЕНТИФИКАТОРЫ ГРУПП (как раньше)
# TFR_nat — сценарий MMM по годам 2021..2050 (вектор длины ny)
# Пример-заглушка: замените на ваш ряд MMM
# TFR_nat <- c(...) # заполняете из MMM
# Мультипликаторы уровней относительно TFR_nat (векторы длины ny)
# (начальные различия и их «схождение» управляются lambda)
# Фактические стартовые и целевые коэффициенты берём из сценария (могут быть переопределены данными)
m_M1_start <- scenario$m_M1_start
m_Mplus_start <- scenario$m_Mplus_start
m_RUS_start <- scenario$m_RUS_start
m_CORE_start <- scenario$m_CORE_start
# Целевые уровни на горизонте (к чему сходятся; по умолчанию 1,1,1)
m_M1_target <- scenario$m_M1_target
m_Mplus_target <- scenario$m_Mplus_target
m_RUS_target <- scenario$m_RUS_target
m_CORE_target <- scenario$m_CORE_target
# Векторы мультипликаторов из сценария
m_M1 <- scenario$m_M1
m_Mplus <- scenario$m_Mplus
m_RUS <- scenario$m_RUS
m_CORE <- scenario$m_CORE
# (3) В ЦИКЛЕ (Блок 8) остаётся как было:
# asfr_nat_t <- asfr_shape * TFR_nat[t]
# TFR_M1_t <- m_M1[t] * TFR_nat[t]; и т.д.Параметры смертности и улучшений
Здесь читается файл с возрастными коэффициентами смертности,
интерпретируются интервалы и формируются базовые вероятности смерти
q_base.
На выходе получаются векторы qF_base и
qM_base, а также функции для их корректного формирования
при отсутствии полных данных.
## ---------------------------
## 5) СМЕРТНОСТЬ: БАЗОВЫЕ q_x ПО ВОЗРАСТУ И ПОЛУ + УЛУЧШЕНИЕ
## ---------------------------
death_path <- file.path(input_dir, "death_age.csv")
parse_death_rates <- function(df, age_vector) {
max_age <- max(age_vector)
required_cols <- c("age", "female", "male")
if (!all(required_cols %in% names(df))) {
stop(
sprintf(
"Файл смертности должен содержать столбцы: %s",
paste(required_cols, collapse = ", ")
)
)
}
clean_numeric <- function(x) {
vals <- gsub("[[:space:]]", "", as.character(x))
vals <- gsub(",", ".", vals, fixed = FALSE)
suppressWarnings(as.numeric(vals))
}
qF <- rep(NA_real_, length(age_vector))
qM <- rep(NA_real_, length(age_vector))
for (i in seq_len(nrow(df))) {
span <- expand_age_span(df$age[i], max_age)
if (!length(span)) {
stop(sprintf("Не удалось интерпретировать возрастной интервал '%s'", df$age[i]))
}
idx <- which(age_vector %in% span)
if (length(idx) == 0) {
next
}
female_val <- clean_numeric(df$female[i])
male_val <- clean_numeric(df$male[i])
if (is.na(female_val) || is.na(male_val)) {
stop(sprintf("В строке %d обнаружены некорректные значения смертности", i))
}
qF[idx] <- female_val / 1000
qM[idx] <- male_val / 1000
}
if (any(is.na(qF)) || any(is.na(qM))) {
missing_ages <- age_vector[which(is.na(qF) | is.na(qM))]
stop(
sprintf(
"Для возрастов %s отсутствуют коэффициенты смертности. Проверьте файл death_age.csv.",
paste(unique(missing_ages), collapse = ", ")
)
)
}
list(female = pmin(qF, 0.99), male = pmin(qM, 0.99))
}
if (file.exists(death_path)) {
possible_delims <- c(",", ";", "\t", "|")
death_raw <- NULL
for (delim in possible_delims) {
death_raw <- tryCatch(
readr::read_delim(
death_path,
delim = delim,
trim_ws = TRUE,
show_col_types = FALSE
),
error = function(e) NULL
)
if (!is.null(death_raw)) break
}
if (is.null(death_raw)) {
stop(sprintf("Не удалось прочитать файл смертности %s", death_path))
}
names(death_raw) <- sanitize_names(names(death_raw))
death_rates <- parse_death_rates(death_raw, ages)
qF_base <- death_rates$female
qM_base <- death_rates$male
} else {
warning(
sprintf(
"Файл %s не найден. Используем сглаженные коэффициенты смертности как приближение.",
death_path
)
)
qF_base <- pmin(0.35, 0.00025 + (ages/100)^6 * 0.30)
qM_base <- pmin(0.40, 0.00035 + (ages/100)^6 * 0.32)
}
q_base <- cbind(qF_base, qM_base) # матрица q_x по возрасту и полуПрофиль распределения чистой миграции
Блок нормирует структуру миграции первого поколения и формирует
матрицу mig_prof, отражающую возрастно-половой состав
притока.
Если данные М1 доступны, используется фактическая пирамида; иначе генерируется сглаженное распределение, гарантируя корректную форму 101×2.
## ---------------------------
## 6) МИГРАЦИЯ: ПРОФИЛЬ РАСПРЕДЕЛЕНИЯ ЧИСТОГО ПРИТОКА ПО ВОЗРАСТ-ПОЛУ
## ---------------------------
## --- Профиль миграции: 101 возраста × 2 пола (сумма по всей матрице = 1)
# Стараемся повторить фактическую структуру М1 на 2021 год; при отсутствии данных используем сглаженное приближение
if (sum(M1_2021_as) > 0) {
mig_prof <- M1_2021_as / sum(M1_2021_as)
} else {
warning("Не удалось использовать базовую структуру М1 для профиля миграции — применяем сглаженное приближение")
stopifnot(abs(sum(M1_share_sex) - 1) < 1e-9)
mig_wF <- norm1(dnorm(ages, mean = 28, sd = 8))
mig_wM <- norm1(1.15 * dnorm(ages, mean = 28, sd = 8))
mig_prof <- cbind(mig_wF * M1_share_sex["F"], mig_wM * M1_share_sex["M"])
mig_prof <- mig_prof / sum(mig_prof)
}
## Нормируем матрицу (итоговая сумма = 1)
mig_prof <- mig_prof / sum(mig_prof)
## Защита: если вдруг mig_prof стал вектором/строкой — развернуть в 101x2
if (is.null(dim(mig_prof)) || !all(dim(mig_prof) == c(length(ages), 2))) {
mig_prof <- matrix(mig_prof, nrow = length(ages), ncol = 2, byrow = FALSE) # страховка: возвращаем форму 101×2
}
stopifnot(all(dim(mig_prof) == c(length(ages), 2)))Инициализация массивов динамики
Этот раздел создаёт трёхмерные хранилища для каждой демографической группы и заполняет первый временной срез исходными данными 2021 года.
Дополнительно формируются векторы former_M1,
Births_tot, Deaths_tot и заготовки долей
идентичности для будущего цикла.
## ---------------------------
## 7) ХРАНИЛИЩА ДИНАМИКИ
## ---------------------------
M1 <- array(0, dim = c(na, ns, ny)) # трёхмерный массив численности мигрантов первого поколения
Mpl <- array(0, dim = c(na, ns, ny)) # массив для потомков мигрантов (M+)
Rus_ar <- array(0, dim = c(na, ns, ny)) # массив для русских (коренное население)
Core_ar <- array(0, dim = c(na, ns, ny)) # массив для прочих коренных народов
RusUnd_ar <- array(0, dim = c(na, ns, ny)) # массив для «Россиян-неопределённых»
M1[,,1] <- M1_2021_as # базовое распределение мигрантов в 2021
Mpl[,,1] <- Mpl_2021_as # базовое распределение потомков мигрантов
Rus_ar[,,1] <- rus_2021_as # базовое распределение русских
Core_ar[,,1] <- core_2021_as # базовое распределение прочих коренных
RusUnd_ar[,,1]<- rus_unknown_2021_as # базовое распределение «Россиян-неопределённых»
former_M1 <- rep(0, ny) # накопленные натурализованные (агрегат)
Births_tot <- rep(0, ny); Deaths_tot <- rep(0, ny) # ежегодное число рождений/смертей
## Для отчёта идентичности (доли в общем населении)
id_RUS_share <- id_IND_share <- id_MIGID_share <- rep(NA_real_, ny) # временные ряды долей идентичностиПроверка размерностей перед моделированием
Блок сверяет формы всех ключевых массивов и таблиц смертности, чтобы убедиться в согласованности размерностей перед основным циклом.
Он не изменяет данные, но вызывает check_tensor_shape
для q_base, M1, Mpl,
Rus_ar, Core_ar и RusUnd_ar,
предотвращая ошибки индексации.
## === ВЫРАВНИВАНИЕ РАЗМЕРОВ ПЕРЕД БЛОКОМ 8 ===
## Получаем фактические размеры из массива M1 (возраст × пол × год)
na_chk <- dim(M1)[1]
ns_chk <- dim(M1)[2]
ny_chk <- dim(M1)[3]
stopifnot(ns_chk == ns) # два пола
stopifnot(na_chk == na)
stopifnot(ny_chk == ny)
check_tensor_shape(q_base, c(na, ns), "q_base") # валидация матриц смертности
check_tensor_shape(M1, c(na, ns, ny), "M1") # проверка формы массива M1
check_tensor_shape(Mpl, c(na, ns, ny), "Mpl") # проверка формы массива M+
check_tensor_shape(Rus_ar, c(na, ns, ny), "Rus_ar") # проверка формы массива русских
check_tensor_shape(Core_ar, c(na, ns, ny), "Core_ar") # проверка формы массива других коренных
check_tensor_shape(RusUnd_ar, c(na, ns, ny), "RusUnd_ar") # проверка формы массива «Россиян-неопределённых»Основной годовой цикл модели
Раздел описывает итерацию по годам, где рассчитываются смертность, миграция, рождаемость и обновления идентичностей на каждом шаге горизонта.
На выходе обновляются массивы состояний, агрегаты рождений и смертей, накопленные натурализации и временные ряды долей идентичности.
## ---------------------------
## 8) ЦИКЛ ПО ГОДАМ (t = 2021..2050)
## ---------------------------
is_fert_age <- ages >= 15 & ages <= 49 # логический вектор фертильного возраста женщин
srb <- c(m = 0.512, f = 1 - 0.512) # отношение полов при рождении (share of male/female)
for (t in 1:ny) {
## --- ТЕКУЩИЕ СТАНЫ И ПРОВЕРКИ ФОРМ ---
states_t <- list(
M1 = M1[,,t],
Mpl = Mpl[,,t],
RUS = Rus_ar[,,t],
IND = Core_ar[,,t],
RUS_UND = RusUnd_ar[,,t]
) # срез массивов на начало года t
for (state_name in names(states_t)) {
check_tensor_shape(states_t[[state_name]], c(na, ns), sprintf("%s_t", state_name))
}
Pop_t <- states_t$M1 + states_t$Mpl + states_t$RUS + states_t$IND # население на начало года (для смертности)
## --- СМЕРТНОСТЬ И СТАРЕНИЕ ---
q_t <- q_base * scenario$mort_impr[t] # [na × ns] годовые вероятности смерти
Sur <- make_survival_matrix(q_t) # вероятность выживания по возрасту и полу
states_survived <- advance_with_mortality(states_t, Sur) # численность после старения/смертности
## --- МЕЖДУНАРОДНАЯ МИГРАЦИЯ (вся → M1) ---
states_survived$M1 <- apply_migration(states_survived$M1, scenario$net_mig[t], mig_prof) # добавляем чистый миграционный поток
## --- РОЖДЕНИЯ И НОВОРОЖДЁННЫЕ ---
birth_info <- compute_births( # расчёт рождений и добавление новорождённых
states = states_survived,
asfr_shape = asfr_shape,
TFR_nat_t = scenario$TFR_nat[t],
multipliers = c(
M1 = scenario$m_M1[t],
Mpl = scenario$m_Mplus[t],
RUS = scenario$m_RUS[t],
IND = scenario$m_CORE[t],
RUS_UND = scenario$m_RUS[t]
),
rho_t = scenario$rho[t],
fert_age_mask = is_fert_age,
srb = srb,
idxF = idxF,
idxM = idxM
)
## --- НАТУРАЛИЗАЦИЯ (бывшие мигранты) ---
former_M1[t] <- (ifelse(t == 1, 0, former_M1[t-1])) + scenario$kappa_nat[t] * sum(states_t$M1) # накопление натурализованных
## --- ЗАПИСЬ СОСТОЯНИЯ НА КОНЕЦ ГОДА t ---
M1[,,t] <- birth_info$states$M1
Mpl[,,t] <- birth_info$states$Mpl
Rus_ar[,,t] <- birth_info$states$RUS
Core_ar[,,t] <- birth_info$states$IND
RusUnd_ar[,,t] <- birth_info$states$RUS_UND
## --- АГРЕГАТЫ ДЛЯ ОТЧЁТА ---
Births_tot[t] <- birth_info$births$total # итоговое число рождений в году t
Deaths_tot[t] <- sum(Pop_t * q_t) # смерти в году t (по началу года)
identity_shares <- update_identity_shares(
birth_info$states,
alpha = scenario$alpha_assim[t],
beta = scenario$beta_M1_core[t],
fallback_share_rus_core = share_RUS_2021
)
id_RUS_share[t] <- identity_shares[["RUS"]]
id_IND_share[t] <- identity_shares[["IND"]]
id_MIGID_share[t] <- identity_shares[["MIGID"]]
## --- ПОДГОТОВКА t+1 (копия нужна, если вне цикла будут модификации до следующей итерации)
if (t < ny) {
M1[,,t+1] <- M1[,,t]
Mpl[,,t+1] <- Mpl[,,t]
Rus_ar[,,t+1] <- Rus_ar[,,t]
Core_ar[,,t+1] <- Core_ar[,,t]
RusUnd_ar[,,t+1] <- RusUnd_ar[,,t]
}
}Постобработка и агрегирование показателей
Блок суммирует трёхмерные массивы в годовые ряды, формирует таблицу
out с ключевыми демографическими индикаторами и готовит
округлённое представление для вывода.
Результаты включают показатели численности, естественного и миграционного прироста, долей идентичностей и траекторий TFR.
## ---------------------------
## 9) АГРЕГИРОВАННЫЙ ВЫХОД, КАК РАНЬШЕ (без colSums)
## ---------------------------
## Суммы по каждому году (возраст × пол суммируем осью 3 -> вектор длины ny)
pop_M1 <- apply(M1, 3, sum) # численность M1 по годам
pop_Mpl <- apply(Mpl, 3, sum) # численность M+ по годам
pop_RUS_total <- apply(Rus_ar, 3, sum) # численность русских (включая неопределённую национальность)
pop_RUS_und <- apply(RusUnd_ar, 3, sum) # диагностический массив «Россиян-неопределённых» (ожидается ноль)
pop_CORE <- apply(Core_ar, 3, sum) # численность прочих коренных по годам
## Общая численность по годам — просто векторная сумма
pop_total <- pop_M1 + pop_Mpl + pop_RUS_total + pop_CORE # совокупное население
## Дополнительные агрегаты для анализа
natural_increase <- Births_tot - Deaths_tot # естественный прирост
cumulative_migration <- cumsum(scenario$net_mig) # накопленный миграционный приток
share_M1 <- ifelse(pop_total > 0, pop_M1 / pop_total, NA_real_) # доля M1 в населении
share_Mpl <- ifelse(pop_total > 0, pop_Mpl / pop_total, NA_real_) # доля M+
share_RUS_total <- ifelse(pop_total > 0, pop_RUS_total / pop_total, NA_real_) # доля русских в общем населении
share_RUS_und <- ifelse((pop_total + pop_RUS_und) > 0, pop_RUS_und / (pop_total + pop_RUS_und), NA_real_) # доля «Россиян-неопределённых» (диагностика)
share_CORE_b <- ifelse(pop_total > 0, pop_CORE / pop_total, NA_real_) # доля прочих коренных
# Основной агрегированный результат (без округления) для аналитики и экспорта
out <- tibble(
year = years,
pop_total = pop_total,
M1 = pop_M1,
Mplus = pop_Mpl,
RUS_core_total = pop_RUS_total,
RUS_core_undetermined = pop_RUS_und,
IND_core = pop_CORE,
former_M1 = former_M1,
births_total = Births_tot,
deaths_total = Deaths_tot,
natural_increase = natural_increase,
net_migration = scenario$net_mig,
cumulative_migration = cumulative_migration,
share_M1 = share_M1,
share_Mplus = share_Mpl,
share_RUS_total = share_RUS_total,
share_RUS_undetermined = share_RUS_und,
share_IND = share_CORE_b,
TFR_nat = scenario$TFR_nat,
TFR_M1 = scenario$m_M1 * scenario$TFR_nat,
TFR_Mplus = scenario$m_Mplus * scenario$TFR_nat,
TFR_RUS = scenario$m_RUS * scenario$TFR_nat,
TFR_IND = scenario$m_CORE * scenario$TFR_nat,
id_RUS_share = id_RUS_share,
id_IND_share = id_IND_share,
id_MIGID_share = id_MIGID_share
)
# Представление с округлением для вывода в консоль
out_print <- out %>%
mutate(
across(c(pop_total, M1, Mplus, RUS_core_total, RUS_core_undetermined, IND_core, former_M1, births_total, deaths_total,
natural_increase, net_migration, cumulative_migration), round),
across(starts_with("TFR"), ~round(.x, 3)),
across(starts_with("share_"), ~round(.x, 3)),
across(starts_with("id_"), ~round(.x, 3))
)
print(head(out_print, 5))## # A tibble: 5 × 26
## year pop_total M1 Mplus RUS_core_total RUS_core_undetermined IND_core
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2021 146485234 7032889 1383105 124076087 0 13993152
## 2 2022 145390370 7175402 1705476 122558884 0 13950608
## 3 2023 144340873 7323470 2014555 121091418 0 13911430
## 4 2024 143309043 7476874 2310540 119648920 0 13872709
## 5 2025 142289941 7635443 2596897 118221405 0 13836197
## # ℹ 19 more variables: former_M1 <dbl>, births_total <dbl>, deaths_total <dbl>,
## # natural_increase <dbl>, net_migration <dbl>, cumulative_migration <dbl>,
## # share_M1 <dbl>, share_Mplus <dbl>, share_RUS_total <dbl>,
## # share_RUS_undetermined <dbl>, share_IND <dbl>, TFR_nat <dbl>, TFR_M1 <dbl>,
## # TFR_Mplus <dbl>, TFR_RUS <dbl>, TFR_IND <dbl>, id_RUS_share <dbl>,
## # id_IND_share <dbl>, id_MIGID_share <dbl>print(tail(out_print, 5))## # A tibble: 5 × 26
## year pop_total M1 Mplus RUS_core_total RUS_core_undetermined IND_core
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2046 122887696 11436442 9.40e6 89803297 0 12244777
## 2 2047 121926383 11584889 9.80e6 88435232 0 12110754
## 3 2048 120982786 11726708 1.02e7 87087771 0 11975938
## 4 2049 120060698 11861600 1.06e7 85762227 0 11840607
## 5 2050 119149258 11989079 1.10e7 84444214 0 11704547
## # ℹ 19 more variables: former_M1 <dbl>, births_total <dbl>, deaths_total <dbl>,
## # natural_increase <dbl>, net_migration <dbl>, cumulative_migration <dbl>,
## # share_M1 <dbl>, share_Mplus <dbl>, share_RUS_total <dbl>,
## # share_RUS_undetermined <dbl>, share_IND <dbl>, TFR_nat <dbl>, TFR_M1 <dbl>,
## # TFR_Mplus <dbl>, TFR_RUS <dbl>, TFR_IND <dbl>, id_RUS_share <dbl>,
## # id_IND_share <dbl>, id_MIGID_share <dbl>Визуализация динамики идентичностей
Этот раздел преобразует таблицу out в формат для
построения столбчатой диаграммы и визуализирует изменение
самоидентификации населения.
На выходе формируется объект графика ggplot p, который
выводится в документ и сохраняется в файл
identity_stacked_2021_2050.png.
## ==== Стековая диаграмма «численность по идентичности, ежегодно» ====
# 1) Подготовка данных (берём из data.frame `out`)
horizon_years <- length(unique(out$year)) # продолжительность горизонта (для динамической графики)
LABEL_EVERY <- max(1, floor(horizon_years / 15)) # шаг подписей по оси X
PERC_MIN <- if (horizon_years > 25) 0.04 else 0.03 # минимальная доля для вывода подписи
stack_df <- out %>%
transmute(
year,
pop_total,
RUS = pop_total * id_RUS_share,
IND = pop_total * id_IND_share,
MIGID = pop_total * id_MIGID_share
) %>%
pivot_longer(c(RUS, IND, MIGID), names_to = "group", values_to = "pop") %>%
mutate(group = factor(group,
levels = c("RUS","IND","MIGID"),
labels = c("Русские", "Другие коренные", "Мигрантская идентичность"))) %>%
group_by(year) %>%
mutate(share = pop / sum(pop),
label = if_else(share >= PERC_MIN, percent(share, accuracy = 1), "")) %>%
ungroup()
year_breaks <- as.character(seq(min(stack_df$year), max(stack_df$year), by = LABEL_EVERY)) # подписи по оси X
# 2) Построение графика
p <- ggplot(stack_df, aes(x = factor(year), y = pop, fill = group)) + # основной столбчатый график идентичности
geom_col(width = 0.9) +
geom_text(aes(label = label),
position = position_stack(vjust = 0.5),
color = "white", size = 3.1, fontface = "bold") +
scale_y_continuous(labels = label_number(scale = 1e-6, accuracy = 0.1, suffix = " млн"),
expand = expansion(mult = c(0, 0.04))) +
scale_x_discrete(breaks = year_breaks) +
scale_fill_manual(values = c("Русские" = "#4C78A8",
"Другие коренные" = "#59A14F",
"Мигрантская идентичность" = "#F28E2B")) +
labs(
x = NULL, y = "Численность, млн человек", fill = "Идентичность",
title = "Структура населения по самоидентификации (ежегодно)",
subtitle = "Высота столбика = общая численность; заливка и подписи = доли идентичностей"
) +
theme_minimal(base_size = 12) +
theme(
panel.grid.minor = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1, size = 8),
legend.position = "top",
legend.box = "horizontal",
plot.margin = margin(10, 15, 10, 10)
)
print(p)
# 3) Сохранить картинку покрупнее (чтобы ничего не наезжало)
ggsave("identity_stacked_2021_2050.png", p, width = 18, height = 9, dpi = 200) # экспорт рисунка в PNGЗаготовка для дополнительных экспериментов
Заключительный блок оставлен пустым для будущих расширений сценария и может использоваться для временных расчётов или отладочных команд.
Сейчас он не выполняет операций, но служит понятным маркером конца документа и места для добавления новых шагов.