Анализ данных, лабораторная работа №6. Реализация бутстрапа на языках R и Python

Задание 1

Я возьму датасет mtcars. Установлю необходимые библиотеки и импортирую его.

library(reticulate)

py_install("pandas")
py_install("seaborn")

import pandas as pd
import seaborn as sns

df = pd.read_csv(
    "https://raw.githubusercontent.com/selva86/datasets/master/mtcars.csv"
)

Задание 2

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

Общая информация:

print(df.shape)
print(df.info())

(32, 14)
<class 'pandas.DataFrame'> 
RangeIndex: 32 entries, 0 to 31 
Data columns (total 14 columns):  
 #   Column   Non-Null Count  Dtype 
---  ------   --------------  -----
 0   mpg      32 non-null     float64
 1   cyl      32 non-null     int64    
 2   disp     32 non-null     float64  
 3   hp       32 non-null     int64    
 4   drat     32 non-null     float64  
 5   wt       32 non-null     float64  
 6   qsec     32 non-null     float64  
 7   vs       32 non-null     int64    
 8   am       32 non-null     int64    
 9   gear     32 non-null     int64    
 10  carb     32 non-null     int64    
 11  fast     32 non-null     int64    
 12  cars     32 non-null     str      
 13  carname  32 non-null     str     
dtypes: float64(5), int64(7), str(2) 
memory usage: 3.6 KB 
None

Таким образом, имеем таблицу с 14 столбцами, 5 с вещественными значениями, 7 с целочисленными, 2 со строковыми.

Описательная статистика:

print(df.describe())

             mpg        cyl        disp  ...       gear     carb       fast
count  32.000000  32.000000   32.000000  ...  32.000000  32.0000  32.000000 
mean    4.434770   6.187500  230.721875  ...   3.687500   2.8125   0.687500 
std     0.661137   1.785922  123.938694  ...   0.737804   1.6152   0.470929 
min     3.224903   4.000000   71.100000  ...   3.000000   1.0000   0.000000 
25%     3.927432   4.000000  120.825000  ...   3.000000   2.0000   0.000000 
50%     4.381780   6.000000  196.300000  ...   4.000000   2.0000   1.000000 
75%     4.774935   8.000000  326.000000  ...   4.000000   4.0000   1.000000 
max     5.822371   8.000000  472.000000  ...   5.000000   8.0000   1.000000 

[8 rows x 12 columns]

Всего имеем 32 строки, показаны среднее, стандартное отклонение, квартили, минимумы и максимумы.

Проверка пропусков:

print(df.isnull().sum())

Как видно, пропусков нет.

Проаназирую конкретные показатели.

Исследую распредения расхода топлива mpg (сколько миль может проехать автомобиль на 1 галлоне топлива).

import matplotlib.pyplot as plt

plt.figure(figsize=(8,5))
sns.histplot(df["mpg"], bins=8, kde=True)

plt.title("Распределение расхода топлива")
plt.show()

Вывод: большинство автомобилей проезжают 4,5-4,8 миль за галлон.

Исследую зависимость расхода автомобиля от его массы.

plt.figure(figsize=(8,5))

sns.scatterplot(
    data=df,
    x="wt",
    y="mpg"
)

plt.title("Зависимость расхода топлива от массы автомобиля")
plt.show()

Вывод: с увеличением массы автомобиля ухудшается его расход.

Исследую связь кол-ва циллиндров с расходом топлива:

plt.figure(figsize=(8,5))

sns.boxplot(
    data=df,
    x="cyl",
    y="mpg"
)

plt.title("Расход топлива для автомобилей с разным числом цилиндров")
plt.show()

Вывод: чем меньше циллиндров, тем меньше расход.

Исследую корреляционную матрицу.

plt.figure(figsize=(12,10))

sns.heatmap(
    df.corr(numeric_only=True),
    annot=True,
    cmap="coolwarm",
    fmt=".2f"
)

plt.xticks(rotation=45, ha='right')
plt.yticks(rotation=0)

plt.show()

Вывод: расход (mpg) отрицательно коррелирует с массой (wt), расход (mpg) также отрицательно коррелирует с количеством лошадиных сил (hp), а лошадиных сил (hp) положительно коррелирует с количеством цилиндров (cyl).

Задание 3

Бустрэп - это метод исследования распределения статистик, основанный на многократной генерации выборок на базе имеющейся выборки.

Начну с реализации бутстрэпа на Python.

В моём датасете я разделю автомобили на две группы (с 4 цилиндрами и 8 цилиндрами), сравню расход топлива (mpg), и бутстрапом оценю доверительный интервал для разницы.

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

# Загрузка данных
df = pd.read_csv(
    "https://raw.githubusercontent.com/selva86/datasets/master/mtcars.csv"
)

# Группы для сравнения
values_a = df[df["cyl"] == 4]["mpg"].values
values_b = df[df["cyl"] == 8]["mpg"].values

# Функция построения доверительного интервала
def get_percentile_ci(bootstrap_stats, alpha):
    
    left, right = np.quantile(
        bootstrap_stats,
        [alpha / 2, 1 - alpha / 2]
    )
    
    return left, right
    
# параметры бутстрапа
B = 10000
alpha = 0.05

# разница средних
pe = np.mean(values_b) - np.mean(values_a)

bootstrap_stats = []

for _ in range(B):
    
    sample_a = np.random.choice(
        values_a,
        size=len(values_a),
        replace=True
    )
    
    sample_b = np.random.choice(
        values_b,
        size=len(values_b),
        replace=True
    )
    
    stat = np.mean(sample_b) - np.mean(sample_a)
    
    bootstrap_stats.append(stat)
    
