Введение в анализ данных

Что такое среднее и как с ним правильно работать.

import numpy as np
import pandas as pd
from random import choices, shuffle
import scipy.stats as sps
import matplotlib.pyplot as plt

import seaborn as sns
sns.set(style='whitegrid', font_scale=1.3, palette='Set2')

1. Кошачий университет

Пусть у нас есть 20 групп котиков-студентов.

Вопрос: сколько котиков в среднем в группе?

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

Проведем эксперимент. Сгенерируем 20 групп, причем в каждой будет от 10 до 30 котиков.

groups_count = 20

groups_size = sps.randint(low=10, high=30).rvs(size=groups_count)
groups_size

array([11, 26, 18, 17, 11, 13, 12, 27, 13, 26, 19, 19, 17, 26, 26, 10, 16,
       29, 17, 19])

Ответ получить очень просто:

groups_size.mean()

18.6

Теперь предположим, что мы хотим провести эксперимент согласно тому, как учат в анализе данных:

• возьмем случайную выборку котиков,
• спросим у них, сколько котиков в их группе,
• усредним полученные ответы.

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

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

students_group = []

for i in range(groups_count):
    students_group += [i] * groups_size[i]

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

random_indexes = np.array(choices(students_group, k=30))
groups_size[random_indexes].mean()

20.7

Хм, среднее получилось больше, чем ответ деканата. Может случайно так получилось? Ведь эксперимент случайный.

Давайте повторим эксперимент 10 раз.

for _ in range(10):
    random_indexes = np.array(choices(students_group, k=30))
    print(np.round(groups_size[random_indexes].mean(), 2))

21.67
21.0
20.67
19.37
21.33
19.3
20.6
22.43
21.43
21.3

Что-то все равно не так...

А что будет если опросить вообще всех котиков? Навернято хоть тут должен получиться тот же ответ.

Попробуем

groups_size[np.array(students_group)].mean()

20.50537634408602

Хм, ответ все равно получился не тот, что сказали нам в деканате... В чем же дело?

Оказывается, дело в том, что выбирая случайного котика мы чаще попадаем на котика из большей группы. В самом деле, если

• $n$ — общее количество котиков,
• $n_j$ — количество котиков в группе $j$,

то вероятность того, что случайно выбранный котик учится в группе $j$ равна $n_j/n$.

Посчитаем теперь средний ответ котика. Пусть $\xi$ — количество студентов в группе у случайно выбранного котика. Тогда $$\mathsf{E} \xi = \frac1n \sum_{котик\ i} \sum_{j} n_j \cdot I\{котик\ i\ из\ группы\ j\} = \frac1n \sum_{j} n_j^2.$$

Проверим ответ:

(groups_size**2).sum() / groups_size.sum()

20.50537634408602

Это среднее можно посчитать другим способом: $$\mathsf{E} \xi = \sum_j n_j \cdot \mathsf{P}(случайный\ котик\ из\ группы\ j).$$

Если рассматривать в качестве данных сами группы, а не котиков, то такой тип усреднения называется взвешенным средним. Общая формула взвешенного среднего чисел $x_1, ..., x_n$ с неотрицательными весами $w_1, ..., w_n$, для которых выполнено $\sum_{i=1}^n w_i=1$, имеет вид $$\sum_{i=1}^n w_i x_i.$$

Выводы:

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

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

---

2. Средняя зарплата

Вопрос: какова средняя зарплата населения?

Важно отличать разные "виды средних":

• обычное среднее арифметическое,
• медиана — значение, слева и справа от которого одинаковое количество элементов,
• мода — самое частое значение.

Разницу между ними наглядно показывает иллюстрация из книги Huff D. How To Lie With Statistics. — New York: W.W. Norton & Company, 1954.

Посмотрим на численном примере. Сгенерируем зарплату 10 000 человек в соответствии с некоторым распределением.

count = 10000
salary = np.abs(sps.t(df=2, scale=100).rvs(size=count))

Посмотрим на описательные статистики. Как мы видим, среднее достаточно сильно отличается от медианы.

pd.Series(salary).describe()

count    10000.000000
mean       143.361444
std        274.727742
min          0.014972
25%         35.979807
50%         81.488075
75%        162.664021
max       8320.602082
dtype: float64

Посмотрим на гистограмму распределения. Однако, в данном случае в выборке есть выбросы — сильно выделяющиеся по сравнению с остальными наблюдения. Видимо, выбросами у нас являются миллиардеры. В анализе данных существуют специальные методы для работы с выбросами, однако не редко их просто выбрасывают из анализа.

Выбросы сильно влияют на вид гистограммы. В данном случае видна широкая часть графика справа, в которой скорее какая-то пустота. Для наглядности стоит рисовать гистограмму в логарифмическом масштабе, на которой явно видно имеющиеся выбросы.

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

plt.subplot(121)
plt.hist(salary, bins=50)
plt.xlabel('Зарплата')
plt.ylabel('Количество человек')
plt.title('Простая гистограмма')

plt.subplot(122)
plt.hist(salary, bins=50, log=True)
plt.xlabel('Зарплата')
plt.ylabel('Количество человек')
plt.title('Гистограмма в логарифмическом масштабе')

plt.show()

Попробуем убрать несколько выбросов и посчитать среднее по оставшимся элементам. Как видим, среднее достаточно сильно уменьшилось, хотя мы выкинули около 10 человек.

salary[salary < 3300].mean()

137.01957408968383

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

3. Когда придет мой автобус? Или каково среднее время ожидания автобуса.

