1. Выбор датасета и разведочный анализ (EDA)

Используем встроенный датасет mtcars — характеристики 32 автомобилей (журнал Motor Trend, 1974 г.): расход топлива mpg, масса wt, тип коробки передач am (0 — автомат, 1 — ручная) и др.

Почему именно он. Датасет маленький (32 наблюдения) — а это как раз тот случай, когда классические формулы для доверительных интервалов работают плохо, и бутстрап особенно полезен. Кроме того, в нём есть и естественная задача «сравнить две группы» (для Python-части), и естественная регрессия mpg ~ wt (для R-части).

library(ggplot2)
library(dplyr)
library(boot)
library(gridExtra)
library(reticulate)

cars <- mtcars
cars$transmission <- factor(cars$am, levels = c(0, 1), labels = c("Автомат", "Ручная"))
head(cars)
summary(cars[, c("mpg", "wt", "hp", "am")])
##       mpg              wt              hp              am        
##  Min.   :10.40   Min.   :1.513   Min.   : 52.0   Min.   :0.0000  
##  1st Qu.:15.43   1st Qu.:2.581   1st Qu.: 96.5   1st Qu.:0.0000  
##  Median :19.20   Median :3.325   Median :123.0   Median :0.0000  
##  Mean   :20.09   Mean   :3.217   Mean   :146.7   Mean   :0.4062  
##  3rd Qu.:22.80   3rd Qu.:3.610   3rd Qu.:180.0   3rd Qu.:1.0000  
##  Max.   :33.90   Max.   :5.424   Max.   :335.0   Max.   :1.0000
cat("Размер датасета:", nrow(cars), "строк,", ncol(cars), "столбцов\n")
## Размер датасета: 32 строк, 12 столбцов
cat("Пропущенные значения:", sum(is.na(cars)), "\n")
## Пропущенные значения: 0
table(cars$transmission)
## 
## Автомат  Ручная 
##      19      13
cars %>%
  group_by(transmission) %>%
  summarise(n = n(),
            mean_mpg = round(mean(mpg), 2),
            median_mpg = median(mpg),
            sd_mpg = round(sd(mpg), 2))

Визуализация данных

p1 <- ggplot(cars, aes(x = transmission, y = mpg, fill = transmission)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.1, alpha = 0.5) +
  labs(title = "Расход топлива по типу коробки",
       x = "Коробка передач", y = "mpg (миль на галлон)") +
  theme_minimal() + theme(legend.position = "none")

p2 <- ggplot(cars, aes(x = mpg, fill = transmission)) +
  geom_histogram(bins = 12, alpha = 0.6, position = "identity") +
  labs(title = "Гистограмма mpg", x = "mpg", y = "Частота", fill = "Коробка") +
  theme_minimal()

grid.arrange(p1, p2, ncol = 2)

ggplot(cars, aes(x = wt, y = mpg, color = transmission)) +
  geom_point(size = 2.5, alpha = 0.85) +
  labs(title = "Связь массы автомобиля и расхода топлива",
       x = "Масса (1000 фунтов)", y = "mpg", color = "Коробка") +
  theme_minimal()

r_wt <- cor(cars$mpg, cars$wt)

Вывод по EDA. Датасет содержит 32 наблюдения без пропусков. Автомобили с ручной коробкой в среднем экономичнее (средний mpg 24.4 против 17.1 у автоматов), но группы сильно перекрываются. Между массой и расходом топлива наблюдается сильная отрицательная линейная связь (r = -0.868): чем тяжелее машина, тем меньше миль на галлон. Стоит помнить, что ручные коробки в выборке чаще стоят на более лёгких машинах, поэтому различие по am частично объясняется массой (это уже вопрос причинности, который выходит за рамки бутстрапа).

2. Бутстрап на языке Python

Используем пакет reticulate, чтобы выполнять Python-код прямо в R Markdown. Данные передаются из R в Python через объект r.cars.

# Если Python не найден — укажите путь вручную, например:
# use_python("C:/Users/ИмяПользователя/AppData/Local/Programs/Python/Python311/python.exe")
# и установите пакеты: py_install(c("numpy", "pandas", "matplotlib"))

