Загрузка библиотек

library(ggplot2)
library(dplyr)

Подготовка данных

data(iris)

versicolor <- iris$Petal.Length[iris$Species == "versicolor"]
virginica <- iris$Petal.Length[iris$Species == "virginica"]

get_percentile_ci <- function(bootstrap_stats, pe, alpha) {
  left <- quantile(bootstrap_stats, alpha / 2)
  right <- quantile(bootstrap_stats, 1 - alpha / 2)
  return(c(left, right))
}

n_versicolor <- length(versicolor)
n_virginica <- length(virginica)
B <- 10000
alpha <- 0.05

pe <- median(virginica) - median(versicolor)

Бутстрап

set.seed(123)
bootstrap_versicolor <- matrix(sample(versicolor, B * n_versicolor, replace = TRUE), nrow = B)
bootstrap_virginica <- matrix(sample(virginica, B * n_virginica, replace = TRUE), nrow = B)

bootstrap_median_versicolor <- apply(bootstrap_versicolor, 1, median)
bootstrap_median_virginica <- apply(bootstrap_virginica, 1, median)
bootstrap_stats <- bootstrap_median_virginica - bootstrap_median_versicolor

ci <- get_percentile_ci(bootstrap_stats, pe, alpha)
has_effect <- !(ci[1] < 0 && ci[2] > 0)

Результаты

print("==================================================")
## [1] "=================================================="
print("РЕЗУЛЬТАТЫ БУТСТРАПА (R)")
## [1] "РЕЗУЛЬТАТЫ БУТСТРАПА (R)"
print("==================================================")
## [1] "=================================================="
cat(sprintf("Медиана длины лепестка Versicolor: %.2f см\n", median(versicolor)))
## Медиана длины лепестка Versicolor: 4.35 см
cat(sprintf("Медиана длины лепестка Virginica: %.2f см\n", median(virginica)))
## Медиана длины лепестка Virginica: 5.55 см
cat(sprintf("Изменение медианы (Virginica - Versicolor): %.2f см\n", pe))
## Изменение медианы (Virginica - Versicolor): 1.20 см
cat(sprintf("95%% доверительный интервал: [%.2f, %.2f]\n", ci[1], ci[2]))
## 95% доверительный интервал: [0.90, 1.45]
cat(sprintf("Разница статистически значима: %s\n", has_effect))
## Разница статистически значима: TRUE

График 1: Распределения медиан

par(mfrow = c(1, 2))

hist(bootstrap_median_versicolor, breaks = 50, col = "blue", border = "black",
     main = "Бутстрап распределение медианы Versicolor",
     xlab = "Медиана длины лепестка (см)", ylab = "Частота")
abline(v = median(versicolor), col = "red", lty = 2, lwd = 2)
legend("topright", legend = "Исходная медиана", col = "red", lty = 2, lwd = 2)

hist(bootstrap_median_virginica, breaks = 50, col = "green", border = "black",
     main = "Бутстрап распределение медианы Virginica",
     xlab = "Медиана длины лепестка (см)", ylab = "Частота")
abline(v = median(virginica), col = "red", lty = 2, lwd = 2)
legend("topright", legend = "Исходная медиана", col = "red", lty = 2, lwd = 2)

par(mfrow = c(1, 1))

График 2: Разница медиан

hist(bootstrap_stats, breaks = 50, col = "purple", border = "black",
     main = "Бутстрап распределение разницы медиан (Virginica - Versicolor)",
     xlab = "Разница медиан длины лепестка (см)", ylab = "Частота")
abline(v = pe, col = "red", lty = 2, lwd = 2)
abline(v = ci[1], col = "blue", lty = 3, lwd = 2)
abline(v = ci[2], col = "blue", lty = 3, lwd = 2)
abline(v = 0, col = "black", lty = 1, lwd = 1)

legend("topright",
       legend = c("Исходная разница медиан",
                  sprintf("Нижняя граница 95%% ДИ: %.2f", ci[1]),
                  sprintf("Верхняя граница 95%% ДИ: %.2f", ci[2])),
       col = c("red", "blue", "blue"),
       lty = c(2, 3, 3),
       lwd = 2)

Выводы