Описательная статистика на Python: числовые, нечисловые и смешанные данные

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

Что такое описательная статистика?

Итак, данные собраны, что делать дальше?! Следующий этап предполагает получение из них каких-либо общих характеристик, закономерностей, и их последующую интерпретацию. Исходный массив данных — не самый удобный объект для формулировки предположений, данные необходимо охарактеризовать набором удобно интерпретируемых с точки зрения статистики параметров. Таковыми, например, являются среднее и стандартное отклонение, вариация и др.

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

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

Одномерная задача с числовыми данными

Имеются выборочные данные \(x_1,x_2,\ldots, x_n\), где \(n\) — общее число наблюдений, \(x_i\) — вещественные числа. В качестве базовой описательной статистики для такого одномерного массива данных будем рассматривать следующий набор характеристик:

  1. среднее значение
  2. максимальное и минимальное значения
  3. стандартное отклонение
  4. коэффициент вариации
  5. квантили: 25%, 50%, 75%
  6. моду
  7. коэффициент асимметрии
  8. эксцесс
  9. ядерную оценку плотности распределения вероятностей

Создадим на Python функцию descriptive1d, принимающую в качестве аргумента массив данных \(x_1,x_2,\ldots, x_n\) и возвращающую 9 перечисленных выше параметров (или объектов — как в случае с ядерной оценкой плотности).

import numpy as np
import scipy.stats as st


def descriptive1d(x):
    _x = x  # Для возможности предобработки данных (например, исключения нечисловых значений) 
    result = []
    result.append(np.mean(_x))
    result.append((np.min(_x), np.max(_x)))
    result.append(np.std(_x))
    result.append(result[0]/result[-1])
    result.append((np.percentile(_x, 25), np.percentile(_x, 50), np.percentile(_x, 75)))
    result.append(st.mode(_x)) 
    result.append(st.skew(_x))  # асимметрия 
    result.append(st.kurtosis(_x))  # эксцесс
    _range = np.linspace(0.9 * np.min(_x), 1.1 * np.max(_x), 100)
    result.append(st.gaussian_kde(_x)(_range))  # оценка плотности распределения
    return tuple(result)

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

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

Квартили 25%, 50% (медиана) и 75% вычисляются c помощью стандартной функции percentile из пакета numpy. Она использует по умолчанию алгоритм линейной интерполяции, если значение процентиля не выпадает точно на целый индекс. Например, медиана, или 50%-процентиль для набора [1,2,3,4] будет 2.5 (здесь нецелый индекс тоже равен \(2.5 = \frac{1+4}{100}\cdot 50\))

Коэффициент асимметрии

Эксцесс — количественная характеристика формы распределения вблизи математического ожидания по отношению к нормальному распределению. Коэффициент эксцесса (\(\gamma_2\)) определяется как отношение оценки четвертого центрального момента к четвертой степени стандартного отклонения за вычетом 3: $$ \gamma_2 = \frac{m_4}{\sigma^4} - 3.$$ В случае положительного эксцесса имеет место "островершинность" распределения вблизи математического ожидания, в случае отрицательного его значения — распределение является более пологим, чем нормальное. Нулевое значение соответствует нормальному распределению (корректирующее слагаемое "\(-3\)" присутствует в формуле специально, чтобы обеспечить равенство эксцесса нулю для нормального распределения).

Знание непараметрической оценки плотности дает общее представление о распределении данных. Это в некотором роде обобщение понятия гистограммы. Нормированная, чтобы обеспечить единичную площадь под своим графиком гистограмма, представляет собой кусочно-линейную оценку плотности распределения вероятностей; ядерные оценки  являются гладкими "обобщениями" гистограммы и  дают более адекватное представление о плотности распределения. Более подробно о роли непараметрических оценок плотности можно прочесть в статье "Оценивание плотности распределения вероятностей по выборочным данным".

Гистограмма модельных данных

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

# формирование двухвершинного (двумодового) распределени из пары нормальных 
sample = np.hstack((1.5*np.random.rand(200).ravel(),2+1.2*np.random.rand(300)))

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

#coding: utf-8

from __future__ import print_function
import numpy as np
import scipy.stats as st
from pylab import *


