Как найти ковариацию в python

Чем больше я узнаю людей, тем больше мне нравится моя собака.

 —Марк Твен

В предыдущих сериях постов для начинающих из ремикса книги Генри Гарнера «Clojure для исследования данных» (Clojure for Data Science) на языке Python мы рассмотрели методы описания выборок с точки зрения сводных статистик и методов статистического вывода из них параметров популяции. Такой анализ сообщает нам нечто о популяции в целом и о выборке в частности, но он не позволяет нам делать очень точные утверждения об их отдельных элементах. Это связано с тем, что в результате сведения данных всего к двум статистикам — среднему значению и стандартному отклонению — теряется огромный объем информации.

Нам часто требуется пойти дальше и установить связь между двумя или несколькими переменными либо предсказать одну переменную при наличии другой. И это подводит нас к теме данной серии из 5 постов — исследованию корреляции и регрессии. Корреляция имеет дело с силой и направленностью связи между двумя или более переменными. Регрессия определяет природу этой связи и позволяет делать предсказания на ее основе.

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

О данных

В этой серии постов используются данные, любезно предоставленные компанией Guardian News and Media Ltd., о спортсменах, принимавших участие в Олимпийских Играх 2012 г. в Лондоне. Эти данные изначально были взяты из блога газеты Гардиан.

Обследование данных

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

Файл all-london-2012-athletes.tsv достаточно небольшой. Мы можем обследовать данные при помощи pandas, как мы делали в первой серии постов «Python, исследование данных и выборы», воспользовавшись функцией read_csv:

def load_data():
    return pd.read_csv('data/ch03/all-london-2012-athletes-ru.tsv', 't')
                                            
def ex_3_1():
    '''Загрузка данных об участниках 
       олимпийских игр в Лондоне 2012 г.'''
    return load_data()

Если выполнить этот пример в консоли интерпретатора Python либо в блокноте Jupyter, то вы должны увидеть следующий ниже результат:

Столбцы данных (нам повезло, что они ясно озаглавлены) содержат следующую информацию:

• ФИО атлета

• страна, за которую он выступает

• возраст, лет

• рост, см.

• вес, кг.

• пол «М» или «Ж»

• дата рождения в виде строки

• место рождения в виде строки (со страной)

• число выигранных золотых медалей

• число выигранных серебряных медалей

• число выигранных бронзовых медалей

• всего выигранных золотых, серебряных и бронзовых медалей

• вид спорта, в котором он соревновался

• состязание в виде списка, разделенного запятыми

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

Визуализация данных

В первую очередь мы рассмотрим разброс роста спортсменов на Олимпийских играх 2012 г. в Лондоне. Изобразим эти значения роста в виде гистограммы, чтобы увидеть характер распределения данных, не забыв сначала отфильтровать пропущенные значения:

def ex_3_2():
    '''Визуализация разброса значений 
       роста спортсменов на гистограмме'''
    df = load_data()
    df['Рост, см'].hist(bins=20)
    plt.xlabel('Рост, см.')
    plt.ylabel('Частота')
    plt.show()

Этот пример сгенерирует следующую ниже гистограмму:

Как мы и ожидали, данные приближенно нормально распределены. Средний рост спортсменов составляет примерно 177 см. Теперь посмотрим на распределение веса олимпийских спортсменов:

def ex_3_3():
    '''Визуализация разброса значений веса спортсменов'''
    df = load_data()
    df['Вес'].hist(bins=20)
    plt.xlabel('Вес')
    plt.ylabel('Частота')
    plt.show()

Приведенный выше пример сгенерирует следующую ниже гистограмму:

Данные показывают четко выраженную асимметрию. Хвост с правой стороны намного длиннее, чем с левой, и поэтому мы говорим, что асимметрия — положительная. Мы можем оценить асимметрию данных количественно при помощи функции библиотеки pandas skew:

def ex_3_4():
    '''Вычисление асимметрии веса спортсменов'''
    df = load_data()
    swimmers = df[ df['Вид спорта'] == 'Swimming']
    return swimmers['Вес'].skew()

0.23441459903001483

К счастью, эта асимметрия может быть эффективным образом смягчена путем взятия логарифма веса при помощи функции библиотеки numpy np.log:

def ex_3_5():
    '''Визуализация разброса значений веса спортсменов на
       полулогарифмической гистограмме с целью удаления 
       асимметрии'''
    df = load_data()
    df['Вес'].apply(np.log).hist(bins=20)
    plt.xlabel('Логарифмический вес')
    plt.ylabel('Частота')
    plt.show()

Этот пример сгенерирует следующую ниже гистограмму:

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

Логнормальное распределение

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