Задача. Сравнить медианный расход топлива у автомобилей с ручной и автоматической коробкой: статистика — разность медиан (ручная − автомат). Алгоритм — из слайда 19: бутстрап-выборки с возвращением для каждой группы, разность квантилей, перцентильный доверительный интервал. Чтобы получить 90-й перцентиль, как на слайде, достаточно поменять QUANTILE = 0.5 на 0.9.

Отладка исходного кода. В слайде размер бутстрап-выборки задан произвольно (n = 1000). Это неверно: бутстрап-выборка должна иметь тот же размер, что и исходная (здесь 19 и 13). Иначе интервал получается искусственно узким, а гистограмма — вырожденной. Поэтому размеры взяты как len(values_a) и len(values_b).

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

np.random.seed(42)

# --- Данные из R ---
df = r.cars[["mpg", "wt", "am"]]
values_a = df.loc[df["am"] == 0, "mpg"].values   # автомат
values_b = df.loc[df["am"] == 1, "mpg"].values   # ручная

print(f"Автомат: {len(values_a)} наблюдений, среднее = {values_a.mean():.2f}")
## Автомат: 19 наблюдений, среднее = 17.15
print(f"Ручная:  {len(values_b)} наблюдений, среднее = {values_b.mean():.2f}")
## Ручная:  13 наблюдений, среднее = 24.39
# --- Функция построения ДИ ---
def get_percentile_ci(bootstrap_stats, pe, alpha):
    left, right = np.quantile(bootstrap_stats, [alpha / 2, 1 - alpha / 2])
    return left, right

# --- Параметры ---
QUANTILE = 0.5     # 0.5 — медиана, 0.9 — 90-й перцентиль
B = 10000          # число итераций
alpha = 0.05       # 95% ДИ
n_a, n_b = len(values_a), len(values_b)

# --- Точечная оценка ---
pe = np.quantile(values_b, QUANTILE) - np.quantile(values_a, QUANTILE)
print(f"\nТочечная оценка разности (ручная - автомат): {pe:.3f}")
## 
## Точечная оценка разности (ручная - автомат): 5.500
# --- Бутстрап ---
bootstrap_metrics_a = np.quantile(np.random.choice(values_a, (B, n_a), True), QUANTILE, axis=1)
bootstrap_metrics_b = np.quantile(np.random.choice(values_b, (B, n_b), True), QUANTILE, axis=1)
bootstrap_stats = bootstrap_metrics_b - bootstrap_metrics_a

# --- ДИ и вывод ---
ci = get_percentile_ci(bootstrap_stats, pe, alpha)
has_effect = not (ci[0] < 0 < ci[1])
se = bootstrap_stats.std(ddof=1)

print(f"Стандартная ошибка (бутстрап): {se:.3f}")
## Стандартная ошибка (бутстрап): 3.205
print(f"{(1 - alpha) * 100:.0f}% доверительный интервал: ({ci[0]:.2f}, {ci[1]:.2f})")
## 95% доверительный интервал: (2.20, 14.00)
print(f"Отличия статистически значимые: {has_effect}")
## Отличия статистически значимые: True
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
fig.suptitle("Бутстрап: разность медиан mpg (ручная − автомат), mtcars", fontsize=14)

# 1. Гистограмма бутстрап-распределения
axes[0].hist(bootstrap_stats, bins=40, color="steelblue", alpha=0.75, edgecolor="white")
axes[0].axvline(pe, color="red", linewidth=2, label=f"Оценка: {pe:.2f}")
axes[0].axvline(ci[0], color="orange", linestyle="--", linewidth=1.5,
                label=f"ДИ: [{ci[0]:.2f}; {ci[1]:.2f}]")
axes[0].axvline(ci[1], color="orange", linestyle="--", linewidth=1.5)
axes[0].axvline(0, color="black", linestyle=":", linewidth=1.5, label="Ноль")
axes[0].set_title("Бутстрап-распределение\nразности медиан")
axes[0].set_xlabel("Разность медиан mpg")
axes[0].set_ylabel("Частота")
axes[0].legend(fontsize=8)

# 2. Боксплоты + точки
_ = axes[1].boxplot([values_a, values_b], patch_artist=True,
                    boxprops=dict(facecolor="lightblue", alpha=0.7),
                    medianprops=dict(color="darkred", linewidth=2))