# Нахождение базовых статистических показателей (описательная статистика)
def descriptive1d(x):
    _x = x  # Для возможности предобработки данных (например, исключения нечисловых значений) 
    result = []
    result.append(len(x)) # Чисо элементов выборки
    result.append(np.mean(_x)) # среднее
    result.append((np.min(_x), np.max(_x))) # (min, max)
    result.append(np.std(_x)) # стандартное отклонение
    result.append(100.0 * result[-1]/result[0]) # коэффициент вариации (Пирсона)
    result.append((np.percentile(_x, 25), np.percentile(_x, 50), np.percentile(_x, 75))) # квартили
    result.append(st.mode(_x))  # мода
    result.append(st.skew(_x))  # асимметрия 
    result.append(st.kurtosis(_x))  # эксцесс
    _range = np.linspace(0.9 * np.min(_x), 1.1 * np.max(_x), 100) # область определения для оценки плотности
    result.append((_range, st.gaussian_kde(_x)(_range)))  # оценка плотности распределения
    return tuple(result)



# формирование двухвершинного (двумодового) распределени из пары нормальных 
sample = np.hstack((1.5 * np.random.rand(200).ravel(), 2 + 1.2 * np.random.rand(300)))

# отрисовка гистограммы с 20-ю разбиениями интервала значений выборки
# bins - это число столбцов гистограммы, включая столбцы нулевой высоты
figure()
hist(sample, bins=20) 

# Вычисление важных показателей
n, m, minmax, s, cv, perct, mode, skew, kurt, kde = descriptive1d(sample)

print('Число элементов выборки: {0:d}'.format(n))
print('Среднее значение: {0:.4f}'.format(m))
print('Минимальное и максимальное значения: ({0:.4f}, {1:.4f})'.format(*minmax))
print('Стандартное отклонение: {0:.4f}'.format(s))
print('Коэффициент вариации (Пирсона): {0:.4f}'.format(cv))
print('Квартили: (25%) = {0:.4f}, (50%) = {1:.4f}, (75%) = {2:.4f}'.format(*perct))
print('Коэффициент асимметрии: {0:.4f}'.format(skew))
print('Коэффициент эксцесса: {0:.4f}'.format(kurt))

# Отрисовка оценки плотности распределения
figure()
plot(kde[0], kde[1])
title('pdf-estimation')
show()

А результаты его работы буду приблизительно такими:

Число элементов выборки: 500
Среднее значение: 1.8730
Минимальное и максимальное значения: (0.0039, 3.1957)
Стандартное отклонение: 1.0088
Коэффициент вариации (Пирсона): 0.2018
Квартили: (25%)  = 0.9050, (50%)  = 2.2307, (75%)  = 2.7351
Коэффициент асимметрии: -0.4418
Коэффициент эксцесса: -1.2886
Оценка плотности Python

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

Перейдем к дальнейшему усовершенствованию скрипта.

Нередки ситуации, когда даже в числовых данных имеются пропуски. Например, метеорологическая станция осуществляет измерения температуры воздуха каждые полчаса; в силу каких-то причин на ней происходит сбой, и в базу данных температуры записываются несуществующие показатели. Это могут быть либо специально принятые для этого обозначения, либо невероятные значения вообще (например, температура -999).  Если говорить о представлении числовых массивов в NumPy, для описания данных с пропусками или содержащих неправильные значения существует специальный тип данных (masked array - замаскированный массив).

Приведем простой пример использования masked array в NumPy.

>>> import numpy as np
>>> x=np.array([1,2,3,1,2,99])
>>> x_masked = np.ma.masked_values(x, 99)
>>> x.mean()
18.0
>>> x_masked.mean()
1.8
>>> 

В приведенном примере мы создали masked_array (x_masked), для которого указали, что значение 99 является неправильным и не должно учитываться при вычислениях (в частности, среднего значения). Отметим также, что специализированный пакет для обработки данных Pandas обладает еще более гибким аппаратом работы с неправильными (и\или пропущенными) данными. Фильтрация неправильных данных  может решаться множеством различных способов, а Pandas и NumPy представляют базовый функционал для организации такой работы.

Мы модифицируем скрипт нахождения базовых статистических показателей, чтобы он пропускал специальные нечисловые значения NumPy: np.nan, np.inf. Для этого, заблаговременно, нами уже была введена вспомогательная переменная _x; в этом месте и будем выполнять фильтрацию.

Для этого, достаточно внести следующие изменения в существующий код:

def descriptive1d(x):
    _x = np.array(x, dtype=np.float)  # Для возможности предобработки данных (например, исключения нечисловых значений) 
    _x = _x[~np.isnan(_x)*~np.isinf(_x)]
    result = []

Описание массива нечисловых данных