Логарифм показывает степень, в которую должно быть возведено фиксированное число (основание) для получения данного числа. Изобразив логарифмы на графике в виде гистограммы, мы показали, что эти степени приближенно нормально распределены. Логарифмы обычно берутся по основанию 10 или основанию e, трансцендентному числу, приближенно равному 2.718. В функции библиотеки numpy np.log и ее инверсии np.exp используется основание e. Выражение log_e также называется натуральным логарифмом, или ln, из-за свойств, делающих его особенно удобным в исчислении.

Логнормальное распределение обычно имеет место в процессах роста, где темп роста не зависит от размера. Этот феномен известен как закон Джибрэта, который был cформулирован в 1931 г. Робертом Джибрэтом, заметившим, что он применим к росту фирм. Поскольку темп роста пропорционален размеру, более крупные фирмы демонстрируют тенденцию расти быстрее, чем фирмы меньшего размера.

^{Нормальное распределение случается в ситуациях, где много мелких колебаний, или вариаций, носит суммирующий эффект, тогда как логнормальное распределение происходит там, где много мелких вариаций имеет мультипликативный эффект.}

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

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

Визуализация корреляции

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

def swimmer_data():
    '''Загрузка данных роста и веса только олимпийских пловцов'''
    df = load_data()
    return df[df['Вид спорта'] == 'Swimming'].dropna()

def ex_3_6():
    '''Визуализация корреляции между ростом и весом'''
    df = swimmer_data()
    xs = df['Рост, см']
    ys = df['Вес'].apply( np.log )
    pd.DataFrame(np.array([xs,ys]).T).plot.scatter(0, 1, s=12, grid=True)
    plt.xlabel('Рост, см.')
    plt.ylabel('Логарифмический вес')
    plt.show()

Этот пример сгенерирует следующий ниже график:

Результат ясно показывает, что между этими двумя переменными имеется связь. График имеет характерно смещенную эллиптическую форму двух коррелируемых, нормально распределенных переменных с центром вокруг среднего значения. Следующая ниже диаграмма сравнивает график рассеяния с распределениями вероятностей роста и логарифма веса:

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

Генерирование джиттера

Поскольку каждое значение округлено до ближайшего сантиметра или килограмма, то значение, записанное как 180 см, на самом деле может быть каким угодно между 179.5 и 180.5 см, тогда как значение 80 кг на самом деле может быть каким угодно между 79.5 и 80.5 кг. Для создания случайных искажений, мы можем добавить случайные помехи в каждую точку данных роста в диапазоне между -0.5 и 0.5 и в том же самом диапазоне проделать с точками данных веса (разумеется, это нужно cделать до того, как мы возьмем логарифм значений веса):

def jitter(limit):
    '''Генератор джиттера (произвольного сдвига точек данных)'''
    return lambda x: random.uniform(-limit, limit) + x

def ex_3_7():
    '''Визуализация корреляции между ростом и весом с джиттером'''
    df = swimmer_data()
    xs = df['Рост, см'].apply(jitter(0.5))
    ys = df['Вес'].apply(jitter(0.5)).apply(np.log)
    pd.DataFrame(np.array([xs,ys]).T).plot.scatter(0, 1, s=12, grid=True)
    plt.xlabel('Рост, см.')
    plt.ylabel('Логарифмический вес')
    plt.show()

График с джиттером выглядит следующим образом:

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

Ковариация

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

Если у нас имеется два ряда чисел, X и Y, то их отклонения от среднего значения составляют:

Здесь x_i — это значение X с индексом i, y_i — значение Y с индексом i, x̅ — среднее значение X, и y̅ — среднее значение Y. Если X и Y проявляют тенденцию изменяться вместе, то их отклонения от среднего будет иметь одинаковый знак: отрицательный, если они — меньше среднего, положительный, если они больше среднего. Если мы их перемножим, то произведение будет положительным, когда у них одинаковый знак, и отрицательным, когда у них разные знаки. Сложение произведений дает меру тенденции этих двух переменных отклоняться от среднего значения в одинаковом направлении для каждой заданной выборки.

Ковариация определяется как среднее этих произведений:

На чистом Python ковариация вычисляется следующим образом:

def covariance(xs, ys):
    '''Вычисление ковариации (несмещенная, т.е. n-1)'''
    dx = xs - xs.mean() 
    dy = ys - ys.mean()
    return (dx * dy).sum() / (dx.count() - 1)

В качестве альтернативы, мы можем воспользоваться функцией pandas cov:

df['Рост, см'].cov(df['Вес'])

1.3559273321696459

Ковариация роста и логарифма веса для наших олимпийских пловцов равна 1.356, однако это число сложно интерпретировать. Единицы измерения здесь представлены произведением единиц на входе.

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

^{Стандартная оценка, англ. standard score, также z-оценка — это относительное число стандартных отклонений, на которые значение переменной отстоит от среднего значения. Положительная оценка показывает, что переменная находится выше среднего, отрицательная — ниже среднего. Это безразмерная величина, получаемая при вычитании популяционного среднего из индивидуальных значений и деления разности на популяционное стандартное отклонение.}

