Анализ данных, лабораторная работа №6 “Реализация бутстрапа на языках R и Python”, отчет
Резников Евгений Александрович 231-333
2026-04-15
Бутстрап (bootstrap)
Бутстрап - это компьютерный метод, который позволяет оценить точность ваших статистических выводов, не зная теоретического распределения данных.
Говоря совсем просто: это метод “имитации” повторных экспериментов на основе всего одной имеющейся у вас выборки.
Работа с датасетом iris и разведочный анализ данных на языках R и Python.
options(repos = c(CRAN = "https://cran.rstudio.com/"))Для работы был выбран датасет iris доступный и в R и в Python.
Выполним разведочный анализ данных на Python. Выведем базовые статистики, парные диаграммы рассеяния и диаграмму корреляции признаков.
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
# Загружаем данные
iris = load_iris()
df = pd.DataFrame(data=iris.data, columns=iris.feature_names)
df['species'] = iris.target
df['species'] = df['species'].map({0: 'setosa', 1: 'versicolor', 2: 'virginica'})
df = df.rename(columns={
'sepal length (cm)': 'Длина чашелистика (см)',
'sepal width (cm)': 'Ширина чашелистика (см)',
'petal length (cm)': 'Длина лепестка (см)',
'petal width (cm)': 'Ширина лепестка (см)'
})
df['species'] = df['species'].map({
'setosa': 'Щетинистый ирис',
'versicolor': 'Разноцветный ирис',
'virginica': 'Виргинский ирис'
})
print("Основные статистики:")## Основные статистики:print(df.describe())## Длина чашелистика (см) ... Ширина лепестка (см)
## count 150.000000 ... 150.000000
## mean 5.843333 ... 1.199333
## std 0.828066 ... 0.762238
## min 4.300000 ... 0.100000
## 25% 5.100000 ... 0.300000
## 50% 5.800000 ... 1.300000
## 75% 6.400000 ... 1.800000
## max 7.900000 ... 2.500000
##
## [8 rows x 4 columns]# Создаем фигуру с несколькими подграфиками
fig = plt.figure(figsize=(15, 12))
# ГРАФИК 1: Ящик с усами (Box plot)
plt.subplot(2, 2, 1)
sns.boxplot(x='species', y='Длина чашелистика (см)', data=df, palette='Set2')
plt.title('Распределение длины чашелистика по видам', fontsize=12, fontweight='bold')
plt.xlabel('Вид ириса', fontsize=10)
plt.ylabel('Длина чашелистика (см)', fontsize=10)
plt.xticks(rotation=45)## ([0, 1, 2], [Text(0, 0, 'Щетинистый ирис'), Text(1, 0, 'Разноцветный ирис'), Text(2, 0, 'Виргинский ирис')])# ГРАФИК 2: Скрипичный график (Violin plot)
plt.subplot(2, 2, 2)
sns.violinplot(x='species', y='Длина лепестка (см)', data=df, palette='Set3')
plt.title('Распределение длины лепестка по видам', fontsize=12, fontweight='bold')
plt.xlabel('Вид ириса', fontsize=10)
plt.ylabel('Длина лепестка (см)', fontsize=10)
plt.xticks(rotation=45)## ([0, 1, 2], [Text(0, 0, 'Щетинистый ирис'), Text(1, 0, 'Разноцветный ирис'), Text(2, 0, 'Виргинский ирис')])# ГРАФИК 3: Ящик с усами + точки
plt.subplot(2, 2, 3)
sns.boxplot(x='species', y='Ширина чашелистика (см)', data=df, palette='pastel')
sns.swarmplot(x='species', y='Ширина чашелистика (см)', data=df, color='black', alpha=0.5, size=3)
plt.title('Ширина чашелистика: боксплот с точками', fontsize=12, fontweight='bold')
plt.xlabel('Вид ириса', fontsize=10)
plt.ylabel('Ширина чашелистика (см)', fontsize=10)
plt.xticks(rotation=45)## ([0, 1, 2], [Text(0, 0, 'Щетинистый ирис'), Text(1, 0, 'Разноцветный ирис'), Text(2, 0, 'Виргинский ирис')])# ГРАФИК 4: Корреляционная матрица
plt.subplot(2, 2, 4)
# Переименовываем для корреляционной матрицы
corr_df = df[['Длина чашелистика (см)', 'Ширина чашелистика (см)',
'Длина лепестка (см)', 'Ширина лепестка (см)']].corr()
sns.heatmap(corr_df, annot=True, cmap='coolwarm', center=0,
square=True, linewidths=0.5, fmt='.2f')
plt.title('Корреляция между признаками', fontsize=12, fontweight='bold')
plt.suptitle('Анализ данных Iris (Ирисы Фишера)', fontsize=16, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()
Выполним аналогичный разведочный анализ данных на R. Выведем базовые статистики, парные диаграммы рассеяния и диаграмму корреляции признаков.
install.packages("corrplot")## Устанавливаю пакет в 'C:/Users/rezni/AppData/Local/R/win-library/4.4'
## (потому что 'lib' не определено)## пакет 'corrplot' успешно распакован, MD5-суммы проверены
##
## Скачанные бинарные пакеты находятся в
## C:\Users\rezni\AppData\Local\Temp\Rtmpa4ahB9\downloaded_packages# Загружаем данные
data(iris)
# Переименовываем всё на русский
names(iris) <- c("Длина_чашелистика", "Ширина_чашелистика",
"Длина_лепестка", "Ширина_лепестка", "Вид")
# Меняем названия видов
levels(iris$Вид) <- c("Щетинистый", "Разноцветный", "Виргинский")
# Статистики
summary(iris)## Длина_чашелистика Ширина_чашелистика Длина_лепестка Ширина_лепестка
## 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
## Вид
## Щетинистый :50
## Разноцветный:50
## Виргинский :50
##
##
## # Парные диаграммы
pairs(iris[,1:4], # ("Щетинистый", "Разноцветный", "Виргинский")
col = c("red", "green3", "blue")[iris$Вид],
pch = 19,
main = "Диаграммы рассеяния для Iris",
labels = c("Длина чашелистика", "Ширина чашелистика",
"Длина лепестка", "Ширина лепестка"))
# Корреляционная матрица
library(corrplot)## corrplot 0.95 loadedcor_matrix <- cor(iris[,1:4])
rownames(cor_matrix) <- colnames(cor_matrix) <-
c("Длина чашелистика", "Ширина чашелистика", "Длина лепестка", "Ширина лепестка")
corrplot(cor_matrix, method = "color", addCoef.col = "black",
tl.col = "black", tl.srt = 45)
Реализация алгоритма бутстрапа на Python.
# Импорт необходимых библиотек
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
def get_percentile_ci(bootstrap_stats, pe, alpha=0.05):
"""Функция для вычисления процентильного доверительного интервала"""
left, right = np.quantile(bootstrap_stats, [alpha / 2, 1 - alpha / 2])
return left, right
# Загружаем данные iris
data = sns.load_dataset("iris")
n = len(data) # Размер выборки
B = 5000 # Количество бутстрап-итераций
alpha = 0.05 # Уровень значимости
# Выделяем длину чашелистика для двух видов ирисов
setosa = data[data.species == "setosa"]["sepal_length"]
versicolor = data[data.species == "versicolor"]["sepal_length"]
# Точечная оценка - разница средних между видами
pe = setosa.mean() - versicolor.mean()
# Бутстрап для setosa: генерируем B выборок с возвращением и считаем средние
bootstrap_setosa = np.random.choice(setosa, (B, n), replace=True).mean(axis=1)
# Бутстрап для versicolor
bootstrap_versicolor = np.random.choice(versicolor, (B, n), replace=True).mean(axis=1)
# Бутстрап-распределение разницы средних
bootstrap_stats = bootstrap_setosa - bootstrap_versicolor
# Вычисляем доверительный интервал
ci = get_percentile_ci(bootstrap_stats, pe, alpha)
# Проверяем статистическую значимость (если 0 не входит в ДИ)
has_effect = not (ci[0] < 0 < ci[1])
# Вывод результатов
print("=" * 60)## ============================================================print("РЕЗУЛЬТАТЫ БУТСТРАП-АНАЛИЗА")## РЕЗУЛЬТАТЫ БУТСТРАП-АНАЛИЗАprint("=" * 60)## ============================================================print(f"Разница средних (Щетинистый - Разноцветный): {pe:.3f} см")## Разница средних (Щетинистый - Разноцветный): -0.930 смprint(f"{((1 - alpha) * 100)}% доверительный интервал: ({ci[0]:.3f}, {ci[1]:.3f}) см")## 95.0% доверительный интервал: (-1.031, -0.831) смprint(f"Различия статистически значимы (0 не входит в ДИ): {has_effect}")## Различия статистически значимы (0 не входит в ДИ): Trueprint("=" * 60)## ============================================================# СОЗДАЕМ 2 ГРАФИКА
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
# ГРАФИК 1: Гистограмма с плотностью
sns.histplot(bootstrap_stats, bins=50, kde=True, color='purple',
alpha=0.6, edgecolor='white', ax=ax1)
# Закрашиваем область доверительного интервала
ax1.axvspan(ci[0], ci[1], alpha=0.2, color='orange', label='95% доверительный интервал')
# Вертикальные линии
ax1.axvline(pe, color='darkred', linestyle='-', linewidth=2.5,
label=f'Точечная оценка: {pe:.3f} см')
ax1.axvline(ci[0], color='orange', linestyle='--', linewidth=2,
label=f'Нижняя граница: {ci[0]:.3f} см')
ax1.axvline(ci[1], color='orange', linestyle='--', linewidth=2,
label=f'Верхняя граница: {ci[1]:.3f} см')
ax1.axvline(0, color='black', linestyle=':', linewidth=2,
label='Нулевое значение')
ax1.set_title('Бутстрап-распределение разницы средних', fontsize=12, fontweight='bold')
ax1.set_xlabel('Разница средних (Щетинистый - Разноцветный), см', fontsize=10)
ax1.set_ylabel('Частота', fontsize=10)
ax1.legend(loc='upper right', fontsize=9)
ax1.grid(alpha=0.3, linestyle='--')
# Добавляем текст на первый график
if has_effect:
ax1.text(0.02, 0.95, 'Статистически значимое различие',
transform=ax1.transAxes, fontsize=10,
bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.7))
# ГРАФИК 2: Ящик с усами (Boxplot)
# Создаем данные для boxplot (сравнение исходных выборок)
box_data = [setosa, versicolor]
bp = ax2.boxplot(box_data, labels=['Щетинистый', 'Разноцветный'],
patch_artist=True, widths=0.6)## <string>:1: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11.# Настройка цветов для boxplot
colors = ['lightcoral', 'lightblue']
for patch, color in zip(bp['boxes'], colors):
patch.set_facecolor(color)
patch.set_alpha(0.7)
# Настройка линий
bp['medians'][0].set_color('darkred')
bp['medians'][0].set_linewidth(2)
bp['medians'][1].set_color('darkblue')
bp['medians'][1].set_linewidth(2)
# Добавляем точки (все наблюдения)
ax2.scatter(np.random.normal(1, 0.04, len(setosa)), setosa,
alpha=0.5, color='darkred', s=20, label='Щетинистый')
ax2.scatter(np.random.normal(2, 0.04, len(versicolor)), versicolor,
alpha=0.5, color='darkblue', s=20, label='Разноцветный')
ax2.set_title('Сравнение исходных выборок\n(длина чашелистика)', fontsize=12, fontweight='bold')
ax2.set_xlabel('Вид ириса', fontsize=10)
ax2.set_ylabel('Длина чашелистика (см)', fontsize=10)
ax2.legend(loc='upper right', fontsize=9)
ax2.grid(alpha=0.3, axis='y', linestyle='--')
# Общий заголовок
plt.suptitle('Бутстрап-анализ: сравнение длины чашелистика у двух видов ирисов',
fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
# Дополнительная статистика на русском
print("\n" + "=" * 60)##
## ============================================================print("ДОПОЛНИТЕЛЬНАЯ СТАТИСТИКА")## ДОПОЛНИТЕЛЬНАЯ СТАТИСТИКАprint("=" * 60)## ============================================================print(f"Стандартная ошибка (SE): {bootstrap_stats.std():.3f} см")## Стандартная ошибка (SE): 0.051 смprint(f"Смещение (bias): {bootstrap_stats.mean() - pe:.3f} см")## Смещение (bias): 0.001 смprint(f"Коэффициент вариации: {bootstrap_stats.std() / abs(pe):.3f}")## Коэффициент вариации: 0.055print(f"Медиана бутстрап-распределения: {np.median(bootstrap_stats):.3f} см")## Медиана бутстрап-распределения: -0.929 см# Статистика по исходным данным
print("\nСТАТИСТИКА ПО ИСХОДНЫМ ВЫБОРКАМ:")##
## СТАТИСТИКА ПО ИСХОДНЫМ ВЫБОРКАМ:print(f"Щетинистый ирис: среднее = {setosa.mean():.3f} см, std = {setosa.std():.3f} см")## Щетинистый ирис: среднее = 5.006 см, std = 0.352 смprint(f"Разноцветный ирис: среднее = {versicolor.mean():.3f} см, std = {versicolor.std():.3f} см")## Разноцветный ирис: среднее = 5.936 см, std = 0.516 смРеализация алгоритма бутстрапа на R.
# Устанавливаем seed для воспроизводимости результатов
set.seed(123)
# Загружаем необходимые библиотеки
library(boot)
library(ggplot2)
# Загружаем данные iris
data(iris)
# Переименовываем столбцы
iris$Species <- factor(iris$Species,
levels = c("setosa", "versicolor", "virginica"),
labels = c("Щетинистый", "Разноцветный", "Виргинский"))
names(iris)[names(iris) == "Sepal.Length"] <- "Длина.чашелистика"
names(iris)[names(iris) == "Sepal.Width"] <- "Ширина.чашелистика"
names(iris)[names(iris) == "Petal.Length"] <- "Длина.лепестка"
names(iris)[names(iris) == "Petal.Width"] <- "Ширина.лепестка"
# Оставляем только два вида ирисов для сравнения
iris_subset <- subset(iris, Species %in% c("Щетинистый", "Разноцветный"))
# Функция для бутстрапа: вычисляет разницу средних длин чашелистика между видами
bootstrap_function <- function(data, indices) {
d <- data[indices, ]
mean_setosa <- mean(d[d$Species == "Щетинистый", "Длина.чашелистика"])
mean_versicolor <- mean(d[d$Species == "Разноцветный", "Длина.чашелистика"])
return(mean_setosa - mean_versicolor)
}
# Количество бутстрап-итераций
B <- 5000
# Выполняем бутстрап с помощью функции boot()
bootstrap_results <- boot(data = iris_subset,
statistic = bootstrap_function,
R = B)
# Вычисляем доверительный интервал методом BCa
ci <- boot.ci(bootstrap_results, type = "bca")
# Выводим доверительный интервал
print(ci)## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 5000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = bootstrap_results, type = "bca")
##
## Intervals :
## Level BCa
## 95% (-1.1035, -0.7540 )
## Calculations and Intervals on Original Scale# Создаем датафрейм с бутстрап-статистиками для визуализации
bootstrap_df <- data.frame(Diff = bootstrap_results$t)
# Точечная оценка (разница средних в исходной выборке)
point_estimate <- bootstrap_function(iris_subset, 1:nrow(iris_subset))
# Визуализация распределения бутстрап-статистик
ggplot(bootstrap_df, aes(x = Diff)) +
# Гистограмма
geom_histogram(aes(y = after_stat(density)),
bins = 50,
fill = "orange",
color = "white",
alpha = 0.7) +
# Ядерная оценка плотности (новый цвет: темно-красный)
geom_density(color = "darkred", size = 1.2) +
# Вертикальная линия для точечной оценки (новый цвет: синий)
geom_vline(xintercept = point_estimate,
color = "blue",
linetype = "solid",
size = 1.2) +
# Вертикальная линия для нижней границы ДИ (новый цвет: оранжевый)
geom_vline(xintercept = ci$bca[4],
color = "orange",
linetype = "dashed",
size = 1.2) +
# Вертикальная линия для верхней границы ДИ (новый цвет: оранжевый)
geom_vline(xintercept = ci$bca[5],
color = "orange",
linetype = "dashed",
size = 1.2) +
# Вертикальная линия для нулевого значения
geom_vline(xintercept = 0,
color = "black",
linetype = "dotted",
size = 0.8) +
# Подписи и заголовок (на русском)
labs(title = "Бутстрап-распределение разницы средних длин чашелистика",
subtitle = paste0("Щетинистый vs Разноцветный (B = ", B, " итераций)"),
x = "Разница средних (Щетинистый - Разноцветный), см",
y = "Плотность") +
# Тема оформления
theme_minimal() +
# Добавляем подписи к линиям (на русском)
annotate("text",
x = point_estimate + 0.05,
y = 2,
label = "Точечная оценка",
color = "blue",
hjust = 0,
size = 3.5) +
annotate("text",
x = ci$bca[4] - 0.05,
y = 1.8,
label = "Нижняя граница ДИ",
color = "orange",
hjust = 1,
size = 3.5) +
annotate("text",
x = ci$bca[5] + 0.05,
y = 1.8,
label = "Верхняя граница ДИ",
color = "orange",
hjust = 0,
size = 3.5)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Вывод дополнительной статистики
cat("РЕЗУЛЬТАТЫ БУТСТРАП-АНАЛИЗА\n")## РЕЗУЛЬТАТЫ БУТСТРАП-АНАЛИЗАcat("Точечная оценка (разница средних):", round(point_estimate, 3), "см\n")## Точечная оценка (разница средних): -0.93 смcat("Стандартная ошибка (SE):", round(sd(bootstrap_results$t), 3), "см\n")## Стандартная ошибка (SE): 0.088 смcat("Смещение (bias):", round(mean(bootstrap_results$t) - point_estimate, 3), "см\n")## Смещение (bias): -0.001 смcat("95% доверительный интервал (BCa): [", round(ci$bca[4], 3), ", ", round(ci$bca[5], 3), "] см\n")## 95% доверительный интервал (BCa): [ -1.103 , -0.754 ] см# Проверка статистической значимости
has_effect <- !(ci$bca[4] < 0 & ci$bca[5] > 0)
cat("Статистически значимое различие (0 не входит в ДИ):", has_effect, "\n")## Статистически значимое различие (0 не входит в ДИ): TRUEВывод
В результате выполнения данной лабораторной работы были закреплены навыки в работе с языками R и Python. Также была проведена работа с бутстрапом.