Valiotti Analytics — построение аналитики для мобильных и digital-стартапов
    DataMarathon.ru — семидневный интенсив в области аналитики для начинающих
Позднее Ctrl + ↑

Обнаружение статистических выбросов в Python

Время чтения текста – 7 минут

Параллельно с выходом материала «Обнаружение выбросов в R» предлагаем посмотреть, как те же методы обнаружения выбросов реализовать в Python.

Данные

Для наглядности эксперимента возьмём тот же пакет данных mpg — скачать его в виде csv-таблицы можно с GitHub. Импортируем библиотеки и читаем таблицу в DataFrame:

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

df = pd.read_csv('mpg.csv')

Минимальные и максимальные значения

Тут всё просто. Выводим описание всего датасета методом describe():

df.describe()

Гистограмма

Такой график тоже можно построить в одну строку, используя внутренние средства библиотеки pandas:

df.hwy.plot(kind='hist', density=1, bins=20, stacked=False, alpha=.5, color='grey')

Box plot

В случае ящика с усами далеко идти тоже не приходится — в pandas есть метод и для этого:

_, bp = df.hwy.plot.box(return_type='both')

Получим точки с графика и выведем их в таблице, используя объект bp:

outliers = [flier.get_ydata() for flier in bp["fliers"]][0]
df[df.hwy.isin(outliers)]

Процентили

При помощи метода quantile получаем соответствующую нижнюю и верхнюю границы, а затем выводим всё, что выходит за их рамки:

lower_bound = df.hwy.quantile(q=0.025)
upper_bound = df.hwy.quantile(q=0.975)
df[(df.hwy < lower_bound) | (df.hwy > upper_bound)]

Фильтр Хэмпеля

Мы используем реализацию фильтра Хэмпеля, найденную на StackOverflow

Опишем функцию, которая заменяет на nan все значения, у которых разница с медианой больше, чем три медианных абсолютных отклонения.

def hampel(vals_orig):
    vals = vals_orig.copy()    
    difference = np.abs(vals.median()-vals)
    median_abs_deviation = difference.median()
    threshold = 3 * median_abs_deviation
    outlier_idx = difference > threshold
    vals[outlier_idx] = np.nan
    return(vals)

И применим к нашему набору данных:

hampel(df.hwy)

0      29.0
1      29.0
2      31.0
3      30.0
4      26.0
       ... 
229    28.0
230    29.0
231    26.0
232    26.0
233    26.0
Name: hwy, Length: 234, dtype: float64

В выводе нет nan-значений, а значит и выбросов фильтр Хэмпеля не обнаружил.

Тест Граббса

Автор реализации теста Граббса и теста Рознера для Python

Опишем три функции: первая находит значение критерия Граббса и максимальное значение в наборе данных, вторая — критическое значение с учётом объёма выборки и уровня значимости, а третья проверяет, является ли значение с максимальным индексом выбросом:

import numpy as np
from scipy import stats

def grubbs_stat(y):
    std_dev = np.std(y)
    avg_y = np.mean(y)
    abs_val_minus_avg = abs(y - avg_y)
    max_of_deviations = max(abs_val_minus_avg)
    max_ind = np.argmax(abs_val_minus_avg)
    Gcal = max_of_deviations / std_dev
    print(f"Grubbs Statistics Value: {Gcal}")
    return Gcal, max_ind

def calculate_critical_value(size, alpha):
    t_dist = stats.t.ppf(1 - alpha / (2 * size), size - 2)
    numerator = (size - 1) * np.sqrt(np.square(t_dist))
    denominator = np.sqrt(size) * np.sqrt(size - 2 + np.square(t_dist))
    critical_value = numerator / denominator
    print(f"Grubbs Critical Value: {critical_value}")
    return critical_value

def check_G_values(Gs, Gc, inp, max_index):
    if Gs > Gc:
        print(f"{inp[max_index]} is an outlier")
    else:
        print(f"{inp[max_index]} is not an outlier")

Заменим значение в 34 строке на 212:

df.hwy[34] = 212

И выполним три функции:

Gcritical = calculate_critical_value(len(df.hwy), 0.05)
Gstat, max_index = grubbs_stat(df.hwy)
check_G_values(Gstat, Gcritical, df.hwy, max_index)

Grubbs Critical Value: 3.652090929984981
Grubbs Statistics Value: 13.745808761040397
212 is an outlier

Тест Рознера

Для теста Рознера достаточно дописать одну функцию, которая принимает набор данных, уровень значимости и число потенциальных выбросов:

def ESD_test(input_series, alpha, max_outliers):
    for iteration in range(max_outliers):
        Gcritical = calculate_critical_value(len(input_series), alpha)
        Gstat, max_index = grubbs_stat(input_series)
        check_G_values(Gstat, Gcritical, input_series, max_index)
        input_series = np.delete(input_series, max_index)

Используя функцию на нашем наборе данных получаем, что значение 212 является выбросом, а 44 — нет:

ESD_test(np.array(df.hwy), 0.05, 3)

Grubbs Critical Value: 3.652090929984981
Grubbs Statistics Value: 13.745808761040408
212 is an outlier
Grubbs Critical Value: 3.6508358337727187
Grubbs Statistics Value: 3.455960616168714
44 is not an outlier
Grubbs Critical Value: 3.649574509044683
Grubbs Statistics Value: 3.5561478280392245
44 is not an outlier

Обнаружение статистических выбросов в R

Время чтения текста – 27 минут

Этот материал — перевод статьи «Outliers detection in R». А ещё у нас есть материал про обнаружение выбросов в Python.

Выбросы — значения или наблюдения, отклоняющиеся от других данных. Всегда нужно сравнивать наблюдение с другими значениями, полученными тем же способом, прежде чем называть их выбросами. Действительно, человек с ростом 200 см, скорее всего, будет считаться отклонением по сравнению с остальным населением, но этот же человек не будет считаться статистическим выбросом, если мы измерим рост баскетболистов.