axes[1].set_xticks([1, 2])
## [<matplotlib.axis.XTick object at 0x000001FD7A44AAD0>, <matplotlib.axis.XTick object at 0x000001FD7A44A210>]
axes[1].set_xticklabels(["Автомат", "Ручная"])
## [Text(1, 0, 'Автомат'), Text(2, 0, 'Ручная')]
for i, vals in enumerate([values_a, values_b], start=1):
    axes[1].scatter(np.random.normal(i, 0.05, len(vals)), vals,
                    color="black", alpha=0.5, s=18, zorder=3)
## <matplotlib.collections.PathCollection object at 0x000001FD7A5486E0>
## <matplotlib.collections.PathCollection object at 0x000001FD7A4DF610>
axes[1].axhline(np.quantile(values_a, QUANTILE), color="blue", linestyle="--", alpha=0.5,
                label="медиана (автомат)")
## <matplotlib.lines.Line2D object at 0x000001FD7A4770E0>
axes[1].axhline(np.quantile(values_b, QUANTILE), color="red", linestyle="--", alpha=0.5,
                label="медиана (ручная)")
## <matplotlib.lines.Line2D object at 0x000001FD7A548AD0>
axes[1].set_title("Расход топлива по типу коробки")
## Text(0.5, 1.0, 'Расход топлива по типу коробки')
axes[1].set_ylabel("mpg (миль на галлон)")
## Text(0, 0.5, 'mpg (миль на галлон)')
axes[1].legend(fontsize=8)
## <matplotlib.legend.Legend object at 0x000001FD7A548C20>
# 3. Эмпирическая функция распределения
sorted_stats = np.sort(bootstrap_stats)
cdf = np.arange(1, len(sorted_stats) + 1) / len(sorted_stats)
axes[2].plot(sorted_stats, cdf, color="purple", linewidth=2)
## [<matplotlib.lines.Line2D object at 0x000001FD7A5492B0>]
axes[2].axvline(ci[0], color="orange", linestyle="--", label=f"2.5% = {ci[0]:.2f}")
## <matplotlib.lines.Line2D object at 0x000001FD7A549010>
axes[2].axvline(ci[1], color="orange", linestyle="--", label=f"97.5% = {ci[1]:.2f}")
## <matplotlib.lines.Line2D object at 0x000001FD7A549160>
axes[2].axvline(0, color="black", linestyle=":", label="Ноль")
## <matplotlib.lines.Line2D object at 0x000001FD7A549400>
axes[2].set_title("Кумулятивное распределение\nбутстрап-статистики")
## Text(0.5, 1.0, 'Кумулятивное распределение\nбутстрап-статистики')
axes[2].set_xlabel("Разность медиан")
## Text(0.5, 0, 'Разность медиан')
axes[2].set_ylabel("Вероятность")
## Text(0, 0.5, 'Вероятность')
axes[2].legend(fontsize=8)
## <matplotlib.legend.Legend object at 0x000001FD7A549550>
plt.tight_layout()
plt.savefig("bootstrap_python.png", dpi=150, bbox_inches="tight")
plt.show()

Вывод по Python-бутстрапу. Медианный расход топлива у автомобилей с ручной коробкой выше на 5.5 mpg. 95% перцентильный доверительный интервал: (2.2; 14). Ноль не входит в интервал, поэтому различие статистически значимо. Интервал довольно широкий (бутстрап-стандартная ошибка ≈ 3.21) — это ожидаемо при 13 и 19 наблюдениях в группах, а бутстрап-распределение медианы неровное (дискретное), потому что медиана берётся из небольшого числа реальных значений.

3. Бутстрап на языке R

Алгоритм из слайда 21: пакет boot. Задача: оценить устойчивость коэффициента детерминации R² модели mpg ~ wt и построить его доверительный интервал.

# Установка пакетов (один раз)
install.packages(c("boot", "ggplot2", "gridExtra", "dplyr"))
# Функция-статистика: принимает данные и индексы бутстрап-выборки, возвращает R²
rsq_function <- function(data, indices) {
  d   <- data[indices, ]
  fit <- lm(mpg ~ wt, data = d)
  summary(fit)$r.squared
}