Корреляция Пирсона

Корреляция Пирсона часто обозначается переменной r и вычисляется следующим образом, где отклонения от среднего dx_i и dy_i вычисляются как и прежде:

Поскольку для переменных X и Y стандартные отклонения являются константными, уравнение может быть упрощено до следующего, где σ_x и σ_y — это стандартные отклонения соответственно X и Y:

В таком виде формула иногда упоминается как коэффициент корреляции смешанных моментов Пирсона или попросту коэффициент корреляции и, как правило, обозначается буквой r.

Ранее мы уже написали функции для вычисления стандартного отклонения. В сочетании с нашей функцией с вычислением ковариации получится следующая ниже имплементация корреляции Пирсона:

def variance(xs):
    '''Вычисление корреляции,
       несмещенная дисперсия при n <= 30'''
    x_hat = xs.mean()
    n = xs.count()
    n = n - 1 if n in range( 1, 30 ) else n  
    return sum((xs - x_hat) ** 2) / n

def standard_deviation(xs):
    '''Вычисление стандартного отклонения'''
    return np.sqrt(variance(xs))

def correlation(xs, ys): 
    '''Вычисление корреляции'''
    return covariance(xs, ys) / (standard_deviation(xs) * 
                                 standard_deviation(ys))

В качестве альтернативы мы можем воспользоваться функцией pandas corr:

df['Рост, см'].corr(df['Вес'])

Поскольку стандартные оценки безразмерны, то и коэффициент корреляции r тоже безразмерен. Если r равен -1.0 либо 1.0, то переменные идеально антикоррелируют либо идеально коррелируют.

Правда, если r = 0, то с необходимостью вовсе не следует, что переменные не коррелируют. Корреляция Пирсона измеряет лишь линейные связи. Как продемонстрировано на следующих графиках, между переменными может существовать еще некая нелинейная связь, которую r не объясняет:

Отметим, что корреляция центрального примера не определена, потому что стандартное отклонение y = 0. Поскольку наше уравнение для r содержало бы деление ковариации на 0, то результат получается бессмысленным. В этом случае между переменными не может быть никакой корреляции; y всегда будет иметь среднее значение. Простое обследование стандартных отклонений это подтвердит.

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

def ex_3_8():
    '''Вычисление корреляции средствами pandas
       на примере данных роста и веса'''
    df = swimmer_data()
    return df['Рост, см'].corr( df['Вес'].apply(np.log))

0.86748249283924894

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

Выборочный r и популяционный ρ

Аналогично среднему значению и стандартному отклонению, коэффициент корреляции является сводной статистикой. Он описывает выборку; в данном случае, выборку спаренных значений: роста и веса. Коэффициент корреляции известной выборки обозначается буквой r, тогда как коэффициент корреляции неизвестной популяции обозначается греческой буквой ρ (рхо).

Как мы убедились в предыдущей серии постов о тестировании гипотез, мы не должны исходить из того, что результаты, полученные в ходе измерения нашей выборки, применимы к популяции в целом. К примеру, наша популяция может состоять из всех пловцов всех недавних Олимпийских игр. И будет совершенно недопустимо обобщать, например, на другие олимпийские виды спорта, такие как тяжелая атлетика или фитнес-плавание.

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

• Размера выборки

• Величины r

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

Проверка статистических гипотез

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

В первую очередь, мы должны сформулировать две гипотезы, нулевую гипотезу и альтернативную:

H₀ — это гипотеза, что корреляция в популяции нулевая. Другими словами, наше консервативное представление состоит в том, что измеренная корреляция целиком вызвана случайной ошибкой при отборе.

H₁ — это альтернативная возможность, что корреляция в популяции не нулевая. Отметим, что мы не определяем направление корреляции, а только что она существует. Это означает, что мы выполняем двустороннюю проверку.

Стандартная ошибка коэффициента корреляции r по выборке задается следующей формулой:

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

Мы можем снова воспользоваться t-распределением и вычислить t-статистику:

В приведенной формуле df — это степень свободы наших данных. Для проверки корреляции степень свободы равна n — 2, где n — это размер выборки. Подставив это значение в формулу, получим:

В итоге получим t-значение 102.21. В целях его преобразования в p-значение мы должны обратиться к t-распределению. Библиотека scipy предоставляет интегральную функцию распределения (ИФР) для t-распределения в виде функции stats.t.cdf, и комплементарной ей (1-cdf) функции выживания stats.t.sf. Значение функции выживания соответствует p-значению для односторонней проверки. Мы умножаем его на 2, потому что выполняем двустороннюю проверку:

def t_statistic(xs, ys):
    '''Вычисление t-статистики'''
    r = xs.corr(ys)  # как вариант, correlation(xs, ys)
    df = xs.count() - 2
    return r * np.sqrt(df / 1 - r ** 2)