Выбросы могут быть вызваны изменчивостью, присущей наблюдаемому явлению. Например, при сборе данных о заработной плате часто возникают выбросы, поскольку некоторые люди зарабатывают гораздо больше остальных. Выбросы также могут возникать из-за экспериментальной ошибки, ошибки измерения или кодирования. Например, вес человека 786 кг явно является ошибкой при кодировании веса объекта. Её или его вес, скорее всего, составляет 78,6 кг или 7,86 кг в зависимости от того, был измерен вес взрослого человека или ребёнка.

По этой причине иногда имеет смысл формально выделять два класса выбросов: экстремальные значения и ошибки. Экстремальные значения интереснее, потому что они возможны, но маловероятны.

В этой статье я представлю несколько подходов к обнаружению выбросов в R от простых методов, таких как описательная статистика (включая минимальные, максимальные значения, гистограмму, прямоугольную диаграмму и процентили), до более формальных методов, таких как фильтр Хэмпеля, тесты Граббса, Диксона и Рознера.

Не существует строгого и однозначного правила насчет того, следует ли удалять выбросы из набора данных перед проведением статистического анализа. Довольно часто переменные значения, вызванные экспериментальной ошибкой или ошибкой измерения удаляются или заменяются на новые значения. Некоторые статистические тесты требуют их отсутствия, чтобы сделать обоснованные выводы, но удаление выбросов рекомендуется не во всех случаях и должно выполняться с осторожностью.

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

  1. Область / контекст вашего анализа и вопрос исследования. В некоторых областях обычно удаляют посторонние значения, поскольку они часто возникают из-за сбоев в процессе. В других областях отклонения сохраняются, потому что они содержат ценную информацию. Также бывает, что анализ выполняется дважды, один раз с посторонними значениями и один раз без них, чтобы оценить их влияние на результаты. Если результаты резко изменятся из-за некоторых определяющих значений, это должно предостеречь исследователя от чрезмерно амбициозных утверждений.
  2. Устойчивость тестов. Например, наклон простой линейной регрессии может значительно варьироваться даже с одним выбросом, тогда как непараметрические тесты, такие как тест Уилкоксона, обычно устойчивы к ним.
  3. Дальность выбросов от других наблюдений. Некоторые наблюдения, рассматриваемые как выбросы, на самом деле не являются экстремальными значениями по сравнению со всеми другими наблюдениями, в то время как другие потенциальные выбросы могут быть действительно отстающими от остальных наблюдений.

Мы будем использовать набор данных mpg из библиотеки ggplot2, чтобы проиллюстрировать различные подходы к обнаружению выбросов в R, и в частности, мы сосредоточимся на работе с переменной hwy (пробег в милях на галлон израсходованного топлива).

Минимальные и максимальные значения

Первое, что необходимо для обнаружения выбросов в R — начать с описательной статистики, и, в частности, с минимальных и максимальных значений.

В R это легко сделать с помощью функции summary():

dat <- ggplot2::mpg
summary(dat$hwy)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   12.00   18.00   24.00   23.44   27.00   44.00

Минимум и максимум — первое и последнее значения в выходных данных выше. В качестве альтернативы, их также можно вычислить с помощью функций min() и max():

min(dat$hwy)

## [1] 12

max(dat$hwy)

## [1] 44

Явная ошибка кодирования, такая как, например, человеческий вес в 786 кг уже будет легко обнаружена с помощью этой простой техники.

Гистограмма

Другой базовый способ обнаружения выбросов — построение гистограммы данных.

При помощи внутренних инструментов R:

hist(dat$hwy,
  xlab = "hwy",
  main = "Histogram of hwy",
  breaks = sqrt(nrow(dat))
) # set number of bins

При помощи ggplot2:

library(ggplot2)

ggplot(dat) +
  aes(x = hwy) +
  geom_histogram(bins = 30L, fill = "#0c4c8a") +
  theme_minimal()

Пара полосок справа в отрыве от основного графика — значения, которые больше остальных.

#Box plot
Помимо гистограмм, box plot (ящик с усами) также полезен для обнаружения потенциальных выбросов.

Используя R:

boxplot(dat$hwy,
  ylab = "hwy"
)

или используя ggplot2:

ggplot(dat) +
  aes(x = "", y = hwy) +
  geom_boxplot(fill = "#0c4c8a") +
  theme_minimal()

Box plot помогает визуализировать количественную переменную, отображая пять общих сводных данных (минимальное значение, среднее значение, первый и третий квартили и максимальное значение) и любое значение, которое было классифицировано как предполагаемый выброс с использованием критерия межквартильного размаха (IQR). Критерий межквартильного размаха означает, что все единицы значения больше q₀,₇₅+ 1.5 ⋅ IQR или меньше q₀,₂₅ — 1,5⋅ IQR рассматриваются R, как потенциальные выбросы. Другими словами, все наблюдения за пределами следующего интервала будут рассматриваться как потенциальные выбросы:

I = [q₀,₂₅ — 1.5 * IQR; q₀,₇₅ + 1.5 * IQR]

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

Даже если наблюдение рассматривается как потенциальный выброс по критерию IQR, это не означает, что его следует удалять. Удаление или сохранение выброса зависит от контекста вашего анализа, от того, являются ли тесты, которые вы собираетесь проводить с наборами данных, устойчивыми к выбросам или нет, и насколько далеки выбросы от других наблюдений.

Также возможно извлечь потенциальные выбросы на основе критерия IQR благодаря функции boxplot.stats()$out:

boxplot.stats(dat$hwy)$out

## [1] 44 44 41

Как видите, на самом деле есть 3 точки, которые считаются потенциальными выбросами: две со значением 44 и одна со значением 41.

Благодаря функции which() можно извлечь номер строки, соответствующий этим посторонним значениям:

out <- boxplot.stats(dat$hwy)$out
out_ind <- which(dat$hwy %in% c(out))
out_ind

## [1] 213 222 223

