Бутстрап (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'})
print("Основные статистики:")
## Основные статистики:
print(df.describe())
##        sepal length (cm)  sepal width (cm)  petal length (cm)  petal width (cm)
## count         150.000000        150.000000         150.000000        150.000000
## mean            5.843333          3.057333           3.758000          1.199333
## std             0.828066          0.435866           1.765298          0.762238
## min             4.300000          2.000000           1.000000          0.100000
## 25%             5.100000          2.800000           1.600000          0.300000
## 50%             5.800000          3.000000           4.350000          1.300000
## 75%             6.400000          3.300000           5.100000          1.800000
## max             7.900000          4.400000           6.900000          2.500000
sns.pairplot(df, hue='species')

plt.show()

plt.figure(figsize=(8, 6))
sns.heatmap(df.drop('species', axis=1).corr(), annot=True, cmap='coolwarm')
plt.title("Корреляционная матрица")
plt.show()

Выполним аналогичный разведочный анализ данных на R. Выведем базовые статистики, парные диаграммы рассеяния и диаграмму корреляции признаков.

install.packages("corrplot")
## package 'corrplot' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\alexa\AppData\Local\Temp\RtmpcFKFYk\downloaded_packages
data(iris)
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  
##                 
##                 
## 
pairs(iris[,1:4], col=iris$Species, pch=19)

cor_matrix <- cor(iris[,1:4])
library(corrplot)
## corrplot 0.95 loaded
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 = 10000      # Количество бутстрап-итераций
alpha = 0.05   # Уровень значимости

# Выделяем длину чашелистика (sepal_length) для двух видов ирисов
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(f'Разница средних (setosa - versicolor): {pe:.3f}')
## Разница средних (setosa - versicolor): -0.930
print(f'{((1 - alpha) * 100)}% доверительный интервал: ({ci[0]:.3f}, {ci[1]:.3f})')
## 95.0% доверительный интервал: (-1.030, -0.834)
print(f'Различия статистически значимы (0 не входит в ДИ): {has_effect}')
## Различия статистически значимы (0 не входит в ДИ): True
# Визуализация бутстрап-распределения
plt.figure(figsize=(12, 7))

# Гистограмма с ядерной оценкой плотности
sns.histplot(bootstrap_stats, bins=50, kde=True, color='skyblue', 
             alpha=0.7, label='Бутстрап-распределение')

# Вертикальные линии для точечной оценки и границ ДИ
plt.axvline(pe, color='red', linestyle='-', linewidth=2, 
            label=f'Точечная оценка: {pe:.3f}')
plt.axvline(ci[0], color='green', linestyle='--', linewidth=2, 
            label=f'Нижняя граница 95% ДИ: {ci[0]:.3f}')
plt.axvline(ci[1], color='green', linestyle='--', linewidth=2, 
            label=f'Верхняя граница 95% ДИ: {ci[1]:.3f}')

# Добавляем вертикальную линию для 0 (нулевая гипотеза)
plt.axvline(0, color='black', linestyle=':', linewidth=1.5, 
            label='Нулевое значение (0)')

plt.title('Распределение бутстрап-статистики (разница средних длин чашелистика)\n' + 
          'setosa vs versicolor', fontsize=14)
plt.xlabel('Разница средних (setosa - versicolor)', fontsize=12)
plt.ylabel('Частота', fontsize=12)
plt.legend(loc='upper right')
plt.grid(alpha=0.3)

# Добавляем текст с интерпретацией на график
if has_effect:
    plt.text(0.02, 0.95, 'Статистически значимое различие', 
             transform=plt.gca().transAxes, fontsize=12,
             bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.7))
else:
    plt.text(0.02, 0.95, 'Статистически незначимое различие', 
             transform=plt.gca().transAxes, fontsize=12,
             bbox=dict(boxstyle='round', facecolor='salmon', alpha=0.7))

plt.tight_layout()
plt.show()

# Дополнительная статистика
print(f"\nДополнительная информация:")
## 
## Дополнительная информация:
print(f"Стандартная ошибка (SE): {bootstrap_stats.std():.3f}")
## Стандартная ошибка (SE): 0.050
print(f"Смещение (bias): {bootstrap_stats.mean() - pe:.3f}")
## Смещение (bias): -0.001
print(f"Коэффициент вариации: {bootstrap_stats.std() / abs(pe):.3f}")
## Коэффициент вариации: 0.054