def ex_3_9():
    '''Выполнение двухстороннего t-теста'''
    df = swimmer_data()
    xs = df['Рост, см']
    ys = df['Вес'].apply(np.log)
    t_value = t_statistic(xs, ys)
    df = xs.count() - 2 
    p = 2 * stats.t.sf(t_value, df)  # функция выживания 
    return {'t-значение':t_value, 'p-значение':p}

{'p-значение': 1.8980236317815443e-106, 't-значение': 25.384018200627057}

P-значение настолько мало, что в сущности равно 0, означая, что шанс, что нулевая гипотеза является истинной, фактически не существует. Мы вынуждены принять альтернативную гипотезу о существовании корреляции.

Интервалы уверенности

Установив, что в более широкой популяции, безусловно, существует корреляция, мы, возможно, захотим количественно выразить диапазон значений, внутри которого, как мы ожидаем, будет лежать параметр ρ, вычислив для этого интервал уверенности. Как и в случае со средним значением в предыдущей серии постов, интервал уверенности для r выражает вероятность (выраженную в %), что параметр ρ популяции находится между двумя конкретными значениями.

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

Приведенный выше график показывает отрицательно скошенное распределение r-выборок для параметра ρ, равного 0.6.

К счастью, трансформация под названием z-преобразование Фишера стабилизирует дисперсию r  по своему диапазону. Она аналогична тому, как наши данные о весе спортсменов стали нормально распределенными, когда мы взяли их логарифм.

Уравнение для z-преобразования следующее:

Стандартная ошибка z равна:

Таким образом, процедура вычисления интервалов уверенности состоит в преобразовании r в z с использованием z-преобразования, вычислении интервала уверенности в терминах стандартной ошибки SE_z и затем преобразовании интервала уверенности в r.

В целях вычисления интервала уверенности в терминах SE_z, мы можем взять число стандартных отклонений от среднего, которое дает нам требуемый уровень доверия. Обычно используют число 1.96, так как оно является числом стандартных отклонений от среднего, которое содержит 95% площади под кривой. Другими словами, 1.96 стандартных ошибок от среднего значения выборочного r содержит истинную популяционную корреляцию ρ с 95%-ой определенностью.

Мы можем убедиться в этом, воспользовавшись функцией scipy stats.norm.ppf. Она вернет стандартную оценку, связанную с заданной интегральной вероятностью в условиях односторонней проверки.

Однако, как показано на приведенном выше графике, мы хотели бы вычесть ту же самую величину, т.е. 2.5%, из каждого хвоста с тем, чтобы 95%-й интервал уверенности был центрирован на нуле. Для этого при выполнении двусторонней проверки нужно просто уменьшить разность наполовину и вычесть результат из 100%. Так что, требуемый уровень доверия в 95% означает, что мы обращаемся к критическому значению 97.5%:

def critical_value(confidence, ntails): # ДИ и число хвостов
    '''Расчет критического значения путем
       вычисления квантиля и получения 
       для него нормального значения'''
    lookup = 1 - ((1 - confidence) / ntails) 
    return stats.norm.ppf(lookup, 0, 1)  # mu=0, sigma=1

critical_value(0.95, 2)

1.959963984540054

Поэтому наш 95%-й интервал уверенности в z-пространстве для ρ задается следующей формулой:

Подставив в нашу формулу z_r и SE_z, получим:

Для r=0.867 и n=859 она даст нижнюю и верхнюю границу соответственно 1.137 и 1.722. В целях их преобразования из z-оценок в r-значения, мы используем следующее обратное уравнение z-преобразования:

Преобразования и интервал уверенности можно вычислить при помощи следующего исходного кода:

def z_to_r(z):
    '''Преобразование z-оценки обратно в r-значение'''
    return (np.exp(z*2) - 1) / (np.exp(z*2) + 1)

def r_confidence_interval(crit, xs, ys): 
    '''Расчет интервала уверенности
       для критического значения и данных'''
    r   = xs.corr(ys)
    n   = xs.count()
    zr  = 0.5 * np.log((1 + r) / (1 - r)) 
    sez = 1 / np.sqrt(n - 3)
    return (z_to_r(zr - (crit * sez))), (z_to_r(zr + (crit * sez)))

def ex_3_10():
    '''Расчет интервала уверенности
       на примере данных роста и веса'''
    df = swimmer_data()
    X = df['Рост, см']
    y = df['Вес'].apply(np.log)
    interval = r_confidence_interval(1.96, X, y) 
    print('Интервал уверенности (95%):', interval)

Интервал уверенности (95%): (0.8499088588880347, 0.8831284878884087)

В результате получаем 95%-й интервал уверенности для ρ, расположенный между 0.850 и 0.883. Мы можем быть абсолютно уверены в том, что в более широкой популяции олимпийских пловцов существует сильная положительная корреляция между ростом и весом.