Имея эту информацию, вы теперь можете легко вернуться к определенным строкам в наборе данных, чтобы проверить их, или напечатать все переменные для этих выбросов:

dat[out_ind, ]

## # A tibble: 3 x 11
##   manufacturer model   displ  year   cyl trans   drv     cty   hwy fl    class  
##   <chr>        <chr>   <dbl> <int> <int> <chr>   <chr> <int> <int> <chr> <chr>  
## 1 volkswagen   jetta     1.9  1999     4 manual… f        33    44 d     compact
## 2 volkswagen   new be…   1.9  1999     4 manual… f        35    44 d     subcom…
## 3 volkswagen   new be…   1.9  1999     4 auto(l… f        29    41 d     subcom…

Ещё можно напечатать выбросы прямо на диаграмме размаха с помощью функции mtext():

boxplot(dat$hwy,
  ylab = "hwy",
  main = "Boxplot of highway miles per gallon"
)
mtext(paste("Outliers: ", paste(out, collapse = ", ")))

Процентили

Этот метод обнаружения посторонних значений основан на процентилях. При использовании метода процентилей все наблюдения, выходящие за пределы интервала, образованного 2,5 и 97,5 процентилями будут рассматриваться как потенциальные выбросы. Другие процентили, такие как 1 и 99 или 5 и 95 процентили, тоже могут быть рассмотрены для построения интервала.

Значения нижнего и верхнего процентилей можно вычислить с помощью функции quantile():

lower_bound <- quantile(dat$hwy, 0.025)
lower_bound

## 2.5% 
##   14

upper_bound <- quantile(dat$hwy, 0.975)
upper_bound

##  97.5% 
## 35.175

В соответствии с этим методом, все наблюдения ниже 14 и выше 35,175 будут рассматриваться как потенциальные выбросы. Номера рядов наблюдений за пределами интервала затем могут быть извлечены с помощью функции which():

outlier_ind <- which(dat$hwy < lower_bound | dat$hwy > upper_bound)
outlier_ind

##  [1]  55  60  66  70 106 107 127 197 213 222 223

Можно вывести значение пробега в милях на галлон израсходованного топлива для таких значений:

dat[outlier_ind, "hwy"]


## # A tibble: 11 x 1
##      hwy
##    <int>
##  1    12
##  2    12
##  3    12
##  4    12
##  5    36
##  6    36
##  7    12
##  8    37
##  9    44
## 10    44
## 11    41

В качестве альтернативы можно вывести все переменные для этих выбросов:

dat[outlier_ind, ]

## # A tibble: 11 x 11
##    manufacturer model    displ  year   cyl trans  drv     cty   hwy fl    class 
##    <chr>        <chr>    <dbl> <int> <int> <chr>  <chr> <int> <int> <chr> <chr> 
##  1 dodge        dakota …   4.7  2008     8 auto(… 4         9    12 e     pickup
##  2 dodge        durango…   4.7  2008     8 auto(… 4         9    12 e     suv   
##  3 dodge        ram 150…   4.7  2008     8 auto(… 4         9    12 e     pickup
##  4 dodge        ram 150…   4.7  2008     8 manua… 4         9    12 e     pickup
##  5 honda        civic      1.8  2008     4 auto(… f        25    36 r     subco…
##  6 honda        civic      1.8  2008     4 auto(… f        24    36 c     subco…
##  7 jeep         grand c…   4.7  2008     8 auto(… 4         9    12 e     suv   
##  8 toyota       corolla    1.8  2008     4 manua… f        28    37 r     compa…
##  9 volkswagen   jetta      1.9  1999     4 manua… f        33    44 d     compa…
## 10 volkswagen   new bee…   1.9  1999     4 manua… f        35    44 d     subco…
## 11 volkswagen   new bee…   1.9  1999     4 auto(… f        29    41 d     subco…

Согласно методу процентилей, существует 11 потенциальных выбросов. Чтобы уменьшить это число, вы можете установить процентили от 1 до 99:

lower_bound <- quantile(dat$hwy, 0.01)
upper_bound <- quantile(dat$hwy, 0.99)

outlier_ind <- which(dat$hwy < lower_bound | dat$hwy > upper_bound)

dat[outlier_ind, ]

## # A tibble: 3 x 11
##   manufacturer model   displ  year   cyl trans   drv     cty   hwy fl    class  
##   <chr>        <chr>   <dbl> <int> <int> <chr>   <chr> <int> <int> <chr> <chr>  
## 1 volkswagen   jetta     1.9  1999     4 manual… f        33    44 d     compact
## 2 volkswagen   new be…   1.9  1999     4 manual… f        35    44 d     subcom…
## 3 volkswagen   new be…   1.9  1999     4 auto(l… f        29    41 d     subcom…

Установка процентилей на 1 и 99 дает те же потенциальные выбросы, что и для критерия IQR.

Фильтр Хэмпеля

Другой метод, известный как фильтр Хэмпеля, заключается в том, чтобы рассматривать как выбросы значения вне интервала, которые формируются медианным значением плюс-минус 3 медианы абсолютных отклонений (MAD):

I = [median - 3 * MAD; median + 3 * MAD]

в которых MAD — это медианное абсолютное отклонение и определяется как медиана абсолютных отклонений от медианы данных:

Для этого метода мы сначала устанавливаем пределы интервала с помощью функций median() и mad():

lower_bound <- median(dat$hwy) - 3 * mad(dat$hwy, constant=1)
lower_bound

## [1] 9

upper_bound <- median(dat$hwy) + 3 * mad(dat$hwy, constant=1)
upper_bound

## [1] 39

Все наблюдения меньше 9 и больше 39 будут рассматриваться как потенциальные выбросы. Номера строк наблюдений за пределами интервала затем могут быть извлечены с помощью функции which():

outlier_ind <- which(dat$hwy < lower_bound | dat$hwy > upper_bound)
outlier_ind

## 213 222 223

Согласно фильтру Хэмпеля, для переменной hwy есть 3 потенциальных выброса.

