Введение

Бутстрап — это компьютерный статистический метод, который позволяет оценивать распределение статистик (среднее, медиана, доверительные интервенты) многократной генерацией выборок из имеющихся данных. Преимущество метода — не требуется предполагать нормальность распределения и можно использовать на небольших выборках.

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

Для анализа используется датасет `iris`. Python и R части будут работать с разными признаками и разными видами ирисов, чтобы результаты выглядели уникально.

options(repos = c(CRAN = "https://cran.rstudio.com/"))

# Настройка Python для RMarkdown через reticulate
library(reticulate)
use_python("C:/Users/XataB/AppData/Local/Programs/Python/Python312/python.exe", required=TRUE)

Разведочный анализ данных

Python: Scatter plot и корреляция

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris

# Загружаем данные 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'})

# Используем petal length и petal width
df_subset = df[['petal length (cm)', 'petal width (cm)', 'species']]

print("Основные статистики по выбранным признакам:")
## Основные статистики по выбранным признакам:
print(df_subset.describe())
##        petal length (cm)  petal width (cm)
## count         150.000000        150.000000
## mean            3.758000          1.199333
## std             1.765298          0.762238
## min             1.000000          0.100000
## 25%             1.600000          0.300000
## 50%             4.350000          1.300000
## 75%             5.100000          1.800000
## max             6.900000          2.500000
# Парные диаграммы рассеяния
sns.pairplot(df_subset, hue='species')

plt.show()

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

R: Scatter plot и корреляция

install.packages("ggplot2")
## Устанавливаю пакет в 'C:/Users/XataB/AppData/Local/R/win-library/4.5'
## (потому что 'lib' не определено)
## пакет 'ggplot2' успешно распакован, MD5-суммы проверены
## 
## Скачанные бинарные пакеты находятся в
##  C:\Users\XataB\AppData\Local\Temp\RtmpghhidW\downloaded_packages
install.packages("GGally")
## Устанавливаю пакет в 'C:/Users/XataB/AppData/Local/R/win-library/4.5'
## (потому что 'lib' не определено)
## пакет 'GGally' успешно распакован, MD5-суммы проверены
## 
## Скачанные бинарные пакеты находятся в
##  C:\Users\XataB\AppData\Local\Temp\RtmpghhidW\downloaded_packages
library(ggplot2)
library(GGally)

data(iris)
# Используем другие виды: versicolor и virginica
iris_subset <- subset(iris, Species %in% c("versicolor", "virginica"))
summary(iris_subset)
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.900   Min.   :2.000   Min.   :3.000   Min.   :1.000  
##  1st Qu.:5.800   1st Qu.:2.700   1st Qu.:4.375   1st Qu.:1.300  
##  Median :6.300   Median :2.900   Median :4.900   Median :1.600  
##  Mean   :6.262   Mean   :2.872   Mean   :4.906   Mean   :1.676  
##  3rd Qu.:6.700   3rd Qu.:3.025   3rd Qu.:5.525   3rd Qu.:2.000  
##  Max.   :7.900   Max.   :3.800   Max.   :6.900   Max.   :2.500  
##        Species  
##  setosa    : 0  
##  versicolor:50  
##  virginica :50  
##                 
##                 
## 
# Pair plot через ggpairs
ggpairs(iris_subset[,c("Sepal.Length","Sepal.Width","Petal.Length","Petal.Width","Species")],
        aes(color=Species, alpha=0.5))
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

Бутстрап на Python

import numpy as np

def get_percentile_ci(bootstrap_stats, alpha=0.05):
    left, right = np.quantile(bootstrap_stats, [alpha/2, 1-alpha/2])
    return left, right

# Длина лепестка для versicolor и virginica
versicolor = df[df.species=="versicolor"]['petal length (cm)']
virginica = df[df.species=="virginica"]['petal length (cm)']

# Точечная оценка медианы
pe = np.median(versicolor) - np.median(virginica)

B = 10000
n_v = len(versicolor)
n_vi = len(virginica)

# Генерация бутстрап-выборок
bootstrap_versicolor = np.median(np.random.choice(versicolor, (B, n_v), replace=True), axis=1)
bootstrap_virginica = np.median(np.random.choice(virginica, (B, n_vi), replace=True), axis=1)

bootstrap_stats = bootstrap_versicolor - bootstrap_virginica
ci = get_percentile_ci(bootstrap_stats)

# Проверка статистической значимости
has_effect = not (ci[0]<0<ci[1])

print(f'Разница медиан (versicolor - virginica): {pe:.3f}')
## Разница медиан (versicolor - virginica): -1.200
print(f'95% ДИ: ({ci[0]:.3f}, {ci[1]:.3f})')
## 95% ДИ: (-1.450, -0.850)
print(f'Статистически значимо: {has_effect}')
## Статистически значимо: True
# Визуализация
plt.figure(figsize=(10,6))
sns.histplot(bootstrap_stats, bins=50, kde=True, color='skyblue', alpha=0.7)
plt.axvline(pe, color='red', linestyle='-', label=f'Точечная оценка: {pe:.3f}')
plt.axvline(ci[0], color='green', linestyle='--', label=f'Нижняя граница 95% ДИ: {ci[0]:.3f}')
plt.axvline(ci[1], color='green', linestyle='--', label=f'Верхняя граница 95% ДИ: {ci[1]:.3f}')
plt.axvline(0, color='black', linestyle=':', label='Нулевая гипотеза')
plt.title('Бутстрап-распределение медиан petal length\nversicolor vs virginica')
plt.xlabel('Разница медиан')
plt.ylabel('Частота')
plt.legend()
plt.show()

Бутстрап на R

set.seed(123)
library(boot)
library(ggplot2)

# Используем Petal.Width и виды versicolor, virginica
iris_sub <- subset(iris, Species %in% c("versicolor", "virginica"))

bootstrap_func <- function(data, indices){
  d <- data[indices, ]
  med_v <- median(d[d$Species=="versicolor", "Petal.Width"])
  med_vi <- median(d[d$Species=="virginica", "Petal.Width"])
  return(med_v - med_vi)
}

B <- 10000
boot_res <- boot(data=iris_sub, statistic=bootstrap_func, R=B)
ci <- boot.ci(boot_res, type="bca")

point_estimate <- bootstrap_func(iris_sub, 1:nrow(iris_sub))
boot_df <- data.frame(Diff=boot_res$t)

ggplot(boot_df, aes(x=Diff)) +
  geom_density(fill="skyblue", alpha=0.7) +
  geom_vline(xintercept=point_estimate, color="purple", linewidth=1.2) +
  geom_vline(xintercept=ci$bca[4], color="red", linetype="dashed", linewidth=1.2) +
  geom_vline(xintercept=ci$bca[5], color="darkgreen", linetype="dashed", linewidth=1.2) +
  geom_vline(xintercept=0, color="black", linetype="dotted") +
  labs(title="Бутстрап-распределение разницы медиан Petal.Width",
       subtitle="Versicolor vs Virginica", x="Разница медиан", y="Плотность") +
  theme_minimal()

cat("Точечная оценка:", round(point_estimate,3), "\n")
## Точечная оценка: -0.7
cat("95% ДИ (BCa): [", round(ci$bca[4],3), ",", round(ci$bca[5],3), "]\n")
## 95% ДИ (BCa): [ -0.9 , -0.6 ]
has_effect <- !(ci$bca[4]<0 & ci$bca[5]>0)
cat("Статистически значимо:", has_effect, "\n")
## Статистически значимо: TRUE

Вывод

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