Вы приходите на автобусную остановку. Согласно расписанию автобус ходит каждые 10 минут. Вы засекаете время, и получается, что автобусы обычно приходят через 9-11 минут.

Почему так не везет?

Казалось бы, если автобус ходит каждые 10 минут, то в среднем его нужно ждать 5 минут. Однако, не все так просто. В действительности автобусы приезжают не точно по расписанию, а случайным образом. Оказывается, справедливо следующее утверждение.

Если автобусы приходят с одинаковой интенсивностью в 10 минут и независимо друг от друга, то среднее время ожидания составляет 10 минут.

Это называют парадоксом времени ожидания. Далее мы попробуем экспериментально проверить это утверждение. Теоретическое доказательство этих фактов ожидает вас на 3 курсе.

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

count = 1000000  # количество автобусов
tau = 10  # средний интервал между автобусами

bus_arrival_times = sps.uniform(scale=count*tau).rvs(size=count)
bus_arrival_times = np.sort(bus_arrival_times)

Посмотрим на гистограмму интервалов прибытия автобусов, сравнивая ее с графиком плотность равномерного распределения.

grid = np.linspace(-10, count*tau+10, 1010)

plt.figure(figsize=(8, 5))
plt.hist(bus_arrival_times, bins=50, density=True, label='Данные')
plt.plot(grid, sps.uniform(scale=count*tau).pdf(grid), lw=5, alpha=0.5, 
         label='Плотность $U[0, max]$', color='#FF3300')
plt.xlabel('Время прибытия автобуса')
plt.legend()
plt.show()

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

hours_count = 100000  # количество часов

# количество автобусов за каждый интервал длиной в 1 час
hist = np.histogram(bus_arrival_times, bins=60*np.arange(hours_count))[0]
# количество интервалов, за которые приехало одинаковое количество автобусов
x, y = np.unique(hist, return_counts=True)

Посмотрим на гистограмму распределения количества автобусов за час и сравним его с графиком дискретной плотности пуассоновского распределения с параметром 6. Число 6 есть теоретически ожидаемое количество автобусов в час если их интенсивность движения есть в среднем 1 автобус за 10 минут.

grid = np.arange(25)

plt.figure(figsize=(8, 5))
plt.bar(x, y/y.sum(), label='Данные')
plt.plot(grid, sps.poisson(mu=60/tau).pmf(grid), marker='o', ms=10, lw=5, alpha=0.5, 
         label='Дискр.\nплотность $Pois(60/\\tau)$', color='#FF3300')
plt.xlabel('Количество автобусов в час')
plt.legend()
plt.show()

Теперь рассмотрим средние длины интервалов между прибытиями автобусов. Посчитаем также средний интервал. Как видим, он составляет примерно 10 минут.

intervals = np.diff(bus_arrival_times)
intervals.mean()

9.999999740914484

Посмотрим на гистограмму длин интервалов между прибытиями автобусов и сравним ее с плотностью экспоненциального распределения.

grid = np.linspace(-1, 70, 1000)

plt.figure(figsize=(8, 5))
plt.hist(intervals, bins=50, density=True, label='Данные', range=(0, 70))
plt.plot(grid, sps.expon(scale=tau).pdf(grid), lw=5, alpha=0.5, 
         label='Плотность $Exp(\\tau)$', color='#FF3300')
plt.xlabel('Время между прибытиями автобусов')
plt.legend()
plt.show()

Проведем опрос миллиона пассажиров с целью выяснить, сколько времени они ждали автобус

n_passengers = 1000000  # количество пассажиров

# сгенерируем для каждого пассажира время его прибытия на остановку
passenger_times = sps.uniform(scale=bus_arrival_times.max()).rvs(size=n_passengers)
# найдем время прибытия следующего автобуса поиском по отсортированному массиву
i = np.searchsorted(bus_arrival_times, passenger_times, side='right')
# вычислим интервал ожидания
wait_times = bus_arrival_times[i] - passenger_times

wait_times.mean()

10.000317351938968

Посмотрим на гистограмму времени ожидания прибытия автобуса и сравним ее с плотностью ТОГО ЖЕ САМОГО экспоненциального распределения.

grid = np.linspace(-1, 70, 1000)

plt.figure(figsize=(8, 5))
plt.hist(wait_times, bins=50, density=True, label='Данные', range=(0, 70))
plt.plot(grid, sps.expon(scale=tau).pdf(grid), lw=5, alpha=0.5, 
         label='Плотность $Exp(\\tau)$', color='#FF3300')
plt.xlabel('Время ожидания автобуса')
plt.legend()
plt.show()

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

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

А сколько времени проходит между каждым пятым прибывающим автобусом? Посчитаем эти интервалы времени.

intervals = np.diff(bus_arrival_times[::5])
intervals.mean()

49.99999847766274

Посмотрим на гистограмму длин интервалов между прибытиями каждого 5-го автобуса и сравним ее с плотностью гамма-распределения.

grid = np.linspace(-1, 180, 1000)

plt.figure(figsize=(8, 5))
plt.hist(intervals, bins=50, density=True, label='Данные', range=(0, 180))
plt.plot(grid, sps.gamma(a=5, scale=tau).pdf(grid), lw=5, alpha=0.5, 
         label='Плотность $\\Gamma(5, \\tau)$', color='#FF3300')
plt.xlabel('Время между прибытиями каждого 5-го автобуса')
plt.legend()
plt.show()

Этот пример также показывает, что экпоненциальное распределение является частным случаем гамма-распредедения.