Статистические тесты

В этом разделе мы представим еще 3 формальных метода обнаружения отклонений:

  1. Тест Граббса (Grubbs’s test)
  2. Тест Диксона (Dixon’s test)
  3. Тест Рознера (Rosner’s test)

Эти статистические тесты являются частью формальных методов обнаружения выбросов, поскольку все они включают вычисление с помощью тестовой статистики, которая сравнивается с табличными критическими значениями.

Обратите внимание, что эти тесты подходят только тогда, когда данные распределены нормально. Таким образом, предположение о соответствии нормальности должно быть проверено перед применением этих тестов для выбросов (Как проверить предположение о соответствии нормальному распределению в R).

Тест Граббса (Grubbs’s test)

Тест Граббса позволяет определить, является ли наибольшее или наименьшее значение в наборе данных выбросом. Он обнаруживает по одному выбросу за раз (максимальное или минимальное значение), поэтому нулевая и альтернативная гипотезы проверки максимального значения выглядит так:

  • H₀: Наивысшее значение не является выбросом
  • H₁: Наивысшее значение является выбросом

А минимального — так:

  • H₀: Наименьшее значение не является выбросом
  • H₁: Наименьшее значение является выбросом

Как и в любом статистическом тесте, если значение P меньше порогового уровня статистической значимости (обычно α = 0.05), то нулевая гипотеза отвергается, и мы приходим к выводу, что наименьшее/наибольшее значение является отклонением. Напротив, если значение P больше или равно пороговому уровню значимости, нулевая гипотеза не отвергается, и мы делаем вывод, что на основе данных о том, что наименьшее / наибольшее значение не является выбросом. Обратите внимание на то, что тест Граббса не подходит для выборки объемом 6 или меньше (n <= 6).

Чтобы выполнить тест Граббса в R, используем функцию grubbs.test() из пакетов outliers:

# install.packages("outliers")
library(outliers)
test <- grubbs.test(dat$hwy)
test 

## 
##  Grubbs test for one outlier
## 
## data:  dat$hwy
## G = 3.45274, U = 0.94862, p-value = 0.05555
## alternative hypothesis: highest value 44 is an outlier

Значение P равняется 0,056. На уровне значимости 5% мы не отвергаем гипотезу о том, что наибольшее значение 44 не является выбросом.

По умолчанию тест выполняется на наибольшем значении (как показано в выходных данных R: alternative hypothesis: highest value). Если вы хотите провести тест для наименьшего значения, просто добавьте аргумент opposite = TRUE в функцию grubbs.test():

test <- grubbs.test(dat$hwy, opposite = TRUE)
test

## 
##  Grubbs test for one outlier
## 
## data:  dat$hwy
## G = 1.92122, U = 0.98409, p-value = 1
## alternative hypothesis: lowest value 12 is an outlier

Вывод указывает на то, что тест теперь выполняется при наименьшем значении

Значение P равно 1. На уровне значимости 5% мы не отвергаем гипотезу о том, что наименьшее значение 12 не является выбросом.

Для иллюстрации этого заменим наблюдения более экстремальным значением и выполним тест Граббса для нового набора данных. Давайте заменим 34-ую строку со значением 212:

dat[34, "hwy"] <- 212

Применяем тест Граббса, чтобы проверить, является ли наибольшее значение выбросом:

test <- grubbs.test(dat$hwy)
test

## 
##  Grubbs test for one outlier
## 
## data:  dat$hwy
## G = 13.72240, U = 0.18836, p-value < 2.2e-16
## alternative hypothesis: highest value 212 is an outlier

Значение p < 0,001. На уровне значимости 5% мы делаем вывод, что наивысшее значение 212 является выбросом.

Тест Диксона (Dixon’s test)

Подобно тесту Граббса, тест Диксона используется для того, чтобы проверить, является ли самое высокое или самое низкое значение выбросом. Таким образом, если под сомнением находятся более одного выброса, тест необходимо проводить индивидуально для этих предполагаемых значений.

Обратите внимание на то, что тест Диксона наиболее полезен для выборки небольшого объема (обычно когда n <= 25).

Чтобы выполнить тест Диксона в R, мы используем функцию dixon.test() из пакета outliers. Однако мы ограничиваем наш набор данных 20 первыми наблюдениями, поскольку тест Диксона может быть выполнен только на небольшом размере выборки:

subdat <- dat[1:20, ]
test <- dixon.test(subdat$hwy)
test

## 
##  Dixon test for outliers
## 
## data:  subdat$hwy
## Q = 0.57143, p-value = 0.006508
## alternative hypothesis: lowest value 15 is an outlier

Результаты показывают, что самое наименьшее значение 15 является выбросом (p-значение = 0,007).

Чтобы проверить максимальное значение, просто добавьте аргумент opposite = TRUE к функции dixon.test():

test <- dixon.test(subdat$hwy,
  opposite = TRUE
)
test

## 
##  Dixon test for outliers
## 
## data:  subdat$hwy
## Q = 0.25, p-value = 0.8582
## alternative hypothesis: highest value 31 is an outlier

Результаты показывают, что максимальное значение 31 не является выбросом (p-значение = 0,858).

Рекомендуется всегда сверять результаты статистического теста на выбросы с диаграммой, чтобы убедиться, что мы проверили все потенциальные выбросы:

out <- boxplot.stats(subdat$hwy)$out
boxplot(subdat$hwy,
  ylab = "hwy"
)
mtext(paste("Outliers: ", paste(out, collapse = ", ")))

По box plot заметно, что мы можем применить тест Диксона к значению 20 в дополнение к значению 15, выполненному ранее. Это можно сделать, найдя номер строки минимального значения, исключив этот номер строки из набора данных, а затем применив тест Диксона к этому новому набору данных:

# find and exclude lowest value
remove_ind <- which.min(subdat$hwy)
subsubdat <- subdat[-remove_ind, ]

# Dixon test on dataset without the minimum
test <- dixon.test(subsubdat$hwy)
test