Примеры исходного кода для этого поста находятся в моем репо на Github. Все исходные данные взяты в репозитории автора книги.

В следующем посте, посте №2, будет рассмотрена сама тема серии — регрессия и приемы оценивания ее качества.

В предыдущей части мы уже начали разбирать описательную статистику, а именно: центральные метрики и метрики оценки вариативности. В этой же части мы будем говорить той же описательной статистики, но уже разберем корреляции между парами данных. Корреляция – это статистическая взаимосвязь между двумя или более случайными величинами. И есть две статистические меры, которые характеризуют корреляцию между наборами данных — ковариация и коэффициент корреляции. Но прежде, чем нам начать с ними работать, необходимо установить и импортировать нужные библиотеки, а именно: math, statistics, numpy, scipy.stats, pandas и matplotlib.pyplot. После этого, мы формируем два списка Python, которые будем использовать для получения соответствующих массивов NumPy и серии Pandas:

x = list(range(-10, 11))
y = [0, 2, 2, 2, 2, 3, 3, 6, 7, 4, 7, 6, 6, 9, 4, 5, 5, 10, 11, 12, 14]
x_, y_ = np.array(x), np.array(y)
x__, y__ = pd.Series(x_), pd.Series(y_)

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

Ковариации

Выборочная ковариация — это мера, которая количественно определяет силу и направление взаимосвязи между парой переменных:
• Если корреляция положительная, то и ковариация тоже положительная. Более сильное отношение соответствует более высокой ценности ковариации.
• Если корреляция отрицательна, то и ковариация также отрицательна. Более сильное отношение соответствует более низкому (или более высокому абсолютному) значению ковариации.
• Если корреляция слабая, то ковариация близка к нулю.
Процесс расчета ковариации представлен ниже:

#Ковариация
#Расчет ковариации в чистом Python
mean_x, mean_y = sum(x) / n, sum(y) / n
cov_xy = (sum((x[k] - mean_x) * (y[k] - mean_y) for k in range(n))/ (n - 1))
print(f'Расчет ковариации в чистом Python: {cov_xy}')
#Расчет ковариации с помощью NumPy
cov_matrix = np.cov(x_, y_)
print(f'Расчет ковариационной матрицы с помощью NumPy функцией cov(): {cov_matrix}')
print(f'Проверка, что левый элемент ковариационной матрицы — это ковариация x и x или дисперсия x, '
      f'а правый элемент — ковариация y и y или дисперсия y: {x_.var(ddof=1)} и {y_.var(ddof=1)}')
cov_xy = cov_matrix[0, 1]
cov_xy2 = cov_matrix[1, 0]
print(f'Проверка, что два других элемента ковариационной матрицы равны '
      f'и представляют фактическую ковариацию между x и y: {cov_xy} и {cov_xy2}')
#Расчет ковариации с помощью Pandas
cov_xy = x__.cov(y__)
cov_xy3 = y__.cov(x__)
print(f'Расчет ковариации с помощью Pandas методом .cov(): {cov_xy} и {cov_xy3}')

Первым делом была найдена ковариация в чистом Python, где сначала необходимо найти среднее значение х и у, а затем применить математическую формулу. Но можно применить функцию cov() библиотеки NumPy, которая возвращает ковариационную матрицу, что и было сделано во втором примере. Обратите внимание, cov() имеет необязательные параметры bias (по умолчанию False) и ddof (по умолчанию None). Их значения по умолчанию подходят для получения образца ковариационной матрицы. Верхний левый элемент ковариационной матрицы — это ковариация x и x или дисперсия x. Точно так же нижний правый элемент — y и y или дисперсия y. А два других элемента ковариационной матрицы равны и представляют фактическую ковариацию между x и y. Как проверить и убедиться, что это правда, разобрано в примере. Использовав np.cov() мы получили то же значение ковариации, что и с чистым Python. Помимо этого, можно было использовать метод .cov() библиотеки Pandas, что и было сделано в третьем примере, где для одного объекта Series вызывается .cov() и передает другой объект в качестве первого аргумента.

Коэффициент корреляции

Коэффициент корреляции или коэффициент корреляции Пирсона — произведение, обозначается символом ??. Коэффициент является еще одним показателем корреляции между данными. К нему надо относиться как к стандартизированной ковариации. Вот несколько важных замечаний:
• ? > 0 указывает на положительную корреляцию.
• ? < 0 указывает на отрицательную корреляцию.
• r = 1 является максимально возможным значением ?. Это свидетельство полной линейной зависимости между переменными.
• r = −1 является минимально возможным значением ?. Это свидетельство полного отсутствия линейной зависимости между переменными.
• r ≈ 0 или когда around около нуля, означает, что корреляция между переменными отсутствует.
Процесс расчета коэффициента корреляции представлен ниже:

#Коэффициент корреляции
#Расчет коэффициента корреляции в чистом Python
var_x = sum((item - mean_x)**2 for item in x) / (n - 1)
var_y = sum((item - mean_y)**2 for item in y) / (n - 1)
std_x, std_y = var_x ** 0.5, var_y ** 0.5
r = cov_xy / (std_x * std_y)
print(f'Расчет коэффициента корреляции в чистом Python: {r}')
#Расчет коэффициента корреляции с помощью scipy.stats
r, p = scipy.stats.pearsonr(x_, y_)
print(f'Расчет коэффициента корреляции и p-value, используя функцию pearsonr() в scipy.stats: {r} и {p} ')
scipy.stats.linregress(x_, y_)
print(f'Расчет коэффициента корреляции с помощью scipy.stats.linregress(): {scipy.stats.linregress(x_, y_)}')
result = scipy.stats.linregress(x_, y_)
r = result.rvalue
print(f'Получение доступа к определенным значениям из результата linregress(), включая коэффициент корреляции, используя точечную запись: {r}')
#Расчет коэффициента корреляции с помощью Pandas
r = x__.corr(y__)
r1 = y__.corr(x__)
print(f'Расчет коэффициента корреляции методом .corr() библиотеки Pandas: {r} и {r1}')

В первом примере показано, как рассчитать это коэффициент в чистом Python. Для этого нам понадобятся средние значения (mean_x и mean_y) и стандартные отклонения (std_x, std_y) для наборов данных x и y, а также их ковариация cov_xy. Далее, этот коэффициент и значение p-value (значение, показывающее принимается или отклоняется гипотеза) были рассчитан с помощью функцию pearsonr() в scipy.stats. Первое значение — это коэффициент корреляции между x_ и x_. Второй элемент — это коэффициент корреляции между y_ и y_. Их значения равны 1,0. В третьем примере коэффициент корреляции был найдем с помощью scipy.stats.linregress(). linregress() принимает x_ и y_, вычисляет линейную регрессию и возвращает результаты — наклон и точка пересечения определяют уравнение прямой регрессии, а rvalue — коэффициент корреляции. Чтобы получить доступ к определенным значениям из результата linregress(), включая коэффициент корреляции, необходимо использовать точечную запись, что показано в примере. И в четвертом примере был использован метод .corr() библиотеки Pandas, который вызывается для одного объекта Series и передает другой объект в качестве первого аргумента.

Работа с 2D данными

Axis
В статистике очень часто работают с 2D данными. NumPy и SciPy предоставляют комплексные средства для работы с ними, а Pandas имеет специальный класс DataFrame для обработки 2D данных. И прежде, чем нам работать с ними, их необходимо создать:

#Работа с данными 2D (таблицы)
#Axis
#Создание 2D массива с помощью Numpy
a = np.array(2, 3, 1],[4, 9, 2], [8, 27, 4], [16, 1, 1], [2, 3, 1એ)
print(f'Вывод, созданного с помощью Numpy, 2D массива: {a}')

Теперь у нас есть набор 2D данных, который мы будем использовать в этом разделе. Можно применять к нему статистические функции и методы Python так же, как к данным 1D, используя или не используя при этом необязательный параметр axis, то есть ось. Ось может принимать любое из следующих значений:
• axis = None — расчет статистики по всем данным в массиве, как в приведенном выше примере. Такое поведение часто используется по умолчанию в NumPy.
• axis = 0 — расчет статистики для каждого столбца массива. Такое поведение часто используется по умолчанию для статистических функций SciPy.
• axis = 1 — расчет статистики для каждой строки массива.
Внизу представлен код, в котором я попробовал использовать к нашему 2D данным статистические функции и методы Python, как с параметром axis, так и без:

#Использование статистических функций и методов Python к 2d массиву с необязательным параметров axis
np.mean(a, axis=0)
a.mean(axis=1)
print(f'Вывод среднего значения 2D массива методом NumPy на оси = 0 и на оси = 1 соответствнно: {np.mean(a, axis=0)} и {a.mean(axis=1)}')
np.median(a, axis=0)
np.median(a, axis=1)
print(f'Вывод медианы 2D массива методом NumPy на оси = 0 и на оси = 1 соответствнно: {np.median(a, axis=0)} и {np.median(a, axis=1)}')
scipy.stats.gmean(a)  # Default: axis=0
print(f'Вывод среднего геометрического значения 2D массива функцией SciPy на оси = 0 : {scipy.stats.gmean(a)}')
scipy.stats.gmean(a, axis=None)
print(f'Вывод среднего геометрического значения для всего 2D массива функцией SciPy: {scipy.stats.gmean(a, axis=None)}')

DataFrames
Класс DataFrame является одним из основных типов данных Pandas. С ним очень удобно работать, потому что в нем есть метки для строк и столбцов. Сам код представлен внизу:

#DataFrames
row_names = ['first', 'second', 'third', 'fourth', 'fifth']
col_names = ['A', 'B', 'C']
df = pd.DataFrame(a, index=row_names, columns=col_names)
print(f'Вывод класса DataFrame с ранее созданным 2d массивом: {df}')
df.mean()
df.var()
print(f'Вывод среднего значения и несмещенной дисперсии для всего класса DataFrame с ранее созданным 2d массивом: {df.mean()} и {df.var()}')
df.mean(axis=1)
print(f'Вывод среднего значения класса DataFrame с ранее созданным 2d массивом по оси = 1: {df.mean(axis=1)}')
df['A']
print(f'Пример изоляции класса DataFrame с ранее созданным 2d массивом по столбцу A: {df["A"]}')
df['A'].mean()
print(f'Пример расчета среднего значения класса DataFrame с ранее созданным 2d массивом по столбцу A: {df["A"].mean()}')

Первым делом был создан DataFrame использовав при этом массив a. На практике имена столбцов имеют значение и должны быть описательными. Имена строк иногда указываются автоматически как 0, 1 и т. Д. Вы можете указать их явно с помощью параметра index, хотя вы можете свободно опускать index, если хотите. Методы DataFrame очень похожи на методы Series, хотя их поведение отличается. Если вы вызываете методы статистики Python без аргументов, то DataFrame будет возвращать результаты для каждого столбца. Если же вы хотите получить результаты для каждой строки, просто укажите параметр axis = 1. Также можно изолировать каждый столбец DataFrame и применять к нему соответствующие методы. Процесс применения статистических методов и функций Python к DataFrame представлен в примере.

Визуализация данных

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

Box Plots

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

#Визуализация данных
#Box Plots
np.random.seed(seed=0)
x = np.random.randn(1000)
y = np.random.randn(100)
z = np.random.randn(10)
fig, ax = plt.subplots()
ax.boxplot((x, y, z), vert=False, showmeans=True, meanline=True,
           labels=('x', 'y', 'z'), patch_artist=True,
           medianprops={'linewidth': 2, 'color': 'purple'},
           meanprops={'linewidth': 2, 'color': 'red'})
plt.show()

Первый оператор инициализирует генератора случайных чисел NumPy с помощью seed(). Поэтому при каждом запуске скрипта могут получаться одинаковые результаты. Не нужно устанавливать начальное значение и, если его не указывать, то каждый раз будут получаться разные результаты. Другие операторы создают три массива NumPy с нормально распределенными псевдослучайными числами. x относится к массиву с 1000 элементами, y имеет 100, а z содержит 10 элементов. Параметры .boxplot():
• х ваши данные.
• vert устанавливает горизонтальную ориентацию графика, когда False. Ориентация по умолчанию — вертикальная.
• showmeans показывает среднее значение ваших данных, когда True.
• meanline представляет среднее в виде линии, когда истина. Представлением по умолчанию является точка.
• labels: метки ваших данных.
• patch_artist определяет, как рисовать график.
• medianprops обозначает свойства линии, представляющей медиану.
• meanprops указывает свойства линии или точки, представляющей среднее значение.
Есть и другие параметры, но их анализ выходит за рамки данного руководства. Вы можете увидеть три сюжета. Каждый из них соответствует одному набору данных (x, y или z) и показывает следующее:
• Среднее значение — это красная пунктирная линия.
• Медиана — это фиолетовая линия.
• Первый квартиль — левый край синего прямоугольника.
• Третий квартиль — это правый край синего прямоугольника.
• Межквартильный диапазон — это длина синего прямоугольника.
• Диапазон — всё слева направо.
• Выбросы — точки слева и справа.
Сюжетная диаграмма может показывать столько информации на одном рисунке!
Результат:

Гистограммы

Гистограмма особенно полезна, когда в наборе данных содержится большое количество уникальных значений. Гистограмма делит значения из отсортированного набора данных на интервалы, также называемые ячейками. Часто все лотки имеют одинаковую ширину, но это не обязательно так. Значения нижней и верхней границ ячейки называются ребрами ячейки. Частота представляет собой одно значение, которое соответствует каждому бину. Это количество элементов набора данных со значениями между краями корзины. По договоренности, все корзины, кроме самой правой, наполовину открыты. Они включают значения, равные нижним границам, но исключают значения, равные верхним границам. Крайняя правая корзина закрыта, так как включает обе границы. Если вы разделите набор данных с ребрами 0, 5, 10 и 15, то есть три элемента:
• Первый и самый левый столбец содержит значения, большие или равные 0 и меньшие 5.
• Второй контейнер содержит значения, большие или равные 5 и меньшие 10.
• Третий и самый правый контейнер содержит значения, большие или равные 10 и меньшие или равные 15.
Сам код:

#Гистограммы
hist, bin_edges = np.histogram(x, bins=10)
fig, ax = plt.subplots()
ax.hist(x, bin_edges, cumulative=False)
ax.set_xlabel('x')
ax.set_ylabel('Frequency')
plt.show()

Функция np.histogram() — это удобный способ получить данные для гистограмм. Он берет массив с вашими данными и количеством (или ребрами) бинов и возвращает два массива NumPy:
• hist содержит частоту или количество элементов, соответствующих каждому бину.
• bin_edges содержит ребра или границы корзины.
Что рассчитывает histogram(), граф .hist() может показать графически. Первый аргумент .hist() — это последовательность с вашими данными. Второй аргумент определяет края бинов. Третий отключает возможность создания гистограммы с накопленными значениями. Результат представлен внизу:

Вы можете видеть края корзины на горизонтальной оси и частоты на вертикальной оси.
При значении аргумента cumulative = True в .hist() можно получить гистограмму с совокупным количеством элементов:

Частота первого и самого левого лотка — это количество элементов в этом лотке. Частота второго бина — это сумма количества элементов в первом и втором бинах. Другие контейнеры отрисовываются аналогично. Наконец, частота последнего и самого правого бина — это общее количество элементов в наборе данных (в данном случае 1000).

Pie Charts круговые диаграммы

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

#Pie Charts круговые диаграммы
x, y, z = 128, 256, 1024
fig, ax = plt.subplots()
ax.pie((x, y, z), labels=('x', 'y', 'z'), autopct='%1.1f%%')
plt.show()

Первым делом мы определяем данные, связанные с тремя метками. Далее, Первый аргумент .pie() — данные, а второй — последовательность соответствующих меток. autopct определяет формат относительных частот, показанных на рисунке. Результат:

Круговая диаграмма показывает x как наименьшую часть круга, y как следующую наибольшую, а затем z как самую большую часть. Проценты обозначают относительный размер каждого значения по сравнению с их суммой.

Bar Charts

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

#Bar Charts
x = np.arange(21)
y = np.random.randint(21, size=21)
err = np.random.randn(21)
fig, ax = plt.subplots()
ax.bar(x, y, yerr=err)
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.show()

В данном коде мы использовали np.arange(), чтобы получить x или массив последовательных целых чисел от 0 до 20. Мы использовали это для представления меток. y — это массив равномерно распределенных случайных целых чисел, также от 0 до 20. Этот массив будет представлять частоты. err содержит нормально распределенные числа с плавающей точкой, которые являются ошибками. Эти значения не являются обязательными. Далее создаем гистограмму с помощью .bar(), если вам нужны вертикальные столбцы или .barh(), если вам нужны горизонтальные столбцы. Результат:

X-Y участки

Диаграмма x-y или диаграмма рассеяния представляет пары данных из двух наборов данных. Горизонтальная ось x показывает значения из набора x, а вертикальная ось y показывает соответствующие значения из набора y. При желании вы можете включить линию регрессии и коэффициент корреляции. Код построения представлен внизу:

#X-Y участки
x = np.arange(21)
y = 5 + 2 * x + 2 * np.random.randn(21)
slope, intercept, r, *__ = scipy.stats.linregress(x, y)
line = f'Regression line: y={intercept:.2f}+{slope:.2f}x, r={r:.2f}'
fig, ax = plt.subplots()
ax.plot(x, y, linewidth=0, marker='s', label='Data points')
ax.plot(x, intercept + slope * x, label=line)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.legend(facecolor='white')
plt.show()

Сначала мы генерируем два набора данных и выполняем линейную регрессию с помощью scipy.stats.linregress().Набор данных x снова является массивом с целыми числами от 0 до 20. y вычисляется как линейная функция x, искаженная некоторым случайным шумом. linregress возвращает несколько значений. Вам понадобится наклон и точка пересечения линии регрессии, а также коэффициент корреляции r. Затем применяем .plot(), чтобы получить график x-y. Результат:

Схемы Зоны активности

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

#Схемы Зоны активности
matrix = np.cov(x, y).round(decimals=2)
fig, ax = plt.subplots()
ax.imshow(matrix)
ax.grid(False)
ax.xaxis.set(ticks=(0, 1), ticklabels=('x', 'y'))
ax.yaxis.set(ticks=(0, 1), ticklabels=('x', 'y'))
ax.set_ylim(1.5, -0.5)
for i in range(2):
    for j in range(2):
        ax.text(j, i, matrix[i, j], ha='center', va='center', color='w')
plt.show()

Здесь тепловая карта содержит метки «x» и «y», а также числа из ковариационной матрицы. Результат:

Желтое поле представляет самый большой элемент из матрицы 130.34, а фиолетовое соответствует наименьшему элементу 38.5. Синие квадраты между ними связаны со значением 77,8.

Вывод

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