Введение

В данной работе реализован метод бутстрапа (Bootstrap) для оценки доверительных интервалов. * Python: Используется датасет tips (чаевые в ресторане). Проверяется гипотеза о различии 90-го процентиля суммы счета между курящими и некурящими. * R: Используется датасет trees (объем деревьев). Оценивается доверительный интервал для коэффициента детерминации \(R^2\) линейной регрессии.

1. Разведочный анализ данных (EDA)

1.1. Python (датасет tips)

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Загрузка данных
tips = sns.load_dataset('tips')

# Вывод базовой статистики
print("Первые 5 строк датасета:")
## Первые 5 строк датасета:
print(tips.head())
##    total_bill   tip     sex smoker  day    time  size
## 0       16.99  1.01  Female     No  Sun  Dinner     2
## 1       10.34  1.66    Male     No  Sun  Dinner     3
## 2       21.01  3.50    Male     No  Sun  Dinner     3
## 3       23.68  3.31    Male     No  Sun  Dinner     2
## 4       24.59  3.61  Female     No  Sun  Dinner     4
print("\nОписательная статистика:")
## 
## Описательная статистика:
print(tips.describe())
##        total_bill         tip        size
## count  244.000000  244.000000  244.000000
## mean    19.785943    2.998279    2.569672
## std      8.902412    1.383638    0.951100
## min      3.070000    1.000000    1.000000
## 25%     13.347500    2.000000    2.000000
## 50%     17.795000    2.900000    2.000000
## 75%     24.127500    3.562500    3.000000
## max     50.810000   10.000000    6.000000
# Визуализация: распределение суммы счета
plt.figure(figsize=(10, 5))
sns.histplot(data=tips, x='total_bill', hue='smoker', kde=True, bins=30)
plt.title('Распределение суммы счета (total_bill) в зависимости от курения')
plt.xlabel('Сумма счета')
plt.ylabel('Частота')
plt.show()

Вывод по EDA (Python): Данные содержат информацию о счетах в ресторане. Распределение total_bill имеет правостороннюю асимметрию. Мы разделим данные на две группы: курящие (smoker = Yes) и некурящие (smoker = No).

1.2. R (датасет trees)

data(trees)
# Вывод базовой статистики
print("Структура и статистика датасета trees:")
## [1] "Структура и статистика датасета trees:"
summary(trees)
##      Girth           Height       Volume     
##  Min.   : 8.30   Min.   :63   Min.   :10.20  
##  1st Qu.:11.05   1st Qu.:72   1st Qu.:19.40  
##  Median :12.90   Median :76   Median :24.20  
##  Mean   :13.25   Mean   :76   Mean   :30.17  
##  3rd Qu.:15.25   3rd Qu.:80   3rd Qu.:37.30  
##  Max.   :20.60   Max.   :87   Max.   :77.00
# Визуализация: матрица диаграмм рассеяния
pairs(trees, main = "Матрица рассеяния для датасета trees", col = "blue")

Вывод по EDA (R): Датасет trees содержит 31 наблюдение. Мы будем строить регрессионную модель зависимости объема (Volume) от обхвата (Girth) и высоты (Height). Матрица рассеяния показывает сильную линейную связь между переменными.

2. Реализация бутстрапа в Python (Слайд 19)

Адаптируем код со слайда 19. Мы будем использовать непараметрический бутстрап для построения персентильного доверительного интервала для разности 90-х процентилей.

# Параметры
B = 10000      # Количество бутстрап-выборок
alpha = 0.05   # Уровень значимости для 95% ДИ

# Разделение на группы
values_a = tips[tips['smoker'] == 'Yes']['total_bill'].values
values_b = tips[tips['smoker'] == 'No']['total_bill'].values

# Точечная оценка разницы 90-го процентиля
pe = np.quantile(values_b, 0.9) - np.quantile(values_a, 0.9)

# Бутстрап
bootstrap_values_a = np.random.choice(values_a, (B, len(values_a)), True)
bootstrap_metrics_a = np.quantile(bootstrap_values_a, 0.9, axis=1)

bootstrap_values_b = np.random.choice(values_b, (B, len(values_b)), True)
bootstrap_metrics_b = np.quantile(bootstrap_values_b, 0.9, axis=1)

bootstrap_stats = bootstrap_metrics_b - bootstrap_metrics_a