## 
##  Dixon test for outliers
## 
## data:  subsubdat$hwy
## Q = 0.44444, p-value = 0.1297
## alternative hypothesis: lowest value 20 is an outlier

Результаты показывают, что второе наименьшее значение 20 не является выбросом (p-значение = 0,13).

Тест Рознера (Rosner’s test)

  1. Тест Рознера на выбросы имеет следующие преимущества:
  2. Он используется для одновременного обнаружения нескольких выбросов (в отличие от теста Граббса и Диксона, которые должны выполняться итеративно для выявления нескольких выбросов)
  3. Он разработан чтобы избежать проблемы, когда выброс, близкий по значению к другому выбросу, может остаться незамеченным.

Обратите внимание, что в отличие от теста Диксона, тест Рознера подходит к большому объему выборки (n≥20). Поэтому мы снова используем исходный набор данных dat, который включает 234 наблюдения.

Для выполнения теста Рознера мы используем функцию rosnerTest() из пакета EnvStats. Для этой функции требуется как минимум 2 аргумента: данные и количество предполагаемых выбросов k.

library(EnvStats)
test <- rosnerTest(dat$hwy,
  k = 3
)
test

## $distribution
## [1] "Normal"
## 
## $statistic
##       R.1       R.2       R.3 
## 13.722399  3.459098  3.559936 
## 
## $sample.size
## [1] 234
## 
## $parameters
## k 
## 3 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 
## 3.652091 3.650836 3.649575 
## 
## $n.outliers
## [1] 1
## 
## $alternative
## [1] "Up to 3 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1]  29  29  31  30  26  26  27  26  25  28  27  25  25  25  25  24  25  23
##  [19]  20  15  20  17  17  26  23  26  25  24  19  14  15  17  27 212  26  29
##  [37]  26  24  24  22  22  24  24  17  22  21  23  23  19  18  17  17  19  19
##  [55]  12  17  15  17  17  12  17  16  18  15  16  12  17  17  16  12  15  16
##  [73]  17  15  17  17  18  17  19  17  19  19  17  17  17  16  16  17  15  17
##  [91]  26  25  26  24  21  22  23  22  20  33  32  32  29  32  34  36  36  29
## [109]  26  27  30  31  26  26  28  26  29  28  27  24  24  24  22  19  20  17
## [127]  12  19  18  14  15  18  18  15  17  16  18  17  19  19  17  29  27  31
## [145]  32  27  26  26  25  25  17  17  20  18  26  26  27  28  25  25  24  27
## [163]  25  26  23  26  26  26  26  25  27  25  27  20  20  19  17  20  17  29
## [181]  27  31  31  26  26  28  27  29  31  31  26  26  27  30  33  35  37  35
## [199]  15  18  20  20  22  17  19  18  20  29  26  29  29  24  44  29  26  29
## [217]  29  29  29  23  24  44  41  29  26  28  29  29  29  28  29  26  26  26
## 
## $data.name
## [1] "dat$hwy"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i      SD.i Value Obs.Num     R.i+1 lambda.i+1 Outlier
## 1 0 24.21795 13.684345   212      34 13.722399   3.652091    TRUE
## 2 1 23.41202  5.951835    44     213  3.459098   3.650836   FALSE
## 3 2 23.32328  5.808172    44     222  3.559936   3.649575   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"

Результаты представлены в таблице $all.stats:

test$all.stats

##   i   Mean.i      SD.i Value Obs.Num     R.i+1 lambda.i+1 Outlier
## 1 0 24.21795 13.684345   212      34 13.722399   3.652091    TRUE
## 2 1 23.41202  5.951835    44     213  3.459098   3.650836   FALSE
## 3 2 23.32328  5.808172    44     222  3.559936   3.649575   FALSE

Основываясь на тесте Рознера, мы видим, что существует только один выброс (см. Столбец Outlier), и что это наблюдение 34 (см. Obs.Num) со значением 212 (см. Value).

Итоги

Обратите внимание, что некоторые преобразования могут «естественным образом» устранить выбросы. Например, если взять натуральный логарифм или квадратный корень из значения, отклонение станет меньше. Я надеюсь, статья помогла вам обнаружить выбросы в R с помощью нескольких методов описательной статистики (включая минимум, максимум, гистограмму, диаграмму размаха и процентили) или благодаря более формальным методам обнаружения выбросов (включая фильтр Хампеля, тест Граббса, Диксона и Рознера). Следующим этапом проверьте эти значения, и если они действительно являются выбросами — решите, как с ними поступить (сохранить, удалить или изменить), прежде чем проводить анализ.

Обзор дашборда Yandex DataLens

Время чтения текста – 2 минуты

Два года назад Яндекс выпустил собственный инструмент для визуализации данных — Yandex DataLens, работающий на базе Yandex Cloud. В блоге уже выходил обзор инструмента — но тогда сервис был на стадии Preview, и за два года функционал инструмента расширили. Сервис тарифицируемый и без привязки платёжного аккаунта поработать в нём не получится, но помимо платного тарифа есть и бесплатный.

Подробнее о тарифах Yandex DataLens можно почитать в документации

В сегодняшнем обзоре BI-систем мы посмотрим, как зарегистрировать аккаунт в DataLens, подключить датасет и создать дополнительные таблицы на основе SQL-запросов, построить визуализации, связать их с фильтрами и добавить на дашборд согласно макету, а затем опубликовать результат.

Внутри команды мы оценили дашборд в DataLens и получили следующие средние оценки (1 — худшая оценка, 10 — лучшая):

Отвечает ли заданным вопросам — 7,0
Порог входа в инструмент — 8,0
Функциональность инструмента — 7,0
Удобство пользования — 8,3
Соответствие результата макету — 7,5
Визуальная составляющая — 8,5
Итог: дашборд получает 8 баллов из 10. Посмотрите на полученный результат.

Бот для преобразования данных из Coinkeeper

Время чтения текста – 6 минут