set.seed(42)
reps <- boot(data = cars, statistic = rsq_function, R = 2000)
print(reps)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = cars, statistic = rsq_function, R = 2000)
## 
## 
## Bootstrap Statistics :
##      original      bias    std. error
## t1* 0.7528328 0.004827619  0.05883449
# Доверительные интервалы: перцентильный и BCa
ci_result <- boot.ci(reps, type = c("perc", "bca"))
print(ci_result)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = reps, type = c("perc", "bca"))
## 
## Intervals : 
## Level     Percentile            BCa          
## 95%   ( 0.6277,  0.8570 )   ( 0.5843,  0.8392 )  
## Calculations and Intervals on Original Scale
ci_low  <- ci_result$bca[4]
ci_high <- ci_result$bca[5]

cat("\nОригинальное R²:", round(reps$t0, 4), "\n")
## 
## Оригинальное R²: 0.7528
cat("95% ДИ (BCa):", round(ci_low, 4), "-", round(ci_high, 4), "\n")
## 95% ДИ (BCa): 0.5843 - 0.8392
boot_df <- data.frame(r_squared = reps$t)

# 1. Гистограмма бутстрап-распределения R²
g1 <- ggplot(boot_df, aes(x = r_squared)) +
  geom_histogram(bins = 40, fill = "steelblue", alpha = 0.7, color = "white") +
  geom_vline(xintercept = reps$t0, color = "red", linewidth = 1.2) +
  geom_vline(xintercept = c(ci_low, ci_high), color = "orange",
             linewidth = 1, linetype = "dashed") +
  labs(title = "Бутстрап-распределение R²",
       subtitle = paste0("95% ДИ (BCa): [", round(ci_low, 3), "; ", round(ci_high, 3), "]"),
       x = "R²", y = "Частота") +
  theme_minimal()

# 2. Q-Q график
g2 <- ggplot(boot_df, aes(sample = r_squared)) +
  stat_qq(color = "steelblue", alpha = 0.5) +
  stat_qq_line(color = "red", linewidth = 1) +
  labs(title = "Q-Q график", subtitle = "Нормальность бутстрап-распределения",
       x = "Теоретические квантили", y = "Выборочные квантили") +
  theme_minimal()

# 3. Данные и линия регрессии
g3 <- ggplot(cars, aes(x = wt, y = mpg)) +
  geom_point(aes(color = transmission), alpha = 0.8, size = 2.5) +
  geom_smooth(method = "lm", se = TRUE, color = "black", linewidth = 1) +
  labs(title = "Регрессия: mpg ~ wt",
       subtitle = paste("R² =", round(reps$t0, 3)),
       x = "Масса (1000 фунтов)", y = "mpg", color = "Коробка") +
  theme_minimal()

grid.arrange(g1, g2, g3, ncol = 3)

Вывод по R-бутстрапу. Коэффициент детерминации R² ≈ 0.753: масса автомобиля объясняет около 75% вариации расхода топлива. 95% BCa-интервал: [0.584; 0.839], бутстрап-стандартная ошибка ≈ 0.059. Интервал заметно шире, чем для iris (там R² = 0,93 при 150 наблюдениях): малый объём выборки даёт большую неопределённость. На Q-Q графике видно, что бутстрап-распределение R² несимметрично (левый «хвост» тяжелее, так как R² ограничен единицей сверху), поэтому BCa-интервал, корректирующий смещение и асимметрию, предпочтительнее обычного нормального интервала.

Общие выводы

  1. Бутстрап — непараметрический метод оценки точности статистик через многократную выборку с возвращением из исходных данных; он не требует предположений о распределении и подходит для «неудобных» статистик (медиана, перцентили, R²), для которых нет простых формул стандартной ошибки.
  2. Python: разность медианного mpg между ручной и автоматической коробкой равна 5.5, 95% ДИ = (2.2; 14) — эффект значим, но интервал широкий из-за малых групп.
  3. R: R² модели mpg ~ wt равен 0.753 с BCa-интервалом [0.584; 0.839] — связь сильная, но на 32 наблюдениях оценка менее точна, чем на больших выборках.
  4. Важное практическое замечание: размер бутстрап-выборки должен совпадать с размером исходной выборки. Именно эту ошибку пришлось исправить при адаптации кода со слайда.
  5. Оба языка дают одинаковую логику метода: Python удобен векторизованной реализацией «вручную» (numpy), R — готовым пакетом boot с встроенными методами построения ДИ (percentile, BCa, basic, normal).