Реализация алгоритма бутстрапа на R.

# Устанавливаем seed для воспроизводимости результатов
set.seed(123)

# Загружаем необходимые библиотеки
library(boot)
library(ggplot2)

# Загружаем данные iris
data(iris)

# Оставляем только два вида ирисов для сравнения
iris_subset <- subset(iris, Species %in% c("setosa", "versicolor"))

# Функция для бутстрапа: вычисляет разницу средних длин чашелистика между видами
# data - исходные данные
# indices - индексы для бутстрап-выборки (с возвращением)
bootstrap_function <- function(data, indices) {
  # Формируем бутстрап-выборку по переданным индексам
  d <- data[indices, ]
  
  mean_setosa <- mean(d[d$Species == "setosa", "Sepal.Length"])
  
  mean_versicolor <- mean(d[d$Species == "versicolor", "Sepal.Length"])

  return(mean_setosa - mean_versicolor)
}

# Количество бутстрап-итераций
B <- 10000

# Выполняем бутстрап с помощью функции boot()
# statistic - наша функция, R - количество повторений
bootstrap_results <- boot(data = iris_subset, 
                          statistic = bootstrap_function, 
                          R = B)

# Вычисляем доверительный интервал методом BCa (bias-corrected and accelerated)
# Это один из самых точных методов для бутстрап-доверительных интервалов
ci <- boot.ci(bootstrap_results, type = "bca")

# Выводим доверительный интервал
print(ci)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = bootstrap_results, type = "bca")
## 
## Intervals : 
## Level       BCa          
## 95%   (-1.1002, -0.7600 )  
## 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 = "skyblue", 
                 color = "black", 
                 alpha = 0.7) +
  # Ядерная оценка плотности
  geom_density(color = "darkblue", size = 1.2) +
  # Вертикальная линия для точечной оценки
  geom_vline(xintercept = point_estimate, 
             color = "purple", 
             linetype = "solid", 
             size = 1.2) +
  # Вертикальная линия для нижней границы ДИ
  geom_vline(xintercept = ci$bca[4], 
             color = "red", 
             linetype = "dashed", 
             size = 1.2) +
  # Вертикальная линия для верхней границы ДИ
  geom_vline(xintercept = ci$bca[5], 
             color = "darkgreen", 
             linetype = "dashed", 
             size = 1.2) +
  # Вертикальная линия для нулевого значения (проверка гипотезы)
  geom_vline(xintercept = 0, 
             color = "black", 
             linetype = "dotted", 
             size = 0.8) +
  # Подписи и заголовок
  labs(title = "Бутстрап-распределение разницы средних длин чашелистика",
       subtitle = paste0("Setosa vs Versicolor (B = ", B, " итераций)"),
       x = "Разница средних (Setosa - Versicolor)", 
       y = "Плотность") +
  # Тема оформления
  theme_minimal() +
  # Добавляем подписи к линиям через annotate
  annotate("text", 
           x = point_estimate + 0.05, 
           y = 2, 
           label = "Точечная оценка", 
           color = "purple", 
           hjust = 0, 
           size = 3.5) +
  annotate("text", 
           x = ci$bca[4] - 0.05, 
           y = 1.8, 
           label = "Нижняя граница", 
           color = "red", 
           hjust = 1, 
           size = 3.5) +
  annotate("text", 
           x = ci$bca[5] + 0.05, 
           y = 1.8, 
           label = "Верхняя граница", 
           color = "darkgreen", 
           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("Точечная оценка (разница средних):", round(point_estimate, 3), "\n")
## Точечная оценка (разница средних): -0.93
cat("Стандартная ошибка (SE):", round(sd(bootstrap_results$t), 3), "\n")
## Стандартная ошибка (SE): 0.087
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.1 ,  -0.76 ]
# Проверка статистической значимости
has_effect <- !(ci$bca[4] < 0 & ci$bca[5] > 0)
cat("Статистически значимое различие (0 не входит в ДИ):", has_effect, "\n")
## Статистически значимое различие (0 не входит в ДИ): TRUE

Вывод

В результате выполнения данной лабораторной работы были закреплены навыки в работе с языками R и Python. Также был изучен процесс бутстрапа.