Работа с массивом нечисловых данных также может быть проведена средствами NumPy. Для этих целей в пакете представлен специальных тип элементов  numpy.object, позволяющих оперировать с такими данными.

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

  1. Определение количетсва уникальных представителей массива
  2. Определения частоты встречаемости этих представителей
  3. Определение наиболее часто встречающегося представителя (мода распределения)
  4. Определение наиболее редко встречающегося представителя
  5. Количественная характеристика неоднородности частот встречаемости элементов (например, энтропия по Шеннону)

Оформим вычисление показателей 1-5 (и нескольких еще дополнительных) в виде функции, реализованной на Python:

def descriptive1d_qualitative(x):
    from collections import Counter
    c = Counter(x) # используется для подсчета уникальных элементов
    result = []
    result.append(len(x))  # Число элементов в выборке
    result.append(dict(c))  # Число уникальных элементов представителей
    result.append({x: c[x]/float(result[0])  for x in dict(c)}) # Встречаемость элементов
    result.append(filter(lambda x: x[1]==c.most_common()[0][1], c.most_common())) # Наиболее часто встречающиеся элементы
    result.append(filter(lambda x: x[1]==c.most_common()[-1][1], c.most_common())) # Наиболее редко встречающиеся элементы
    result.append(st.entropy(np.array(c.values()))) # Энтропия по Шеннону 
    result.append(np.log(len(x))) # Максимальная энтропия
    return tuple(result)

Для тестирования создадим случайный массив нечисловых данных — букв латинского алфавита:

import string, random

data = map(lambda x: random.choice(string.ascii_lowercase), xrange(100))

n, uniques, freq, most_common, rare, entropy, maxentropy = descriptive1d_qualitative(data)

print('Число элементов выборки: {0:d}'.format(n))
print('Количество объектов по видам: {0}'.format(uniques))
print('Встречаемость объектов: {0}'.format(freq))
print('Часто встречающиеся: {0}'.format(most_common))
print('Редко встречающиеся: {0}'.format(rare))
print('Энтропия по Шеннону: {0:.4f}'.format(entropy))
print('Максимальная энтропия: {0:.4f}'.format(maxentropy))

В итоге, получим что-то подобное (вариации возможны из-за случайной выборки):

Число элементов выборки: 100
Количество объектов по видам: {'a': 4, 'c': 5, 'e': 6, 'd': 4, 'g': 3, 'f': 5, 'i': 1, 'h': 3, 'k': 4, 'j': 4, 'm': 2, 'l': 3, 'o': 2, 'n': 4, 'q': 6, 'p': 5, 's': 1, 'r': 6, 'u': 4, 't': 7, 'v': 8, 'y': 5, 'x': 4, 'z': 4}
Встречаемость объектов: {'a': 0.04, 'c': 0.05, 'e': 0.06, 'd': 0.04, 'g': 0.03, 'f': 0.05, 'i': 0.01, 'h': 0.03, 'k': 0.04, 'j': 0.04, 'm': 0.02, 'l': 0.03, 'o': 0.02, 'n': 0.04, 'q': 0.06, 'p': 0.05, 's': 0.01, 'r': 0.06, 'u': 0.04, 't': 0.07, 'v': 0.08, 'y': 0.05, 'x': 0.04, 'z': 0.04}
Часто встречающиеся: [('v', 8)]
Редко встречающиеся: [('i', 1), ('s', 1)]
Энтропия по Шеннону: 3.0880
Максимальная энтропия: 4.6052

Смешанные данные

Случай смешанных данных во многом определяется особенностью решаемой задачи. Нечисловые данные могут быть отфильтрованы из общего массива числовых данных с использованием, например, маскированных массивов (masked array, numpy). Таким образом, задача может быть сведена к анализу массива числовых и нечисловых данных. Это подход, не учитывающий специфики задачи обработки данных, вполне применим если необходимо составить общую их характеристику. 

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

def separate_data(x):
    from numbers import Number
    nums = filter(lambda x: isinstance(x, Number), x)
    nans = filter(lambda x: not isinstance(x, Number), x)
    return nums, nans

Функция separate_data позволяет разделить данные на числовые и нечисловые, т.е. ситуация смешанных данных сводится к рассмотренным выше случаям.

>>> x = [1, 2, 'a', 'b', 4, 3.43]
>>> separate_data(x)
([1, 2, 4, 3.43], ['a', 'b'])

Заметим, что функция разделения на числовые и нечисловые данные строку, представляющую десятичную запись числа (например,  '1'), отнесет к нечисловым данным. Можно модифицировать функцию разделения, чтобы она учитывала и такой формат представления числовых данных. Для этого удобней всего воспользоваться методом isdigit() для строкового класса в Python.
 

blog comments powered by Disqus