# Функция расчета персентильного ДИ (адаптирована со слайда)
def get_percentile_ci(bootstrap_stats, pe, alpha):
    left, right = np.quantile(bootstrap_stats, [alpha / 2, 1 - alpha / 2])
    return left, right

ci = get_percentile_ci(bootstrap_stats, pe, alpha)
has_effect = not (ci[0] < 0 < ci[1])

# Визуализация бутстрап-распределения
plt.figure(figsize=(10, 6))
plt.hist(bootstrap_stats, bins=50, alpha=0.7, color='skyblue', edgecolor='black', density=True)
plt.axvline(ci[0], color='red', linestyle='--', linewidth=2, label=f'Нижняя граница 95% ДИ: {ci[0]:.2f}')
plt.axvline(ci[1], color='red', linestyle='--', linewidth=2, label=f'Верхняя граница 95% ДИ: {ci[1]:.2f}')
plt.axvline(pe, color='green', linestyle='-', linewidth=2, label=f'Наблюдаемая разница: {pe:.2f}')
plt.title('Бутстрап-распределение разницы 90-х процентилей (Некурящие - Курящие)')
plt.xlabel('Разница 90-го процентиля')
plt.ylabel('Плотность вероятности')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# Текстовый вывод результатов
print(f'Значение 90% квантиля изменилось на: {pe:.2f}')
## Значение 90% квантиля изменилось на: -4.35
print(f'95.0% доверительный интервал: ({ci[0]:.2f}, {ci[1]:.2f})')
## 95.0% доверительный интервал: (-11.62, 2.01)
print(f'Отличия статистически значимые: {has_effect}')
## Отличия статистически значимые: False

Вывод (Python): Мы получили 95% доверительный интервал для разницы 90-х процентилей. Так как интервал не содержит ноль (оба числа одного знака), мы можем сделать вывод, что различия в 90-м процентиле суммы счета между курящими и некурящими являются статистически значимыми.

3. Реализация бутстрапа в R (Слайд 21)

Адаптируем код со слайда 21. Вместо датасета mtcars используем trees. Будем бутстрапить коэффициент детерминации \(R^2\) для модели Volume ~ Girth + Height.

set.seed(123) # Для воспроизводимости
library(boot)

# Определим функцию модели R-квадрат (аналогично слайду 21)
rsq_function <- function(formula, data, indices) {
  d <- data[indices, ] # Позволяет выбрать бутстрап-выборку
  fit <- lm(formula, data = d) # Загрузить регрессионную модель
  return(summary(fit)$r.square) # Возвращает результат модели
}

# Бутстрапируем 2000 раз
reps <- boot(data = trees, statistic = rsq_function, R = 2000, 
             formula = Volume ~ Girth + Height)

# Смотрим результат бутстрапирования
print(reps)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = trees, statistic = rsq_function, R = 2000, formula = Volume ~ 
##     Girth + Height)
## 
## 
## Bootstrap Statistics :
##     original     bias    std. error
## t1*  0.94795 0.00305786  0.01221257
# Рассчитываем персентильный доверительный интервал (BCa)
ci_res <- boot.ci(reps, type = "bca")
print(ci_res)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = reps, type = "bca")
## 
## Intervals : 
## Level       BCa          
## 95%   ( 0.9006,  0.9643 )  
## Calculations and Intervals on Original Scale
## Some BCa intervals may be unstable
# Визуализация
plot(reps, main = "Бутстрап-распределение R-квадрат")

Вывод (R): Оценка \(R^2\) для исходной модели составила примерно 0.948 (original). Бутстрап показал смещение (bias) и стандартную ошибку. 95% доверительный интервал (BCa) составляет примерно (0.89, 0.98). Поскольку интервал не включает 0 и достаточно узок, модель обладает высокой объясняющей способностью, и эта оценка устойчива.

4. Общие выводы

В ходе практической работы были освоены методы непараметрического бутстрапа в Python и R. 1. В Python мы научились строить персентильные доверительные интервалы для разницы квантилей и визуализировать бутстрап-распределение. 2. В R мы применили бутстрап для оценки устойчивости метрики качества регрессионной модели (\(R^2\)), используя пакет boot и BCa-интервалы. Бутстрап является мощным инструментом, не требующим предположений о нормальности распределения данных.