Lab4

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

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

import seaborn as sns
import matplotlib.pyplot as plt

Загружаем датасет iris

iris = sns.load_dataset('iris')

Выводим основные характеристики данных

print(iris.describe())

##        sepal_length  sepal_width  petal_length  petal_width
## 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(iris, hue='species')

plt.show()

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

library(ggplot2)

Загружаем датасет iris

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  
##                 
##                 
## 

Строим графики для анализа данных

ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width, color = Species)) + geom_point() + theme_minimal()

Задание 2: Реализация бутстрапа

Бутстрап на Python

import numpy as np
import matplotlib.pyplot as plt

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

Параметры

n = 1000
B = 10000
alpha = 0.05

Генерация данных

values_a = np.random.normal(90, 20, n)
values_b = np.random.normal(95, 15, n)

Расчет перцентилей

pe = np.quantile(values_b, 0.9)
bootstrap_values_a = np.random.choice(values_a, (B, n), replace=True)
bootstrap_metrics_a = np.quantile(bootstrap_values_a, 0.9, axis=1)
bootstrap_values_b = np.random.choice(values_b, (B, n), replace=True)
bootstrap_metrics_b = np.quantile(bootstrap_values_b, 0.9, axis=1)

Разница в перцентилях

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

Визуализация

plt.hist(bootstrap_stats, bins=50, edgecolor='k')
plt.axvline(ci[0], color='r', linestyle='dashed', linewidth=1)
plt.axvline(ci[1], color='r', linestyle='dashed', linewidth=1)
plt.title('Bootstrap Confidence Interval')
plt.xlabel('Difference in 90th Percentile')
plt.ylabel('Frequency')
plt.show()

print(f"Значение 90% квантиля изменилось на: {pe:.2f}")

## Значение 90% квантиля изменилось на: 115.76

print(f"95% доверительный интервал: ({ci[0]:.2f}, {ci[1]:.2f})")

## 95% доверительный интервал: (-2.97, 2.77)

print(f"Отличия статистически значимы: {has_effect}")

## Отличия статистически значимы: False

Бутстрап на R

library(boot)
set.seed(0)

Определим функцию модели R-квадрат

rsq_function <- function(formula, data, indices) {
  d <- data[indices, ]
  fit <- lm(formula, data = d)
  return(summary(fit)$r.squared)
}

Бутстрапируем 2000 раз

reps <- boot(data = mtcars, statistic = rsq_function, R = 2000, formula = mpg ~ disp)
print(reps)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = mtcars, statistic = rsq_function, R = 2000, formula = mpg ~ 
##     disp)
## 
## 
## Bootstrap Statistics :
##      original      bias    std. error
## t1* 0.7183433 0.003272215  0.06390141

Визуализация

plot(reps)

Персентильный доверительный интервал

boot_ci <- boot.ci(reps, type = "bca")
print(boot_ci)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = reps, type = "bca")
## 
## Intervals : 
## Level       BCa          
## 95%   ( 0.5437,  0.8149 )  
## Calculations and Intervals on Original Scale