Coinkeeper — кроссплатформенное приложение для учёта финансов. Внутри можно выпустить виртуальную банковскую карту Visa с бесплатным годовым обслуживанием, которая будет присылать уведомления, если вы тратите больше, чем запланировали. Помимо уведомлений, приложение ведёт историю трат и позволяет выгрузить сводный отчёт в формате csv. Данные, которое выгружает приложение ещё не готовы к анализу и выглядят так:

Азат Шарипов сделал скрипт обработки данных в пригодный для Tableau вид и подготовил Tableau Public книгу, а Рома Бунин в рамках своего проекта «Переверстка» переработал дашборд.

Мы решили тоже поучаствовать, и с нашей стороны Елизавета Мазурова сделала чат-бота.

Чат-бот крутой! Помимо того, что он может как и прежде отдавать обратно .csv-файл, он позволяет автоматизировать рутину по обновлению отчета через Google-таблицы. Как, наверное, многие помнят, Tableau Public может работать на гугл-таблицах или csv файлах, но не разрешает подключение к данным. Бот умный: он создаст за вас гугл-таблицу и когда вы повторно отправите ему новый файл обновит ее.

Использование бота

Перейдите в диалог с ботом и введите команду /start — в ответе бот расскажет немного о себе. Для продолжения работы нажмите на кнопку «Начать».

Сразу после можно отправить csv-файл, выгруженный из Coinkeeper:

Выберите тип файла — csv или таблицу в Google Spreadsheets.

В случае выбора csv-файла бот пришлёт его:

А в случае ссылки в первый раз нужно будет пройти небольшую регистрацию — указать почту и наименование для файла.

Затем бот пришлёт ссылку на файл:

Скрипт преобразовал данные, и таблицу можно указать в качестве источника данных в Tableau. А благодаря тому, что в случае загрузки нового файла создаётся не новая таблица, а обновляется старая, отчёт в Tableau тоже обновится. В результате открывается возможность еженедельно присылать боту новую таблицу и сразу переходить в обновлённый отчёт.

Python и тексты нового альбома Земфиры: анализируем суть песен

Время чтения текста – 18 минут

Неделю назад вышёл первый за 8 лет студийный альбом Земфиры «Бордерлайн». К работе помимо рок-певицы приложили руку разные люди, в том числе и её родственники — рифф для песни «таблетки» написал её племянник из Лондона. Альбом получился разнообразным: например, песня «остин» посвящена главному персонажу игры Homescapes российской студии Playrix (кстати, посмотрите свежие Бизнес-секреты с братьями Бухманами, там они тоже про это рассказывают) — Земфире нравится игра, и для трека она связалась со студией. А сингл «крым» был написан в качестве саундтрека к новой картине соратницы Земфиры — Ренаты Литвиновой.

Послушать альбом в Apple Music / Яндекс.Музыке / Spotify

Тем не менее, дух всего альбома довольно мрачен — в песнях часто повторяются слова «боль», «ад», «бесишь» и прочие по смыслу. Мы решили провести разведочный анализ нового альбома, а затем при помощи модели Word2Vec и косинусной меры посмотреть на семантическую близость песен между собой и вычислить общее настроение альбома.

Для тех, кому скучно читать про подготовку данных и шаги анализа можно перейти сразу к результатам.

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

Для начала работы напишем скрипт обработки данных. Цель скрипта — из множества текстовых файлов, в каждом из которых лежит по песне, собрать единую csv-таблицу. При этом текст треков очищаем от знаков пунктуации и ненужных слов.

import pandas as pd
import re
import string
import pymorphy2
from nltk.corpus import stopwords

Создаём морфологический анализатор и расширяем список всего, что нужно отбросить:

morph = pymorphy2.MorphAnalyzer()
stopwords_list = stopwords.words('russian')
stopwords_list.extend(['куплет', 'это', 'я', 'мы', 'ты', 'припев', 'аутро', 'предприпев', 'lyrics', '1', '2', '3', 'то'])
string.punctuation += '—'

Названия песен приведены на английском — создадим словарь для перевода на русский и словарь, из которого позднее сделаем таблицу:

result_dict = dict()

songs_dict = {
    'snow':'снег идёт',
    'crimea':'крым',
    'mother':'мама',
    'ostin':'остин',
    'abuse':'абьюз',
    'wait_for_me':'жди меня',
    'tom':'том',
    'come_on':'камон',
    'coat':'пальто',
    'this_summer':'этим летом',
    'ok':'ок',
    'pills':'таблетки'
}

Опишем несколько функций. Первая читает целиком песню из файла и удаляет переносы строки, вторая очищает текст от ненужных символов и слов, а третья при помощи морфологического анализатора pymorphy2 приводит слова к нормальной форме. Модуль pymorphy2 не всегда хорошо справляется с неоднозначностью — для слов «ад» и «рай» потребуется дополнительная обработка.

def read_song(filename):
    f = open(f'{filename}.txt', 'r').read()
    f = f.replace('\n', ' ')
    return f

def clean_string(text):
    text = re.split(' |:|\.|\(|\)|,|"|;|/|\n|\t|-|\?|\[|\]|!', text)
    text = ' '.join([word for word in text if word not in string.punctuation])
    text = text.lower()
    text = ' '.join([word for word in text.split() if word not in stopwords_list])
    return text

def string_to_normal_form(string):
    string_lst = string.split()
    for i in range(len(string_lst)):
        string_lst[i] = morph.parse(string_lst[i])[0].normal_form
        if (string_lst[i] == 'аду'):
            string_lst[i] = 'ад'
        if (string_lst[i] == 'рая'):
            string_lst[i] = 'рай'
    string = ' '.join(string_lst)
    return string

Проходим по каждой песне и читаем файл с соответствующим названием:

name_list = []
text_list = []
for song, name in songs_dict.items():
    text = string_to_normal_form(clean_string(read_song(song)))
    name_list.append(name)
    text_list.append(text)

Затем объединяем всё в DataFrame и сохраняем в виде csv-файла.