bootstrap_stats = np.array(bootstrap_stats)

ci = get_percentile_ci(
    bootstrap_stats,
    alpha
)

has_effect = not (ci[0] < 0 < ci[1])

print(f"Среднее mpg для 4 цилиндров:, {np.mean(values_a):.2f}")
print(f"Среднее mpg для 8 цилиндров: {np.mean(values_b):.2f}")
print(f"Разница средних mpg: {pe:.2f}")
print(f"95% ДИ: ({ci[0]:.2f}; {ci[1]:.2f})")
print(f"Различия статистически значимы: {has_effect}")

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

Среднее mpg для 4 цилиндров:, 5.15
Среднее mpg для 8 цилиндров: 3.87
Разница средних mpg: -1.27
95% ДИ: (-1.57; -0.99)
Различия статистически значимы: True

Дополнительная информация:
Стандартная ошибка (SE): 0.151
Смещение (bias): 0.000
Коэффициент вариации: 0.118

plt.figure(figsize=(8,6))

sns.boxplot(
    data=df[df["cyl"].isin([4, 8])],
    x="cyl",
    y="mpg"
)

plt.title(
    "Расход топлива у автомобилей с 4 и 8 цилиндрами"
)

plt.xlabel("Число цилиндров")

plt.ylabel("MPG")

plt.show()

Вывод: В качестве исследовалась переменная mpg. Выполнено 10000 бутстрап-перевыборок с возвращением, после чего для каждой выборки вычислялась разность средних значений.

Полученная оценка разности средних составила -1.27. Построенный 95%-й доверительный интервал равен (-1.57; -0.99). Поскольку доверительный интервал не содержит нулевое значение, различия между группами являются статистически значимыми.

Стандартная ошибка оценки составила 0.151, что свидетельствует о достаточно высокой устойчивости результата. Смещение (bias) отсутствует и равно 0.000, следовательно бутстрап-оценка не демонстрирует систематического отклонения. Коэффициент вариации составил 0.118, что указывает на умеренную изменчивость полученных оценок.

Задание 4

Реализую теперь бутстрэп на R.

Линейная регрессия позволяет оценить зависимость одной переменной от другой.

Коэффициент детерминации R^2 показывает, какую долю вариации зависимой переменной объясняет модель.

Бутстрап позволяет оценить устойчивость оценки R^2, стандартную ошибку, а также доверительный интервал.

Производить буду бутстрапирование коэффициента детерминации для модели mpg~disp.

set.seed(0)

library(boot)

# функция расчета R^2
rsq_function <- function(formula, data, indices){

  d <- data[indices, ]

  fit <- lm(formula, data = d)

  return(summary(fit)$r.squared)
}

# бутстрап
reps <- boot(
  data = mtcars,
  statistic = rsq_function,
  R = 5000,
  formula = mpg ~ disp
)

reps

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = mtcars, statistic = rsq_function, R = 5000, formula = mpg ~ 
    disp)


Bootstrap Statistics :
     original      bias    std. error
t1* 0.7183433 0.003270593  0.06413799

ORDINARY NONPARAMETRIC BOOTSTRAP   

Call: 
boot(data = mtcars, statistic = rsq_function, R = 5000, formula = mpg ~
    disp)   

Bootstrap Statistics :
      original      bias    std. error
 t1* 0.7183433 0.003270593  0.06413799

График №1. Распределение бутстрап-оценок

hist(
  reps$t,
  breaks = 30,
  col = "lightblue",
  main = "Bootstrap distribution of R²",
  xlab = "R²"
)

abline(
  v = reps$t0,
  col = "red",
  lwd = 3
)

Красная линия показывает исходное значение R^2.

График №2. QQ-plot

qqnorm(reps$t)

qqline(
  reps$t,
  col = "red",
  lwd = 2
)

График показывает близость распределения бутстрап-оценок к нормальному.

Выведу доверительный интервал и дополнительные показатели.

boot.ci(
  reps,
  type = "bca"
)

bias <- mean(reps$t) - reps$t0

se <- sd(reps$t)

cv <- se / mean(reps$t)

cat("Bias =", bias, "\n")
cat("SE =", se, "\n")
cat("CV =", cv, "\n")

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 5000 bootstrap replicates  

CALL :  
boot.ci(boot.out = reps, type = "bca")  

Intervals :  
Level       BCa           
95%   ( 0.5499,  0.8195 ) 

Calculations and Intervals on Original Scale

Bias = 0.003270593  
SE = 0.06413799  
CV = 0.08888131 

Вывод

Для модели mpg∼disp выполнено 5000 бутстрап-перевыборок. Получен 95%-й BCa-доверительный интервал для коэффициента детерминации R^2: (0.5499; 0.8195). Смещение мало (Bias = 0.0033), стандартная ошибка составляет 0.0641. Это подтверждает устойчивость оценки и наличие достаточно сильной связи между объёмом двигателя и расходом топлива.

Выводы по работе

Проведён разведочный анализ датасета автомобилей и реализован метод бутстрапа в Python и R. В Python подтверждены статистически значимые различия между исследуемыми группами автомобилей. В R оценена устойчивость коэффициента детерминации регрессионной модели и построен доверительный интервал для R^2. Результаты показывают, что бутстрап позволяет надёжно оценивать статистические характеристики и их доверительные интервалы без предположения о нормальности исходных данных.