Практическая работа №6: Бутстрап на R и Python
Ярослав Гринко
2026-05-30
1. Выбор датасета и разведочный анализ (EDA)
Мы используем встроенный датасет iris —
это данные об ирисах: 150 цветков, 3 вида, 4 числовых признака
(длина/ширина лепестков и чашелистиков).
# Загружаем нужные библиотеки
library(ggplot2)
library(dplyr)
library(boot)
# Смотрим на первые строки датасета
head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa# Краткая статистика по всем столбцам
summary(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
## # Размер датасета: сколько строк и столбцов
cat("Размер датасета:", nrow(iris), "строк,", ncol(iris), "столбцов\n")## Размер датасета: 150 строк, 5 столбцов# Проверяем, есть ли пропуски
cat("Пропущенные значения:", sum(is.na(iris)), "\n")## Пропущенные значения: 0Визуализация данных
# График 1: Распределение длины лепестка по видам
p1 <- ggplot(iris, aes(x = Species, y = Petal.Length, fill = Species)) +
geom_boxplot(alpha = 0.7) +
labs(title = "Длина лепестка по видам ириса",
x = "Вид", y = "Длина лепестка (см)") +
theme_minimal() +
theme(legend.position = "none")
# График 2: Гистограмма длины лепестка
p2 <- ggplot(iris, aes(x = Petal.Length, fill = Species)) +
geom_histogram(bins = 20, alpha = 0.6, position = "identity") +
labs(title = "Гистограмма длины лепестка",
x = "Длина лепестка (см)", y = "Частота") +
theme_minimal()
# Выводим оба графика рядом
gridExtra::grid.arrange(p1, p2, ncol = 2)
# График 3: Диаграмма рассеяния — длина vs ширина лепестка
ggplot(iris, aes(x = Petal.Length, y = Petal.Width, color = Species)) +
geom_point(size = 2, alpha = 0.8) +
labs(title = "Связь между длиной и шириной лепестка",
x = "Длина лепестка", y = "Ширина лепестка") +
theme_minimal()
Вывод по EDA: Датасет iris содержит 150 наблюдений без пропусков. Видно, что setosa явно отличается от двух других видов по размеру лепестков. Существует сильная линейная связь между длиной и шириной лепестка.
2. Бутстрап на языке Python
Используем библиотеку
reticulate, которая позволяет запускать Python прямо внутри R Markdown!
# Подключаем библиотеку для запуска Python в R
library(reticulate)
# Если Python не настроен, раскомментируй и укажи путь:
# use_python("C:/Users/ИмяПользователя/AppData/Local/Programs/Python/Python311/python.exe")# ============================================================
# БУТСТРАП НА PYTHON (адаптация слайда 19)
# Задача: сравниваем 90-й перцентиль длины лепестка
# у двух видов ириса: versicolor и virginica
# ============================================================
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from sklearn.datasets import load_iris
import pandas as pd
# --- Загружаем данные ---
iris_data = load_iris()
df = pd.DataFrame(iris_data.data, columns=iris_data.feature_names)
df['species'] = iris_data.target
# Длина лепестка: индекс 2 (petal length)
# Вид 1 = versicolor, Вид 2 = virginica
values_a = df[df['species'] == 1]['petal length (cm)'].values # versicolor
values_b = df[df['species'] == 2]['petal length (cm)'].values # virginica
print(f"Versicolor: {len(values_a)} наблюдений, среднее = {values_a.mean():.2f}")## Versicolor: 50 наблюдений, среднее = 4.26print(f"Virginica: {len(values_b)} наблюдений, среднее = {values_b.mean():.2f}")## Virginica: 50 наблюдений, среднее = 5.55# --- Функция построения доверительного интервала ---
def get_percentile_ci(bootstrap_stats, pe, alpha):
"""Строит перцентильный доверительный интервал."""
left, right = np.quantile(bootstrap_stats, [alpha / 2, 1 - alpha / 2])
return left, right
# --- Параметры бутстрапа ---
n = 1000 # размер каждой бутстрап-выборки
B = 10000 # количество итераций бутстрапа
alpha = 0.05 # уровень значимости (95% ДИ)
# --- Вычисляем точечную оценку (разница 90-х перцентилей) ---
pe = np.quantile(values_b, 0.9) - np.quantile(values_a, 0.9)
print(f"\nРазница 90-х перцентилей (virginica - versicolor): {pe:.4f}")##
## Разница 90-х перцентилей (virginica - versicolor): 1.5100# --- Бутстрап-цикл ---
# Для каждой итерации: берём случайные выборки с возвращением
# и считаем разницу перцентилей
bootstrap_metrics_a = np.quantile(
np.random.choice(values_a, (B, n), True), 0.9, axis=1
)
bootstrap_metrics_b = np.quantile(
np.random.choice(values_b, (B, n), True), 0.9, 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]) # 0 не попадает в ДИ → эффект значим
print(f"\nЗначение 90% квантиля изменилось на: {pe:.2f}")##
## Значение 90% квантиля изменилось на: 1.51print(f"{(1 - alpha) * 100:.0f}% доверительный интервал: ({ci[0]:.2f}, {ci[1]:.2f})")## 95% доверительный интервал: (1.40, 1.61)print(f"Отличия статистически значимые: {has_effect}")## Отличия статистически значимые: True# ============================================================
# ВИЗУАЛИЗАЦИЯ
# ============================================================
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
fig.suptitle('Бутстрап: сравнение 90-го перцентиля длины лепестка', fontsize=14)
# --- График 1: Гистограмма бутстрап-распределения ---
axes[0].hist(bootstrap_stats, bins=50, color='steelblue', alpha=0.7, 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('Разница 90-х перцентилей')
axes[0].set_ylabel('Частота')
axes[0].legend(fontsize=8)
# --- График 2: Боксплоты исходных данных ---
axes[1].boxplot([values_a, values_b], labels=['Versicolor', 'Virginica'], patch_artist=True,
boxprops=dict(facecolor='lightblue', alpha=0.7))## {'whiskers': [<matplotlib.lines.Line2D object at 0x0000022E56B01810>, <matplotlib.lines.Line2D object at 0x0000022E56B01950>, <matplotlib.lines.Line2D object at 0x0000022E56B020D0>, <matplotlib.lines.Line2D object at 0x0000022E56B02210>], 'caps': [<matplotlib.lines.Line2D object at 0x0000022E56B01A90>, <matplotlib.lines.Line2D object at 0x0000022E56B01BD0>, <matplotlib.lines.Line2D object at 0x0000022E56B02350>, <matplotlib.lines.Line2D object at 0x0000022E56B02490>], 'boxes': [<matplotlib.patches.PathPatch object at 0x0000022E549F63C0>, <matplotlib.patches.PathPatch object at 0x0000022E56B01F90>], 'medians': [<matplotlib.lines.Line2D object at 0x0000022E56B01D10>, <matplotlib.lines.Line2D object at 0x0000022E56B025D0>], 'fliers': [<matplotlib.lines.Line2D object at 0x0000022E56B01E50>, <matplotlib.lines.Line2D object at 0x0000022E56B02710>], 'means': []}axes[1].set_title('Длина лепестка\nпо видам')
axes[1].set_ylabel('Длина лепестка (см)')
axes[1].axhline(np.quantile(values_a, 0.9), color='blue', linestyle='--', alpha=0.5, label='90-й %')
axes[1].axhline(np.quantile(values_b, 0.9), color='red', linestyle='--', alpha=0.5)
axes[1].legend(fontsize=8)
# --- График 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)
axes[2].axvline(ci[0], color='orange', linestyle='--', label=f'2.5% = {ci[0]:.2f}')
axes[2].axvline(ci[1], color='orange', linestyle='--', label=f'97.5% = {ci[1]:.2f}')
axes[2].axvline(0, color='black', linestyle=':', label='Ноль')
axes[2].set_title('Кумулятивное распределение\nбутстрап-статистики')
axes[2].set_xlabel('Разница перцентилей')
axes[2].set_ylabel('Вероятность')
axes[2].legend(fontsize=8)
plt.tight_layout()
plt.savefig('bootstrap_python.png', dpi=150, bbox_inches='tight')
plt.show()
print("График сохранён!")## График сохранён!Вывод по Python-бутстрапу: У virginica 90-й перцентиль длины лепестка значительно больше, чем у versicolor. Доверительный интервал не включает 0, поэтому разница статистически значима.
3. Бутстрап на языке R
# Устанавливаем нужные пакеты (только один раз!)
install.packages("boot")
install.packages("ggplot2")
install.packages("gridExtra")# ============================================================
# БУТСТРАП НА R (адаптация слайда 21)
# Задача: бутстрап для коэффициента детерминации R²
# модели: Длина лепестка ~ Ширина лепестка (данные iris)
# ============================================================
library(boot)
library(ggplot2)
library(gridExtra)
# --- Шаг 1: Определяем функцию-статистику для бутстрапа ---
# Эта функция принимает данные и индексы бутстрап-выборки,
# обучает линейную регрессию и возвращает R²
rsq_function <- function(data, indices) {
d <- data[indices, ] # берём подвыборку по индексам
fit <- lm(Petal.Length ~ Petal.Width, data = d) # линейная регрессия
return(summary(fit)$r.square) # возвращаем R²
}
# --- Шаг 2: Запускаем бутстрап ---
set.seed(42) # фиксируем случайность для воспроизводимости
reps <- boot(data = iris,
statistic = rsq_function,
R = 2000) # 2000 бутстрап-повторений
# --- Шаг 3: Смотрим результаты ---
print(reps)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = iris, statistic = rsq_function, R = 2000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.9271098 0.0005393181 0.00994126# --- Шаг 4: Доверительный интервал методом BCa ---
ci_result <- boot.ci(reps, type = "bca")
print(ci_result)## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = reps, type = "bca")
##
## Intervals :
## Level BCa
## 95% ( 0.9035, 0.9435 )
## Calculations and Intervals on Original Scalecat("\n--- Интерпретация ---\n")##
## --- Интерпретация ---cat("Оригинальное R²:", round(reps$t0, 4), "\n")## Оригинальное R²: 0.9271cat("95% ДИ (BCa):", round(ci_result$bca[4], 4), "-", round(ci_result$bca[5], 4), "\n")## 95% ДИ (BCa): 0.9035 - 0.9435Визуализация бутстрапа на R
# Преобразуем бутстрап-оценки в датафрейм
boot_df <- data.frame(r_squared = reps$t)
# Извлекаем границы доверительного интервала
ci_low <- ci_result$bca[4]
ci_high <- ci_result$bca[5]
# --- График 1: Гистограмма бутстрап-распределения R² ---
p1 <- ggplot(boot_df, aes(x = r_squared)) +
geom_histogram(bins = 50, fill = "steelblue", alpha = 0.7, color = "white") +
geom_vline(xintercept = reps$t0, color = "red", linewidth = 1.2,
linetype = "solid") +
geom_vline(xintercept = ci_low, color = "orange", linewidth = 1,
linetype = "dashed") +
geom_vline(xintercept = ci_high, color = "orange", linewidth = 1,
linetype = "dashed") +
annotate("text", x = reps$t0, y = Inf, label = paste("R²=", round(reps$t0,3)),
vjust = 2, hjust = -0.1, color = "red", size = 3.5) +
labs(title = "Бутстрап-распределение R²",
subtitle = paste0("95% ДИ (BCa): [", round(ci_low,3), "; ", round(ci_high,3), "]"),
x = "R² (коэффициент детерминации)", y = "Частота") +
theme_minimal()
# --- График 2: Q-Q plot (нормальность бутстрап-распределения) ---
p2 <- 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: Исходные данные с линией регрессии ---
p3 <- ggplot(iris, aes(x = Petal.Width, y = Petal.Length, color = Species)) +
geom_point(alpha = 0.6, size = 2) +
geom_smooth(method = "lm", se = TRUE, color = "black", linewidth = 1) +
labs(title = "Регрессия: Длина ~ Ширина лепестка",
subtitle = paste("R² =", round(reps$t0, 3)),
x = "Ширина лепестка (см)", y = "Длина лепестка (см)") +
theme_minimal()
# Выводим все три графика
grid.arrange(p1, p2, p3, ncol = 3)
Вывод по R-бутстрапу: Коэффициент детерминации R² ≈ 0.927, что означает: ширина лепестка объясняет около 93% вариации длины лепестка. Доверительный интервал BCa очень узкий и не включает значения ниже 0.9 — модель стабильна и надёжна.
Общие выводы
Бутстрап — мощный непараметрический метод для оценки точности статистик без предположений о распределении данных.
Python-реализация показала, что 90-й перцентиль длины лепестка у virginica статистически значимо больше, чем у versicolor (ДИ не включает 0).
R-реализация подтвердила высокую линейную связь между шириной и длиной лепестка (R² ≈ 0.93). Бутстрап подтвердил устойчивость этой оценки.
Оба языка дают схожие результаты, но R имеет удобный пакет
bootс встроенными методами построения ДИ.