df = pd.DataFrame()
df['name'] = name_list
df['text'] = text_list
df['time'] = [290, 220, 187, 270, 330, 196, 207, 188, 269, 189, 245, 244]
df.to_csv('borderline.csv', index=False)

Результат:

Облако слов по всему альбому

Начнём анализ с построения облака слов — оно отобразит, какие слова чаще всего встречаются в песнях. Импортируем нужные библиотеки, читаем csv-файл и устанавливаем конфигурации:

import nltk
from wordcloud import WordCloud
import pandas as pd
import matplotlib.pyplot as plt
from nltk import word_tokenize, ngrams

%matplotlib inline
nltk.download('punkt')
df = pd.read_csv('borderline.csv')

Теперь создаём новую фигуру, устанавливаем параметры оформления и при помощи библиотеки wordcloud отображаем слова с размером прямо пропорциональным частоте упоминания слова. Над каждым графиком дополнительно указываем название песни.

fig = plt.figure()
fig.patch.set_facecolor('white')
plt.subplots_adjust(wspace=0.3, hspace=0.2)
i = 1
for name, text in zip(df.name, df.text):
    tokens = word_tokenize(text)
    text_raw = " ".join(tokens)
    wordcloud = WordCloud(colormap='PuBu', background_color='white', contour_width=10).generate(text_raw)
    plt.subplot(4, 3, i, label=name,frame_on=True)
    plt.tick_params(labelsize=10)
    plt.imshow(wordcloud)
    plt.axis("off")
    plt.title(name,fontdict={'fontsize':7,'color':'grey'},y=0.93)
    plt.tick_params(labelsize=10)
    i += 1

EDA текстов альбома

Теперь проанализируем тексты песен — импортируем библиотеки для работы с данными и визуализации:

import plotly.graph_objects as go
import plotly.figure_factory as ff
from scipy import spatial
import collections
import pymorphy2
import gensim

morph = pymorphy2.MorphAnalyzer()

Сначала посчитаем число слов в каждой песне, число уникальных слов и процентное соотношение:

songs = []
total = []
uniq = []
percent = []

for song, text in zip(df.name, df.text):
    songs.append(song)
    total.append(len(text.split()))
    uniq.append(len(set(text.split())))
    percent.append(round(len(set(text.split())) / len(text.split()), 2) * 100)

А теперь составим из этого DataFrame и дополнительно посчитаем число слов в минуту для каждой песни:

df_words = pd.DataFrame()
df_words['song'] = songs
df_words['total words'] = total
df_words['uniq words'] = uniq
df_words['percent'] = percent
df_words['time'] = df['time']
df_words['words per minute'] = round(total / (df['time'] // 60))
df_words = df_words[::-1]

Данные хорошо бы визуализировать — построим две столбиковые диаграммы: одну для числа слов в песне, а другую для числа слов в минуту.

colors_1 = ['rgba(101,181,205,255)'] * 12
colors_2 = ['rgba(62,142,231,255)'] * 12

fig = go.Figure(data=[
    go.Bar(name='📝 Всего слов',
           text=df_words['total words'],
           textposition='auto',
           x=df_words.song,
           y=df_words['total words'],
           marker_color=colors_1,
           marker=dict(line=dict(width=0)),),
    go.Bar(name='🌀 Уникальных слов',
           text=df_words['uniq words'].astype(str) + '<br>'+ df_words.percent.astype(int).astype(str) + '%' ,
           textposition='inside',
           x=df_words.song,
           y=df_words['uniq words'],
           textfont_color='white',
           marker_color=colors_2,
           marker=dict(line=dict(width=0)),),
])

fig.update_layout(barmode='group')

fig.update_layout(
    title = 
        {'text':'<b>Соотношение числа уникальных слов к общему количеству</b><br><span style="color:#666666"></span>'},
    showlegend = True,
    height=650,
    font={
        'family':'Open Sans, light',
        'color':'black',
        'size':14
    },
    plot_bgcolor='rgba(0,0,0,0)',
)
fig.update_layout(legend=dict(
    yanchor="top",
    xanchor="right",
))

fig.show()
colors_1 = ['rgba(101,181,205,255)'] * 12
colors_2 = ['rgba(238,85,59,255)'] * 12

fig = go.Figure(data=[
    go.Bar(name='⏱️ Длина трека, мин.',
           text=round(df_words['time'] / 60, 1),
           textposition='auto',
           x=df_words.song,
           y=-df_words['time'] // 60,
           marker_color=colors_1,
           marker=dict(line=dict(width=0)),
          ),
    go.Bar(name='🔄 Слов в минуту',
           text=df_words['words per minute'],
           textposition='auto',
           x=df_words.song,
           y=df_words['words per minute'],
           marker_color=colors_2,
           textfont_color='white',
           marker=dict(line=dict(width=0)),
          ),
])

fig.update_layout(barmode='overlay')

fig.update_layout(
    title = 
        {'text':'<b>Длина трека и число слов в минуту</b><br><span style="color:#666666"></span>'},
    showlegend = True,
    height=650,
    font={
        'family':'Open Sans, light',
        'color':'black',
        'size':14
    },
    plot_bgcolor='rgba(0,0,0,0)'
)


fig.show()

Работа с Word2Vec моделью

При помощи модуля gensim загружаем модель, указывая на бинарный файл:

model = gensim.models.KeyedVectors.load_word2vec_format('model.bin', binary=True)

Для материала мы использовали готовую обученную на Национальном Корпусе Русского Языка модель от сообщества RusVectōrēs

Модель Word2Vec основана на нейронных сетях и позволяет представлять слова в виде векторов, учитывая семантическую составляющую. Это означает, что если мы возьмём два слова — например, «мама» и «папа», представим их в виде двух векторов и посчитаем косинус, значения будет близко к 1. Аналогично, у двух слов, не имеющих ничего общего по смыслу косинусная мера близка к 0.

Опишем функцию get_vector: она будет принимать список слов, распознавать для каждого часть речи, а затем получать и суммировать вектора — так мы сможем находить вектора не для одного слова, а для целых предложений и текстов.

def get_vector(word_list):
    vector = 0
    for word in word_list:
        pos = morph.parse(word)[0].tag.POS
        if pos == 'INFN':
            pos = 'VERB'
        if pos in ['ADJF', 'PRCL', 'ADVB', 'NPRO']:
            pos = 'NOUN'
        if word and pos:
            try:
                word_pos = word + '_' + pos
                this_vector = model.word_vec(word_pos)
                vector += this_vector
            except KeyError:
                continue
    return vector

Для каждой песни находим вектор и собираем соответствующий столбец в DataFrame:

vec_list = []
for word in df['text']:
    vec_list.append(get_vector(word.split()))
df['vector'] = vec_list

Теперь сравним вектора между собой, посчитав их косинусную близость. Те песни, у которых косинусная метрика выше 0,5 запомним отдельно — так мы получим самые близкие пары песен. Данные о сравнении векторов запишем в двумерный список result.

similar = dict()
result = []
for song_1, vector_1 in zip(df.name, df.vector):
    sub_list = []
    for song_2, vector_2 in zip(df.name.iloc[::-1], df.vector.iloc[::-1]):
        res = 1 - spatial.distance.cosine(vector_1, vector_2)
        if res > 0.5 and song_1 != song_2 and (song_1 + ' / ' + song_2 not in similar.keys() and song_2 + ' / ' + song_1 not in similar.keys()):
            similar[song_1 + ' / ' + song_2] = round(res, 2)
        sub_list.append(round(res, 2))
    result.append(sub_list)

Самые похожие треки соберём в отдельный DataFrame:

df_top_sim = pd.DataFrame()
df_top_sim['name'] = list(similar.keys())
df_top_sim['value'] = list(similar.values())
df_top_sim.sort_values(by='value', ascending=False)

И построим такой же bar chart:

colors = ['rgba(101,181,205,255)'] * 5

fig = go.Figure([go.Bar(x=df_top_sim['name'],
                        y=df_top_sim['value'],
                        marker_color=colors,
                        width=[0.4,0.4,0.4,0.4,0.4],
                        text=df_top_sim['value'],
                        textfont_color='white',
                        textposition='auto')])

fig.update_layout(
    title = 
        {'text':'<b>Топ-5 схожих песен</b><br><span style="color:#666666"></span>'},
    showlegend = False,
    height=650,
    font={
        'family':'Open Sans, light',
        'color':'black',
        'size':14
    },
    plot_bgcolor='rgba(0,0,0,0)',
    xaxis={'categoryorder':'total descending'}
)

fig.show()

Имея вектор каждой песни, давайте посчитаем вектор всего альбома — сложим вектора песен. Затем для такого вектора при помощи модели получим самые близкие по духу и смыслу слова.

def get_word_from_tlist(lst):
    for word in lst:
        word = word[0].split('_')[0]
        print(word, end=' ')

vec_sum = 0
for vec in df.vector:
    vec_sum += vec
sim_word = model.similar_by_vector(vec_sum)
get_word_from_tlist(sim_word)

небо тоска тьма пламень плакать горе печаль сердце солнце мрак

Наверное, это ключевой результат и описание альбома Земфиры всего лишь в 10 словах.

Наконец, построим общую тепловую карту, каждая ячейка которой — результат сравнения косинусной мерой текстов двух треков.

colorscale=[[0.0, "rgba(255,255,255,255)"],
            [0.1, "rgba(229,232,237,255)"],
            [0.2, "rgba(216,222,232,255)"],
            [0.3, "rgba(205,214,228,255)"],
            [0.4, "rgba(182,195,218,255)"],
            [0.5, "rgba(159,178,209,255)"],
            [0.6, "rgba(137,161,200,255)"],
            [0.7, "rgba(107,137,188,255)"],
            [0.8, "rgba(96,129,184,255)"],
            [1.0, "rgba(76,114,176,255)"]]

font_colors = ['black']
x = list(df.name.iloc[::-1])
y = list(df.name)
fig = ff.create_annotated_heatmap(result, x=x, y=y, colorscale=colorscale, font_colors=font_colors)
fig.show()

Результаты анализа и интерпретация данных

Давайте ещё раз посмотрим на всё, что у нас получилось — начнём с облака слов. Нетрудно заметить, что у слов «боль», «невозможно», «сорваться», «растерзаны», «сложно», «терпеть», «любить» размер весьма приличный — всё потому, что такие слова встречаются часто на протяжении всего текста песен:

Одной из самых «разнообразных» песен оказался сингл «крым» — в нём 74% уникальных слов. А в песне «снег идёт» слов совсем мало, поэтому большинство — 82% уникальны. Самой большой песней в альбоме получился трек «таблетки» — суммарно там около 150 слов.

Как было выяснено на прошлом графике, самый «динамичный» трек — «таблетки», целых 37 слов в минуту — практически по слову на каждые две секунды. А самый длинный трек — «абъюз», в нём же и согласно предыдущему графику практически самый низкий процент уникальных слов — 46%.

Топ-5 самых семантически похожих пар текстов:

Ещё мы получили вектор всего альбома и подобрали самые близкие слова. Только посмотрите на них — «тьма», «тоска», «плакать», «горе», «печаль», «сердце» — это же ведь и есть тот перечень слов, который характеризует лирику Земфиры!

небо тоска тьма пламень плакать горе печаль сердце солнце мрак

Финал — тепловая карта. По визуализации заметно, что практически все песни достаточно схожи между собой — косинусная мера у многих пар превышает значение в 0.4.

Выводы

В материале мы провели EDA всего текста нового альбома и при помощи предобученной модели Word2Vec доказали гипотезу — большинство песен «бордерлайна» пронизывают довольно мрачные и тексты. И это нормально, ведь Земфиру мы любим именно за искренность и прямолинейность.

Ранее Ctrl + ↓