Numpy

Несколько лет назад, занимался изучением теории музыки, продавал и писал аудио-инструментал для аренды или заказов. Изначально, процесс явно творческий, но вскоре, мой интерес к коммерческой части превысил и возник вопрос: В каком же темпе создавать ритм музыки?.

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

BPM [в музыке] показатель, для определения скорости исполнения композиции, путём измерения количества тактовых долей в минуту.

1: Пролог

Устанавливаем Matplotlibи Pandas с необходимыми зависимостями через pip-менеджер в консоли/терминале.

Создаём директорию, а потом виртуальное окружение для проекта. После, подключаем библиотеки в IDE [в моём случае: PyCharm].

2: BPM

BPM будем вычислять через функцию Detect tempo в FL Studio и через сайт tunebat.com

3: DataSet

Начинаем создание DataSetа [выборки-коллекции данных] в Excel, для каждого жанра. Экспортируем в CSV-формат с настройками разделителя запятой. Следующие CSV-файлы создавал в IDE, так удобнее. Выборки перемещаем в директорию, где находится файл самой программы.

В первой строке CSV-файлов указываются параметры, которые разделяются запятыми. Следующие строки содержат уже значения этих параметров. При окончательной проверке, DataSet должен последовательно содержать данные: названия трека, BPM и год выхода композиции. Будем использовать информацию выборки в сто песен, для каждого жанра из выбранных пяти.

4: Rap построение точечной диаграммы и гистограммы

Выборка взята здесь: rollingstone.com/100-greatest-hip-hop-songs-of-all-time
Сам CSV-DataSet: github.com/Rap.csv

На основе информации DataSet'а, создаём точечную диаграмму [Scatter Plots] для изучения взаимосвязи между BPM и годом выпуска, а также для отображения концентраций при ранжировании данных.

from matplotlib import pyplot as plt                              # Первый каноничный импортimport pandas as pd                                                    # Второй каноничный импорт для обработки DataSet'аplt.style.use('fivethirtyeight')                                         # Назначаем стилистику визуализацииdata_set = pd.read_csv('Rap.csv')                               # Считываем данные SCV-файла с DataSet'омbpm = data_set['bpm']                                                  # Переменная, для параметра BPM в каждой строкеyear = data_set['year']                                                  # Переменная, для параметра "год релиза" в каждой строкеplt.scatter(                                                                     # Построение точечного графика и его настройкаbpm, year,                                                                   # Данные для осей x и yc=bpm,                                                                        # Привязка цвета к нужной осиs=bpm*1.5,                                                                  # Зависимость размера точкиcmap='gist_heat',                                                        # Цветовая карта графикаedgecolor='black',                                                       # Цвет контура точкиlinewidth=.7                                                                 # Толщина контура точки)bar = plt.colorbar(                                                          # Построение шкалы BPMorientation='horizontal',                                            # Ориентация шкалыshrink=0.8,                                                               # Масштаб шкалыextend='both',                                                           # Скос краёв шкалыextendfrac=.1                                                           # Угол скоса краёв)bar.set_label('Шкала ударов в минуту', fontsize=18)   # Подпись шкалыplt.title('Популярность скорости '                                  # Заголовок графика  'исполнения в Rap\'е ', fontsize=25)plt.xlabel('BPM', fontsize=18)                                         # Ось абсциссplt.ylabel('Год релиза', fontsize=18)                               # Ось ординатplt.tight_layout()                                                              # Настройка параметров подзаголовков в области отображенияplt.show()                                                                        # Вывод на экран

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

import pandas as pdfrom matplotlib import pyplot as pltfrom collections import Counterplt.style.use("fivethirtyeight")data_set = pd.read_csv('Rap end.csv')index = data_set['number']ranges = data_set['bpm_range']counter = Counter()for index in ranges:counter.update(index.split(';'))range_bpm = []value = []for item in counter.most_common(100):range_bpm.append(item[0])value.append(item[1])range_bpm.reverse()value.reverse()plt.barh(range_bpm, value,linewidth=.5,edgecolor='black',color='#e85b45',label='Количество точек на предыдущем графике')plt.legend()plt.title('Популярность интервала значений BPM в rap\'е', fontsize=25)plt.xlabel('Количество песен в диапазоне BPM', fontsize=18)plt.ylabel('Диапазоны BPM', fontsize=18)plt.tight_layout()plt.show()

5: Рок

Выборка взята здесь: rockfm.ru/top100
Сам CSV-DataSet: github.com/Rock.csv

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

from matplotlib import pyplot as plt                              # Первый каноничный импортimport pandas as pd                                                    # Второй каноничный импорт для обработки DataSet'аplt.style.use('fivethirtyeight')                                         # Назначаем стилистику визуализацииdata_set = pd.read_csv('Rock.csv')                             # Считываем данные SCV-файла с DataSet'омbpm = data_set['bpm']                                                  # Переменная, для параметра BPM в каждой строкеyear = data_set['year']                                                  # Переменная, для параметра "год релиза" в каждой строкеplt.scatter(                                                                     # Построение точечного графика и его настройкаbpm, year,                                                                   # Данные для осей x и yc=bpm,                                                                        # Привязка цвета к нужной осиs=bpm*1.5,                                                                  # Зависимость размера точкиcmap='gist_heat',                                                        # Цветовая карта графикаedgecolor='black',                                                       # Цвет контура точкиlinewidth=.7                                                                 # Толщина контура точкиalpha=.7                                                                      # Прозрачность точки)bar = plt.colorbar(                                                          # Построение шкалы BPMorientation='horizontal',                                            # Ориентация шкалыshrink=0.8,                                                               # Масштаб шкалыextend='both',                                                           # Скос краёв шкалыextendfrac=.1                                                           # Угол скоса краёв)bar.set_label('Шкала ударов в минуту', fontsize=18)    # Подпись шкалыplt.title('Популярность скорости '                                   # Заголовок графика  'исполнения в роке', fontsize=25)plt.xlabel('BPM', fontsize=18)                                          # Ось абсциссplt.ylabel('Год релиза', fontsize=18)                                # Ось ординатplt.tight_layout()                                                               # Настройка параметров подзаголовков в области отображенияplt.show()                                                                         # Вывод на экран

import pandas as pdfrom matplotlib import pyplot as pltfrom collections import Counterplt.style.use("fivethirtyeight")data_set = pd.read_csv('Rock end.csv')index = data_set['number']ranges = data_set['bpm_range']counter = Counter()for index in ranges:counter.update(index.split(';'))range_bpm = []value = []for item in counter.most_common(100):range_bpm.append(item[0])value.append(item[1])range_bpm.reverse()value.reverse()plt.barh(range_bpm, value,linewidth=.5,edgecolor='black',color='#e85b45',label='Количество точек на предыдущем графике')plt.legend()plt.title('Популярность интервала значений BPM в роке', fontsize=25)plt.xlabel('Количество песен в диапазоне BPM', fontsize=18)plt.ylabel('Диапазоны BPM', fontsize=18)plt.tight_layout()plt.show()

6: Блюз

Выборка взята здесь: digitaldreamdoor.com/best_bluesong
Сам CSV-DataSet: github.com/Blues.csv

from matplotlib import pyplot as plt                              # Первый каноничный импортimport pandas as pd                                                    # Второй каноничный импорт для обработки DataSet'аplt.style.use('fivethirtyeight')                                         # Назначаем стилистику визуализацииdata_set = pd.read_csv('Blues.csv')                            # Считываем данные SCV-файла с DataSet'омbpm = data_set['bpm']                                                  # Переменная, для параметра BPM в каждой строкеyear = data_set['year']                                                  # Переменная, для параметра "год релиза" в каждой строкеplt.scatter(                                                                     # Построение точечного графика и его настройкаbpm, year,                                                                   # Данные для осей x и yc=bpm,                                                                        # Привязка цвета к нужной осиs=bpm*1.5,                                                                  # Зависимость размера точкиcmap='gist_heat',                                                        # Цветовая карта графикаedgecolor='black',                                                       # Цвет контура точкиlinewidth=.7                                                                 # Толщина контура точки)bar = plt.colorbar(                                                          # Построение шкалы BPMorientation='horizontal',                                            # Ориентация шкалыshrink=0.8,                                                               # Масштаб шкалыextend='both',                                                           # Скос краёв шкалыextendfrac=.1                                                           # Угол скоса краёв)bar.set_label('Шкала ударов в минуту', fontsize=18)    # Подпись шкалыplt.title('Популярность скорости '                                   # Заголовок графика  'исполнения в блюзе', fontsize=25)plt.xlabel('BPM', fontsize=18)                                          # Ось абсциссplt.ylabel('Год релиза', fontsize=18)                                # Ось ординатplt.tight_layout()                                                               # Настройка параметров подзаголовков в области отображенияplt.show()                                                                         # Вывод на экран

import pandas as pdfrom matplotlib import pyplot as pltfrom collections import Counterplt.style.use("fivethirtyeight")data_set = pd.read_csv('Blues end.csv')index = data_set['number']ranges = data_set['bpm_range']counter = Counter()for index in ranges:counter.update(index.split(';'))range_bpm = []value = []for item in counter.most_common(100):range_bpm.append(item[0])value.append(item[1])range_bpm.reverse()value.reverse()plt.barh(range_bpm, value,linewidth=.5,edgecolor='black',color='#e85b45',label='Количество точек на предыдущем графике')plt.legend()plt.title('Популярность интервала значений BPM в блюзе', fontsize=25)plt.xlabel('Количество песен в диапазоне BPM', fontsize=18)plt.ylabel('Диапазоны BPM', fontsize=18)plt.tight_layout()plt.show()

7: Chillout

Выборка взята здесь: open.spotify.com
Сам CSV-DataSet: github.com/Chillout.csv

Много наложений точек друг на друга. К сожалению, не знаю, как это исправить. Пришлось сделать точки более прозрачными, с помощью аргумента alpha функции .scatter.

from matplotlib import pyplot as plt                              # Первый каноничный импортimport pandas as pd                                                    # Второй каноничный импорт для обработки DataSet'аplt.style.use('fivethirtyeight')                                         # Назначаем стилистику визуализацииdata_set = pd.read_csv('Chillout.csv')                         # Считываем данные SCV-файла с DataSet'омbpm = data_set['bpm']                                                  # Переменная, для параметра BPM в каждой строкеyear = data_set['year']                                                  # Переменная, для параметра "год релиза" в каждой строкеplt.scatter(                                                                     # Построение точечного графика и его настройкаbpm, year,                                                                   # Данные для осей x и yc=bpm,                                                                        # Привязка цвета к нужной осиs=bpm*1.5,                                                                  # Зависимость размера точкиcmap='gist_heat',                                                        # Цветовая карта графикаedgecolor='black',                                                       # Цвет контура точкиlinewidth=.7                                                                 # Толщина контура точкиalpha=.5                                                                      # Прозрачность точки)bar = plt.colorbar(                                                          # Построение шкалы BPMorientation='horizontal',                                            # Ориентация шкалыshrink=0.8,                                                               # Масштаб шкалыextend='both',                                                           # Скос краёв шкалыextendfrac=.1                                                           # Угол скоса краёв)bar.set_label('Шкала ударов в минуту', fontsize=18)   # Подпись шкалыplt.title('Популярность скорости '                                  # Заголовок графика  'исполнения в Chillout', fontsize=25)plt.xlabel('BPM', fontsize=18)                                          # Ось абсциссplt.ylabel('Год релиза', fontsize=18)                               # Ось ординатplt.tight_layout()                                                               # Настройка параметров подзаголовков в области отображенияplt.show()                                                        # Вывод на экран

import pandas as pdfrom matplotlib import pyplot as pltfrom collections import Counterplt.style.use("fivethirtyeight")data_set = pd.read_csv('Chillout end.csv')index = data_set['number']ranges = data_set['bpm_range']counter = Counter()for index in ranges:counter.update(index.split(';'))range_bpm = []value = []for item in counter.most_common(100):range_bpm.append(item[0])value.append(item[1])range_bpm.reverse()value.reverse()plt.barh(range_bpm, value,linewidth=.5,edgecolor='black',color='#e85b45',label='Количество точек на предыдущем графике')plt.legend()plt.title('Популярность интервала значений BPM в Chillout', fontsize=25)plt.xlabel('Количество песен в диапазоне BPM', fontsize=18)plt.ylabel('Диапазоны BPM', fontsize=18)plt.tight_layout()plt.show()

8: EDM

Выборка взята здесь: edmcharts.net
Сам CSV-DataSet: github.com/EDM.csv

Здесь также для наглядности пришлось сделать точки ещё более прозрачными. Если кто-то знает, как исправить дефект наложения, прошу написать в комментариях.

from matplotlib import pyplot as plt                              # Первый каноничный импортimport pandas as pd                                                    # Второй каноничный импорт для обработки DataSet'аplt.style.use('fivethirtyeight')                                         # Назначаем стилистику визуализацииdata_set = pd.read_csv('EDM.csv')                             # Считываем данные SCV-файла с DataSet'омbpm = data_set['bpm']                                                  # Переменная, для параметра BPM в каждой строкеyear = data_set['year']                                                  # Переменная, для параметра "год релиза" в каждой строкеplt.scatter(                                                                     # Построение точечного графика и его настройкаbpm, year,                                                                   # Данные для осей x и yc=bpm,                                                                        # Привязка цвета к нужной осиs=bpm*1.5,                                                                  # Зависимость размера точкиcmap='gist_heat',                                                        # Цветовая карта графикаedgecolor='black',                                                       # Цвет контура точкиlinewidth=.7                                                                 # Толщина контура точкиalpha=.2                                                                      # Прозрачность точки)bar = plt.colorbar(                                                          # Построение шкалы BPMorientation='horizontal',                                            # Ориентация шкалыshrink=0.8,                                                               # Масштаб шкалыextend='both',                                                           # Скос краёв шкалыextendfrac=.1                                                           # Угол скоса краёв)bar.set_label('Шкала ударов в минуту', fontsize=18)   # Подпись шкалыplt.title('Популярность скорости '                                  # Заголовок графика  'исполнения в EDM', fontsize=25)plt.xlabel('BPM', fontsize=18)                                          # Ось абсциссplt.ylabel('Год релиза', fontsize=18)                               # Ось ординатplt.tight_layout()                                                               # Настройка параметров подзаголовков в области отображенияplt.show()                                                        # Вывод на экран

import pandas as pdfrom matplotlib import pyplot as pltfrom collections import Counterplt.style.use("fivethirtyeight")data_set = pd.read_csv('EDM end.csv')index = data_set['number']ranges = data_set['bpm_range']counter = Counter()for index in ranges:counter.update(index.split(';'))range_bpm = []value = []for item in counter.most_common(100):range_bpm.append(item[0])value.append(item[1])range_bpm.reverse()value.reverse()plt.barh(range_bpm, value,linewidth=.5,edgecolor='black',color='#e85b45',label='Количество точек на предыдущем графике')plt.legend()plt.title('Популярность интервала значений BPM в EDM', fontsize=25)plt.xlabel('Количество песен в диапазоне BPM', fontsize=18)plt.ylabel('Диапазоны BPM', fontsize=18)plt.tight_layout()plt.show()

9: Заключение

Самым простым графиком сравним количество попаданий в каждый диапазон, композиций, из всех проанализированных ранее жанров*.

* такие жанры как ethnic, ambient, folk, dubstep, reggae и др, не удалось к сожалению разобрать из-за отсутствия качественной выборки...

from matplotlib import pyplot as pltplt.style.use('fivethirtyeight')x = ['<60', '60-80', '80-100', '100-120', '120-140', '140-160', '160-180']y = [2, 32, 144, 119, 135, 39, 29]plt.plot(x, y, label='BPM', c='#e85b45')plt.legend()plt.title('Сравнение всех диапазонов BPM во всех жанрах', fontsize=25)plt.xlabel('Диапазон BPM', fontsize=18)plt.ylabel('Количество треков', fontsize=18)plt.tight_layout()plt.show()

Серия из 5 постов для начинающих представляет собой ремикс первой главы книги 2015 года под названием Clojure для науки о данных (Clojure for Data Science). Автор книги, Генри Гарнер, любезно дал согласие на использование материалов книги для данного ремикса с использованием языка Python.

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

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

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

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

Пост 1 посвящен подготовке среды и данных.

Статистика

Важно не кто голосует, а кто подсчитывает голоса

Иосиф Сталин

Как только перед нами возникает задача проанализировать данные, которые состоят из двух и более чисел, становится содержательным вопрос о том, каким образом эти числа распределены. Вы, наверное, уже слыхали такие выражения, как длинный хвост и правило 80/20. Они касаются разброса чисел по диапазону. В этой главе мы продемонстрируем смысл распределений и познакомим с наиболее полезным из них: нормальным распределением.

При изучении распределений чрезвычайную важность играет наглядная и удобная визуализация данных, и для этого мы воспользуемся Python-овской библиотекой pandas. Мы покажем, как пользоваться ею для загрузки, преобразования и разведывательного анализа реальных данных, а также начнем работать с фундаментальной библиотекой numpy для научных вычислений. Мы проведем сопоставительный анализ результатов двух общенациональных выборов всеобщих выборов в Великобритании 2010 г. и российских выборов депутатов Государственной Думы Федерального Собрания РФ шестого созыва 2011 г. и увидим, каким образом даже элементарный анализ может предъявить подтверждающие данные о потенциальных фальсификациях.

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

В этом посте мы будем пользоваться тремя главными библиотеками экосистемы SciPy: одноименной библиотекой SciPy для выполнения сложных математико-статистических расчетов, библиотекой pandas для загрузки данных из разнообразных источников, управления ими и их визуализации, а также библиотекой NumPy в основном для работы с массивами и матрицами.

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

В основе библиотеки pandas лежит понятие кадра данных, DataFrame, т.е. структуры, состоящей из строк и столбцов, или записей и полей. Если у вас есть опыт работы с реляционными базами данных, то таблицы pandas можно представить, как таблицы базы данных. Каждый столбец в кадре данных поименован, а каждая строка имеет одинаковое число столбцов, как и любая другая. Загрузить данные в кадр данных pandas можно несколькими способами, и тот, которым мы воспользуемся, будет зависеть от того, в каком виде наши данные хранятся:

• Если данные представлены текстовым файлом с разделением полей данных запятыми (.csv) или символами табуляции (.tsv), то мы будем использовать функцию чтения данных read_csv

• Если данные представлены файлом Excel (например, файл .xls или .xlsx), то мы воспользуемся функцией чтения данных read_excel

• Для любого другого источника данных (внешняя база данных, веб-сайт, буфер обмена данными, JSON-файлы, HTML-файлы и т. д.) предусмотрен ряд других функций

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

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

pd.read_excel('data/ch01/UK2010.xls')

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

def load_uk(): '''Загрузить данные по Великобритании''' return pd.read_excel('data/ch01/UK2010.xls') 

Эта функция вернет кадр данных DataFrame библиотеки pandas, содержащий данные по Великобритании. Далее в этой главе, мы определим дополнительные имплементации загрузки этого же и еще одного набора данных.

Первая строка электронной таблицы UK2010.xls содержит имена столбцов. Функция библиотеки pandas read_excel резервирует их в качестве имен столбцов возвращаемого кадра данных. Начнем обследование данных с их проверки атрибут кадра данных columns возвращает имена столбцов в виде списка, при этом адресация атрибутов осуществляется при помощи оператора точки (.):

def ex_1_1(): '''Получить имена полей кадра данных''' return load_uk().columns

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

Index(['Press Association Reference', 'Constituency Name', 'Region', 'Election Year', 'Electorate', 'Votes', 'AC', 'AD', 'AGS', 'APNI', ... 'UKIP', 'UPS', 'UV', 'VCCA', 'Vote', 'Wessex Reg', 'WRP', 'You', 'Youth', 'YRDPL'], dtype='object', length=144)

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

• Информация для Ассоциации прессы: число, идентифицирующее избирательный округ (представленный одним депутатом)

• Название избирательного округа: стандартное название, данное избирательному округу

• Регион: географический район Великобритании, где округ расположен

• Год выборов: год, в котором выборы состоялись

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

• Голосование: общее число проголосовавших

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

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

def ex_1_2(): '''Получить значения поля "Год выборов"''' return load_uk()['Election Year']

В результате будет выведен следующий список:

0 2010.01 2010.02 2010.0...646 2010.0647 2010.0648 2010.0649 2010.0650 NaNName: Election Year, dtype: float64

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

def ex_1_3(): '''Получить значения в поле "Год выборов" без дубликатов''' return load_uk()['Election Year'].unique()

[ 2010. nan]

Значение 2010 еще больше подкрепляет наши ожидания в отношении того, что эти данные относятся к 2010 году. Впрочем, наличие специального значения nan, от англ. not a number, т.е. не число, которое сигнализирует о пропущенных данных, является неожиданным и может свидетельствовать о проблеме с данными.

Мы еще не знаем, в скольких элементах набора данных пропущены значения, и установив их число, мы смогли бы решить, что делать дальше. Простой способ подсчитать такие пустые элементы состоит в использовании подкласса словарей Counter языка Python из модуля collections. Этот словарь трансформирует последовательность значений в коллекцию, где ключам поставлены в соответствие количества появлений элементов данных, т.е. их частоты:

def ex_1_4(): '''Рассчитать частоты в поле "Год выборов"  (количества появлений разных значений)''' return Counter( load_uk()['Election Year'] )

Counter({nan: 1, 2010.0: 650})

Нам не потребуется много времени, чтобы получить подтверждение, что в 2010 г. в Великобритании было 650 избирательных округов. Знание предметной области, как в этом случае, имеет неоценимое значение при проверке достоверности новых данных. Таким образом, весьма вероятно, что значение nan является посторонним, и его можно удалить. Мы увидим, как это сделать, в следующем разделе.

Исправление данных

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

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

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

def ex_1_5(): '''Вернуть отфильтрованную по полю "Год выборов"  запись в кадре данных (в виде словаря)''' df = load_uk() return df[ df['Election Year'].isnull() ]

Выражение dt['Election Year'].isnull() вернет булеву последовательность, в которой все элементы, кроме последнего, равны False, в результате чего будет возвращена последняя запись кадра данных. Если Вы знаете язык запросов SQL, то отметите, что этот метод очень похож на условный оператор WHERE.

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

 df = load_uk() return df[ df[ 'Election Year' ].notnull() ]

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

def load_uk_scrubbed(): '''Загрузить и отфильтровать данные по Великобритании''' df = load_uk() return df[ df[ 'Election Year' ].notnull() ]

Теперь при написании любого фрагмента кода с описанием доступа к набору данных в файле, можно выбирать вариант загрузки: обычный при помощи функции load_uk либо его аналог с очисткой данных load_uk_scrubbed. Приведенный выше пример должен вернуть список из 650 чисел, обозначающих количество избирателей в каждом избирательном округе Великобритании.

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

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

Пост 3 для начинающих посвящен генерированию распределений, их свойствам, а также графикам для их сопоставительного анализа.

Булочник и Пуанкаре

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

В те времена хлебопекарное ремесло регламентировалось государством, и Пуанкаре обнаружил, что, хотя результаты взвешивания буханок хлеба подчинялись нормальному распределению, пик находился не на публично афишируемом 1 кг, а на 950 г. Он сообщил властям о булочнике, у которого он регулярно покупал хлеб, и тот был оштрафован. Такова легенда ;-).

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

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

Генерирование распределений

В целях развития нашего интуитивного понимания относительно нормального распределения и дисперсии, давайте смоделируем честного и нечестного булочников, и для этого воспользуемся функцией генерирования нормально распределенных случайных величин stats.norm.rvs. (rvs от англ. normal variates, т.е. случайные величины). Честного булочника можно смоделировать в виде нормального распределения со средним значением 1000, что соответствует справедливой буханке хлеба весом 1 кг. При этом мы допустим наличие дисперсии в процессе выпекания, которая приводит к стандартному отклонению в 30г.

def honest_baker(mu, sigma): '''Модель честного булочника''' return pd.Series( stats.norm.rvs(loc, scale, size=10000) )def ex_1_18(): '''Смоделировать честного булочника на гистограмме''' honest_baker(1000, 30).hist(bins=25) plt.xlabel('Честный булочник')  plt.ylabel('Частота') plt.show()

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

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

def dishonest_baker(mu, sigma): '''Модель нечестного булочника''' xs = stats.norm.rvs(loc, scale, size=10000)  return pd.Series( map(max, bootstrap(xs, 13)) ) def ex_1_19(): '''Смоделировать нечестного булочника на гистограмме''' dishonest_baker(950, 30).hist(bins=25) plt.xlabel('Нечестный булочник')  plt.ylabel('Частота') plt.show()

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

Совершенно очевидно, что эта гистограмма выглядит не совсем так, как другие, которые мы видели. Среднее значение по-прежнему равно 1 кг, но разброс значений вокруг среднего больше не является симметричным. Мы говорим, что эта гистограмма показывает смещенное нормальное распределение.

Асимметрия

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

Библиотека pandas располагает функцией skew для измерения асимметрии:

def ex_1_20(): '''Получить коэффициент асимметрии нормального распределения''' s = dishonest_baker(950, 30) return { 'среднее' : s.mean(),  'медиана' : s.median(),  'асимметрия': s.skew() }

{'асимметрия': 0.4202176889083849, 'медиана': 998.7670301469957, 'среднее': 1000.059263920949}

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

Графики нормального распределения

Ранее в этой главе мы познакомились с квантилями как средством описания статистического распределения данных. Напомним, что функция quantile принимает число между 0 и 1 и возвращает значение последовательности в этой точке. 0.5-квантиль соответствует значению медианы.

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

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

def qqplot( xs ): '''Квантильный график (график квантиль-квантиль, Q-Q plot)''' d = {0:sorted(stats.norm.rvs(loc=0, scale=1, size=len(xs))), 1:sorted(xs)} pd.DataFrame(d).plot.scatter(0, 1, s=5, grid=True) df.plot.scatter(0, 1, s=5, grid=True) plt.xlabel('Квантили теоретического нормального распределения') plt.ylabel('Квантили данных') plt.title ('Квантильный график', fontweight='semibold')def ex_1_21(): '''Показать квантильные графики  для честного и нечестного булочников''' qqplot( honest_baker(1000, 30) ) plt.show() qqplot( dishonest_baker(950, 30) ) plt.show()

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

Выше показан квантильный график для честного булочника. Далее идет квантильный график для нечестного булочника:

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

^{Надписи: нормально распределенные, тяжелые хвосты, легкие хвосты, скошенность влево, скошенность вправо, раздельные кластеры}

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

Технические приемы сопоставительной визуализации

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

Коробчатые диаграммы

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

def ex_1_22(): '''Показать коробчатую диаграмму с данными честного и нечестного булочников''' d = {'Честный булочник' :honest_baker(1000, 30), 'Нечестный булочник':dishonest_baker(950, 30)}  pd.DataFrame(d).boxplot(sym='o', whis=1.95, showmeans=True) plt.ylabel('Вес буханки (гр.)') plt.show()

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

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

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

Интегральные функции распределения

Интегральные функции распределения (ИФР), также именуемые кумулятивными функциями распределения, от англ. Cumulative Distribution Function (CDF), описывают вероятность, что значение, взятое из распределения, будет меньше x. Как и все распределения вероятностей, их значения лежат в диапазоне между 0 и 1, где 0 это невозможность, а 1 полная определенность. Например, представьте, что я собираюсь бросить шестигранный кубик. Какова вероятность, что выпадет значение меньше 6?

Для уравновешенного кубика вероятность выпадения пятерки или меньшего значения равна 5/6. И наоборот, вероятность, что выпадет единица, равна всего1/6. Тройка или меньше соответствуют равным шансам то есть вероятности 50%.

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

^{ИФР и квантили тесно друг с другом связаны ИФР является инверсией квантильной функции. Если 0.5-квантиль соответствует значению 1000, тогда ИФР для 1000 составляет 0.5.}

Подобно тому, как функция pandas quantile позволяет нам отбирать значения из распределения в конкретных точках, эмпирическая ИФР empirical_cdf позволяет нам внести значение из последовательности и вернуть значение в диапазоне между 0 и 1. Это функция более высокого порядка, т.е. она принимает значение (в данном случае последовательность значений) и возвращает функцию, которую потом можно вызывать, сколько угодно, с различными значениями на входе, и возвращая ИФР для каждого из них.

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

Построим график ИФР одновременно для честного и нечестного булочников. Для этих целей можно воспользоваться функцией библиотеки pandas построения двумерного графика plot для визуализации ИФР, изобразив на графике исходные данные то есть выборки из распределений честного и нечестного булочников в сопоставлении с вероятностями, вычисленными относительно эмпирической ИФР. Функция plot ожидает, что значения x и значения y будут переданы в виде двух раздельных последовательностей значений. Для этих целей мы воспользуемся конструктором кадра данных pandas DataFrame.

Чтобы изобразить оба распределения на одном графике, мы должны передать функции plot несколько серий. Для многих своих графиков pandas предоставляет функции, которые позволяют добавлять дополнительные серии. В случае с функцией plot мы можем присвоить указатель на создаваемый график, присвоив временной переменной (ax) результат первого вызова функции plot, и затем при повторных вызовах указывать эту переменную в именованном аргументе функции (ax=ax). Можно также передать необязательную метку серии. Мы выполним это в следующем ниже примере, чтобы на готовом графике отличить две серии друг от друга. Сначала определим универсальную функцию построения эмпирической ИФР против теоретической, которая получает на вход кортеж из двух серий (tp[1] и tp[3]) и их названий и метки осей, и затем вызовем ее:

def empirical_cdf(x): """Вернуть эмпирическую ИФР для x""" sx = sorted(x) return pd.DataFrame( {0: sx, 1:sp.arange(len(sx))/len(sx)} )def ex_1_23(): '''Показать графики эмпирической ИФР честного булочника в сопоставлении с нечестным''' df = empirical_cdf(honest_baker(1000, 30)) df2 = empirical_cdf(dishonest_baker(950, 30)) ax = df.plot(0, 1, label='Честный булочник')  df2.plot(0, 1, label='Нечестный булочник', grid=True, ax=ax)  plt.xlabel('Вес буханки') plt.ylabel('Вероятность') plt.legend(loc='best') plt.show()

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

Несмотря на то, что этот график выглядит совсем по-другому, он в сущности показывает ту же самую информацию, что и коробчатая диаграмма. Мы видим, что две линии пересекаются примерно в медиане 0.5, соответствующей 1000 гр. Линия нечестного булочника обрезается в нижнем хвосте и удлиняется на верхнем хвосте, что соответствует асимметричному распределению.

Следующая часть, часть 4, серии постов Python, наука о данных и выборы посвящена техническим приемам визуализации данных.

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

Описательные статистики

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

def ex_1_6(): '''Число значений в поле "Электорат"''' return load_uk_scrubbed()['Electorate'].count()

650

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

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

• Среднее значение

Наиболее распространенный способ усреднить набор данных взять его среднее значение. Среднее значение на самом деле представляет собой один из нескольких способов измерения центра распределения данных. Среднее значение числового ряда вычисляется на Python следующим образом:

def mean(xs):  '''Среднее значение числового ряда''' return sum(xs) / len(xs) 

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

def ex_1_7(): '''Вернуть среднее значение поля "Электорат"''' return mean( load_uk_scrubbed()['Electorate'] )

70149.94

На самом деле, библиотека pandas уже содержит функцию mean, которая гораздо эффективнее вычисляет среднее значение последовательности. В нашем случае ее можно применить следующим образом:

load_uk_scrubbed()['Electorate'].mean()

• Медиана

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

def median(xs): '''Медиана числового ряда''' n = len(xs) mid = n // 2 if n % 2 == 1: return sorted(xs)[mid] else: return mean( sorted(xs)[mid-1:][:2] )

Медианное значение электората Великобритании составляет:

def ex_1_8(): '''Вернуть медиану поля "Электорат"''' return median( load_uk_scrubbed()['Electorate'] )

70813.5

Библиотека pandas тоже располагает встроенной функцией для вычисления медианного значения, которая так и называется median.

• Дисперсия

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

Она может содержать целые числа от 1 до 99 или два ряда чисел, состоящих из 49 нулей и 50 девяносто-девяток, а может быть и так, что она содержит ряд из 98 чисел, равных -1 и одно единственное значение 5048, или же вообще все значения могут быть равны 50.

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

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

Выражение

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

def variance(xs): '''Дисперсия числового ряда, несмещенная дисперсия при n <= 30''' mu = mean(xs) n = len(xs) n = n-1 if n in range(1, 30) else n  square_deviation = lambda x : (x - mu) ** 2  return sum( map(square_deviation, xs) ) / n

Для вычисления квадрата выражения используется оператор языка Python возведения в степень **.

Поскольку мы взяли средний квадрат отклонения, т.е. получили квадрат отклонения и затем его среднее, то единицы измерения дисперсии тоже будут в квадрате, т.е. дисперсия электората Великобритании будет измеряться "людьми в квадрате". Несколько неестественно рассуждать об избирателях в таком виде. Единицу измерения можно привести к более естественному виду, снова обозначающему "людей", путем извлечения квадратного корня из дисперсии. В результате получим так называемое стандартное отклонение, или среднеквадратичное отклонение:

def standard_deviation(xs): '''Стандартное отклонение числового ряда''' return sp.sqrt( variance(xs) ) def ex_1_9(): '''Стандартное отклонение поля "Электорат"''' return standard_deviation( load_uk_scrubbed()['Electorate'] )

7672.77

В библиотеке pandas функции для вычисления дисперсии и стандартного отклонения имплементированы соответственно, как var и std. При этом последняя по умолчанию вычисляет несмещенное значение, поэтому, чтобы получить тот же самый результат, нужно применить именованный аргумент ddof=0, который сообщает, что требуется вычислить смещенное значение стандартного отклонения:

load_uk_scrubbed()['Electorate'].std( ddof=0 )

• Квантили

Медиана это один из способов вычислить из списка срединное значение, т.е. находящееся ровно по середине, дисперсия же предоставляет способ измерить разброс данных вокруг среднего значения. Если весь разброс данных представить на шкале от 0 до 1, то значение 0.5 будет медианным.

Для примера рассмотрим следующую ниже последовательность чисел:

[10 11 15 21 22.5 28 30]

Отсортированная последовательность состоит из семи чисел, поэтому медианой является число 21 четвертое в ряду. Его также называют 0.5-квантилем. Мы можем получить более полную картину последовательности чисел, взглянув на 0.0 (нулевой), 0.25, 0.5, 0.75 и 1.0 квантили. Все вместе эти цифры не только показывают медиану, но также обобщают диапазон данных и сообщат о характере распределения чисел внутри него. Они иногда упоминаются в связи с пятичисловой сводкой.

Один из способов составления пятичисловой сводки для данных об электорате Великобритании показан ниже. Квантили можно вычислить непосредственно в pandas при помощи функции quantile. Последовательность требующихся квантилей передается в виде списка.

def ex_1_10(): '''Вычислить квантили: возвращает значение в последовательности xs,  соответствующее p-ому проценту''' q = [0, 1/4, 1/2, 3/4, 1] return load_uk_scrubbed()['Electorate'].quantile(q=q)

0.00 21780.000.25 65929.250.50 70813.500.75 74948.501.00 109922.00Name: Electorate, dtype: float64

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

Группирование данных в корзины

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

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

На приведенном выше рисунке показано 15 значений x, разбитых на 5 равноразмерных корзин. Подсчитав количество точек, попадающих в каждую корзину, мы можем четко увидеть, что большинство точек попадают в корзину по середине, а меньшинство в корзины по краям. Следующая ниже функция Python nbin позволяет добиться того же самого результата:

def nbin(n, xs):  '''Разбивка данных на частотные корзины''' min_x, max_x = min(xs), max(xs) range_x = max_x - min_x fn = lambda x: min( int((abs(x) - min_x) / range_x * n), n-1 ) return map(fn, xs)

Например, мы можем разбить диапазон 0-14 на 5 корзин следующим образом:

list( nbin(5, range(15)) )

[0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4]

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

def ex_1_11(): '''Разбиmь электорат Великобритании на 5 корзин''' series = load_uk_scrubbed()['Electorate'] return Counter( nbin(5, series) )

Counter({2: 450, 3: 171, 1: 26, 0: 2, 4: 1})

Количество точек в крайних корзинах (0 и 4) значительно ниже, чем в корзинах в середине количества, судя по всему, растут по направлению к медиане, а затем снова снижаются. В следующем разделе мы займемся визуализацией формы этих количеств.

Гистограммы

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

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

def ex_1_12(): '''Построить гистограмму частотных корзин        электората Великобритании''' load_uk_scrubbed()['Electorate'].hist() plt.xlabel('Электорат Великобритании') plt.ylabel('Частота') plt.show()

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

Число корзин, на которые данные разбиваются, можно сконфигурировать, передав в функцию при построении гистограммы именованный аргумент bins:

def ex_1_13(): '''Построить гистограмму частотных корзин  электората Великобритании с 200 корзинами''' load_uk_scrubbed()['Electorate'].hist(bins=200) plt.xlabel('Электорат Великобритании') plt.ylabel('Частота') plt.show()

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

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

def ex_1_14(): '''Построить гистограмму частотных корзин  электората Великобритании с 20 корзинами''' load_uk_scrubbed()['Electorate'].hist(bins=20) plt.xlabel('Электорат Великобритании') plt.ylabel('Частота') plt.show()

Ниже показана гистограмма теперь уже из 20 корзин:

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

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

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

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

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

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

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

Центральная предельная теорема

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

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

def ex_1_15(): '''Показать гистограмму равномерного распределения  синтетического набора данных''' xs = stats.uniform.rvs(0, 1, 10000) pd.Series(xs).hist(bins=20) plt.xlabel('Равномерное распределение') plt.ylabel('Частота') plt.show()

Обратите внимание, что в этом примере мы впервые использовали тип Series библиотеки pandas для числового ряда данных.

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

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

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

def bootstrap(xs, n, replace=True):  '''Вернуть список массивов меньших размеров  по n элементов каждый''' return np.random.choice(xs, (len(xs), n), replace=replace) def ex_1_16(): '''Построить гистограмму средних значений''' xs = stats.uniform.rvs(loc=0, scale=1, size=10000) pd.Series( map(sp.mean, bootstrap(xs, 10)) ).hist(bins=20) plt.xlabel('Распределение средних значений')  plt.ylabel('Частота') plt.show()

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

Хотя величина среднего значения близкая к 0 или 1 не является невозможной, она является чрезвычайно невероятной и становится менее вероятной по мере роста числа усредненных чисел и числа выборочных средних. Фактически, на выходе получается результат очень близкий к нормальному распределению.

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

До 20-ого века самого термина еще не существовало, хотя этот эффект был зафиксирован еще в 1733 г. французским математиком Абрахамом де Mуавром (Abraham de Moivre), который использовал нормальное распределение, чтобы аппроксимировать число орлов в результате бросания уравновешенной монеты. Исход бросков монеты лучше всего моделировать при помощи биномиального распределения, с которым мы познакомимся в главе 4, Классификация. В отличие от центральной предельной теоремы, которая позволяет получать выборки из приближенно нормального распределения, библиотека ScyPy содержит функции для эффективного генерирования выборок из самых разнообразных статистических распределений, включая нормальное:

def ex_1_17(): '''Показать гистограмму нормального распределения  синтетического набора данных''' xs = stats.norm.rvs(loc=0, scale=1, size=10000) pd.Series(xs).hist(bins=20) plt.xlabel('Нормальное распределение') plt.ylabel('Частота') plt.show()

Отметим, что в функции sp.random.normal параметр loc это среднее значение, scale дисперсия и size размер выборки. Приведенный выше пример сгенерирует следующую гистограмму нормального распределения:

По умолчанию среднее значение и стандартное отклонение для получения нормального распределения равны соответственно 0 и 1.

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

Заключительный пост 5 для начинающих посвящен сопоставительной визуализации электоральных данных.

Сопоставительная визуализация электоральных данных

Теперь рассмотрим набор данных других всеобщих выборов, на этот раз Российских, проходивших в 2011 г. Россия гораздо более крупная страна, и поэтому данные о проголосовавших на выборах там гораздо объемнее. Для этого мы загрузим в оперативную память один большой TSV-файл с разделением полей данных символом табуляции.

def load_ru(): '''Загрузить данные по России''' return pd.read_csv('data/ch01/Russia2011.tsv', '\t')

Посмотрим, какие имена столбцов имеются в российских данных:

def ex_1_29(): '''Показать список полей электоральных  данных по России''' return load_ru().columns

Будет выведен следующий список столбцов:

Index(['Код ОИК', 'ОИК ', 'Имя участка','Число избирателей, внесенных в список избирателей',...'Политическая партия СПРАВЕДЛИВАЯ РОССИЯ','Политическая партия ЛДПР - Либерально-демократическая партия России','Политическая партия "ПАТРИОТ РОССИИ"','Политическая партия КОММУНИСТИЧЕСКАЯ ПАРТИЯ КОММУНИСТ РОССИИ','Политическая партия "Российская объединенная демократическая партия "ЯБЛОКО"','Политическая партия "ЕДИНАЯ РОССИЯ"','Всероссийская политическая партия "ПАРТИЯ РОСТА"'],dtype='object')

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

Наряду с набором данных функция библиотеки pandas rename ожидает словарь, в котором ключам с текущими именами столбцов поставлены в соответствие значения с новыми именами. Если объединить ее с данными, которые мы уже рассматривали, то мы получим следующее:

def load_ru_victors(): '''Загрузить данные по России,  выбрать, переименовать и вычислить поля''' new_cols_dict = { 'Число избирателей, внесенных в список избирателей':'Электорат', 'Число действительных избирательных бюллетеней': 'Действительные бюллетени', 'Политическая партия "ЕДИНАЯ РОССИЯ"':'Победитель'  } newcols = list(new_cols_dict.values())  df = load_ru().rename( columns=new_cols_dict )[newcols]  df['Доля победителя'] = df['Победитель'] / df['Действительные бюллетени']  df['Явка'] = df['Действительные бюллетени'] / df['Электорат']  return df

Библиотека pandas располагает функцией безопасного деления divide, которая идентична операции /, но защищает от деления на ноль. Она вместо пропущенного значения (nan) в одном из полей подставляет значение, передаваемое в именованном аргументе fill_value. Если же оба значения поля равны nan, то результат будет отсутствовать. Поэтому операцию деления можно было бы переписать следующим образом:

 df[ 'Доля победителя' ] = \ df[ 'Победитель' ].divide( df[ 'Действительные бюллетени' ], \ fill_value=1 )

Визуализация электоральных данных РФ

Мы ранее видели, что гистограмма явки на выборы в Великобритании была приближенно нормально распределенной (хотя и с легкими хвостами). Теперь, когда мы загрузили и преобразовали данные о выборах в России, посмотрим, насколько они сопоставимы:

def ex_1_30(): '''Показать гистограмму  электоральных данных по России''' load_ru_victors()['Явка'].hist(bins=20) plt.xlabel('Явка в России')  plt.ylabel('Частота') plt.show()

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

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

Учитывая ожидания, заданные данными из Британии и центральной предельной теоремой (ЦПТ), такой результат любопытен. Для начала покажем данные на квантильном графике:

def ex_1_31(): '''Показать квантильный график  победителя на выборах в РФ''' qqplot( load_ru_victors()['Доля победителя'].dropna() ) plt.show()

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

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

На самом деле, этот квантильный график дезориентирует, и происходит этот именно потому, что хвост очень тяжелый: плотность точек между 0.5 и 1.0 на гистограмме говорит о том, что пик должен составлять порядка 0.7 с последующим правым хвостом за пределами 1.0. Наличие значения, превышающего 100% явно выходит за рамки логики, но квантильный график не объясняет это (он не учитывает, что речь идет о процентах), так что внезапное отсутствие данных за пределами 1.0 интерпретируется как подрезанный правый хвост.

С учетом центральной предельной теоремы и того, что мы наблюдали в данных выборов в Великобритании, тенденция к 100% явке избирателей на выборы выглядит очень любопытно. Давайте выполним параллельный сопоставительный анализ наборов данных по Великобритании и России.

Сравнительная визуализация

Предположим, мы хотели бы сравнить распределение электоральных данных между Великобританией и Россией. Мы уже видели в этой главе, как использовать ИФР и коробчатые диаграммы, поэтому теперь займемся исследованием альтернативного варианта, который аналогичен гистограмме.

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

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

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

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

Функции массы вероятности

Функция массы вероятности (ФМВ), от англ. Probability Mass Function (PMF), чаще именуемая функцией вероятности дискретной случайной величины, имеет много общего с гистограммой. Однако, вместо того, чтобы показывать количества значений, попадающих в группы, она показывает вероятность, что взятое из распределения число будет в точности равно заданному значению. Поскольку функция закрепляет вероятность за каждым значением, которое может быть возвращено распределением, и поскольку вероятности измеряются по шкале от 0 до 1, (где 1 соответствует полной определенности), то площадь под функцией массы вероятности равна 1.

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

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

def plot_as_pmf(dt, label, ax): '''График функции вероятности дискретной случайной величины (или функции массы вероятности)''' s = pd.cut(dt, bins=40, labels=False) # разбить на 40 корзин pmf = s.value_counts().sort_index() / len(s) # подсчитать кво в корзинах newax = pmf.plot(label=label, grid=True, ax=ax)  return newax

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

def ex_1_32(): '''Сопоставление данных явки по Великобритании и РФ, данные нормализованы на основе функции массы вероятностей''' ax = plot_as_pmf(load_uk_victors()['Явка'], 'Великобритания', None) plot_as_pmf(load_ru_victors()['Явка'], 'Россия', ax) plt.xlabel('Интервальные группы явки') # Частотные корзины plt.ylabel('Вероятность') plt.legend(loc='best') plt.show()

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

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

Данные российских выборов показывают чрезвычайно аномальный результат, хотя и не настолько высокий, как модальный пик в центре распределения, который приблизительно соответствует 50% явке. Исследователь Питер Климек (Peter Klimek) и его коллеги в Венском медицинском университете пошли дальше и предположили, что этот результат является явным признаком подтасовки результатов голосования.

Диаграммы рассеяния

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

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

Заключительный технический прием визуализации, с которым мы познакомим в этой главе, представлен диаграммой рассеяния. Диаграммы рассеяния очень хорошо подходят для визуализации взаимосвязей между двумя переменными: там, где существует линейная взаимосвязь, на графике она будет видна, как диагональная направленность. Библиотека pandas содержит для этого вида графиков функцию scatter с такими же аргументами, что и для функции двумерных графиков plot.

def ex_1_33(): '''Показать диаграмму рассеяния  выборов в Великобритании''' df = load_uk_victors()[ ['Явка', 'Доля победителей'] ] df.plot.scatter(0, 1, s=3) plt.xlabel('Явка') plt.ylabel('Доля победителя') plt.show()

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

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

^{Как отмечалось ранее, британские выборы 2010 г. были далеко необычными: они привели к "подвисшему" парламенту и коалиционному правительству. Фактически, "победители" в данном случае представлены обеими сторонами, которые были противниками, вплоть до дня выборов. И поэтому голосование за любую из партий считается как голосование за победителя.}

Затем, мы создадим такую же диаграмму рассеяния для выборов в России:

def ex_1_34(): '''Показать диаграмму рассеяния выборов в РФ''' df = load_ru_victors()[ ['Явка', 'Доля победителя'] ] df.plot.scatter(0, 1, s=3) plt.xlabel('Явка') plt.ylabel('Доля победителя') plt.show()

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

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

Настройка прозрачности рассеяния

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

Выполнить настройку альфа-канала, регулирующего прозрачность изображаемых на графике pandas точек можно при помощи именованного аргумента alpha в функции scatter в виде числа между 0 и 1, где 1 означает полную непрозрачность, 0 полную прозрачность.

def ex_1_35(): '''Показать диаграмму рассеяния (с прозрачностью) выборов в РФ''' df = load_ru_victors()[ ['Явка', 'Доля победителя'] ] rows = sp.random.choice(df.index.values, 10000) df.loc[rows].plot.scatter(0, 1, s=3, alpha=0.1) plt.xlabel('Явка') plt.ylabel('Доля победителя') plt.axis([0, 1.05, 0, 1.05]) plt.show()

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

Приведенная выше диаграмма рассеяния показывает общую направленность совместного изменения доли победителя и явки на выборы. Мы видим корреляцию между двумя значениями и "горячую точку" в правом верхнем углу графика, которая соответствует явке близкой к 100% и 100%-ому голосованию в пользу побеждающей стороны. Как раз эта особенность в частности является признаком того, что исследователи из Венского медицинского университета обозначили как сигнатура фальсификации выборов. Этот факт также подтверждается результатами других спорных выборов по всему миру, например, таких как президентские выборы 2011 г. в Уганде.

^{Результаты многих других выборов по всему миру на уровне округов доступны здесь. На указанном веб-сайте можно получить ссылку на научно-исследовательскую работу и скачать другие наборы данных, на которых можно попрактиковать то, что Вы узнали в этой главе об исправлении и преобразовании реальных данных.}

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

Выводы

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

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

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

Пост 4 для начинающих посвящен техническим приемам визуализации данных.

Важность визуализации

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

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

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

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

На самом деле статистические распределения для последовательностей, взятых из крупных выборок, могут быть настолько надежными, что даже незначительное отклонение от них может являться свидетельством противоправной деятельности. Закон Бенфорда, или закон первой цифры, показывает любопытную особенность случайных чисел, которые генерируются в широком диапазоне. Единица появляется в качестве ведущей цифры примерно в 30% случаев, в то время как цифры крупнее появляется все реже и реже. Например, девятка появляется в виде первой цифры менее чем в 5% случаев.

^{Закон Бенфорда назван в честь физика Фрэнка Бенфорда (Frank Benford), который сформулировал его в 1938 г., показав его состоятельность на различных источниках данных. Проявление этого закона было ранее отмечено американским астрономом Саймоном Ньюкомом (Simon Newcomb), который еще более 50 лет назад до него обратил внимание на страницы своих логарифмических справочников: страницы с номерами, начинавшихся с цифры 1, имели более потрепанный вид.}

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

Визуализация данных об электорате

Вернемся к данным выборов и сравним электоральную последовательность, которую мы создали ранее, относительно теоретической нормальной ИФР. Для создания нормальной ИФР из последовательности значений можно воспользоваться функцией sp.random.normal библиотеки SciPy, как уже было показано выше. Среднее значение и стандартное отклонение по умолчанию равны соответственно 0 и 1, поэтому нам нужно предоставить измеренные среднее значение и стандартное отклонение, взятые из электоральных данных. Эти значения для наших электоральных данных составляют соответственно 70150 и 7679.

Ранее в этой главе мы уже генерировали эмпирическую ИФР. Следующий ниже пример просто сгенерирует обе ИФР и выведет их на одном двумерном графике:

def ex_1_24(): '''Показать эмпирическую и подогнанную ИФР  электората Великобритании''' emp = load_uk_scrubbed()['Electorate'] fitted = stats.norm.rvs(emp.mean(), emp.std(ddof=0), len(emp)) df = empirical_cdf(emp) df2 = empirical_cdf(fitted) ax = df.plot(0, 1, label='эмпирическая')  df2.plot(0, 1, label='подогнанная', grid=True, ax=ax)  plt.xlabel('Электорат') plt.ylabel('Вероятность') plt.legend(loc='best') plt.show()

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

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

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

def ex_1_25(): '''Показать квантильный график  электората Великобритании''' qqplot( load_uk_scrubbed()['Electorate'] ) plt.show()

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

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

Добавление производных столбцов

В целях выяснения процента электората, который проголосовал за одну из двух партий, требуется вычислить сумму голосов, отданных за каждую из них. Для этого нам понадобится создать новое поле данных Victors (Победители) из данных, которые соответствуют Консервативной (Con) и Либерально-демократической (LD) партиям и заодно проверим, имеются ли пропущенные значения.

def ex_1_26(): '''Вычислить производное поле данных "Победители" и  число имеющихся в нем пропущенных значений''' df = load_uk_scrubbed() df['Победители'] = df['Con'] + df['LD'] freq = Counter(df['Con'].apply( lambda x: x > 0 )) print('Поле "Победители": %d, в т.ч. пропущено %d'  % (freq[True], freq[False]))

Поле "Победители": 631, в т.ч. пропущено 19

Результат показывает, что в 19 случаях данные отсутствуют. Очевидно, что в каком-то из столбцов: столбце Con либо столбце LD (либо обоих), данные отсутствуют, но в каком именно? Снова воспользуемся словарем Counter, чтобы увидеть масштаб проблемы:

'''Проверить пропущенные значения в полях "Консервативная партия" (Con) и   "Либерально-демократическая партия" (LD)'''df = load_uk_scrubbed()Counter(df['Con'].apply(lambda x: x > 0)),  Counter(df['LD'].apply(lambda x: x > 0))

(Counter({False: 19, True: 631}), Counter({False: 19, True: 631}))

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

def ex_1_27(): '''Выборка полей данных по условию, что поля "Консервативная партия" (Con) и  "Либерально-демократическая" (LD) не пустые''' df = load_uk_scrubbed() rule = df['Con'].isnull() & df['LD'].isnull() return df[rule][['Region', 'Electorate', 'Con', 'LD']]

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

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

def load_uk_victors(): '''Загрузить данные по Великобритании,  выбрать поля и отфильтровать''' df = load_uk_scrubbed() rule = df['Con'].notnull() df = df[rule][['Con', 'LD', 'Votes', 'Electorate']]  df['Победители'] = df['Con'] + df['LD']  df['Доля победителей'] = df['Победители'] / df['Votes']  df['Явка'] = df['Votes'] / df['Electorate'] return df

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

def ex_1_28(): '''Показать квантильный график победителей  на выборах в Великобритании''' qqplot( load_uk_victors()['Доля победителей'] ) plt.show()

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

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

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

Следующая заключительная часть, часть 5, серии постов Python, наука о данных и выборы посвящена сопоставительной визуализации электоральных данных.

Удивительное дело, но в русскоязычном сегменте интернета почти нет материала, разъясняющего понятным языком соглашение Эйнштейна о суммировании. Не менее удивительно то, что материалов, позволяющих понять принцип работы функции einsum в русскоязычном интернете ещё меньше. На английском есть довольно развёрнутый ответ о работе einsum на stack overflow, а на русском только некоторое число сайтов, предоставляющих кривой перевод этого самого ответа. Хочу исправить эту проблему с недостатком материалов, и всех, кому интересно приглашаю к прочтению!

Обсуждаем соглашение Эйнштейна

Прежде всего отмечу, что соглашение Эйнштейна чаще всего используется в тензорном анализе и его приложениях, поэтому дальше в статье будет несколько референсов к тензорам.
Когда вы только начинаете работать с тензорами, вас может смутить, что кроме привычных подстрочных индексов, используются также и надстрочные индексы, которые по началу вообще можно принять за возведение в степень. Пример:
"а с верхним индексом i" будет записано как , а "a в квадрате с верхним индексом i" будет записываться . Возможно, по-началу это вводит в заблуждение и кажется неудобным, но со временем можно привыкнуть.

Соглашение: далее в статье объекты вида или я буду называть термами.

О чём вообще соглашение Эйнштейна?

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

Правило 1: Суммирование ведётся по всем индексам, повторяющимся дважды в одном терме.

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

С использованием соглашения Эйнштейна это выражение может быть переписано так:

Таким образом мы избавляемся от знака суммы, и просто пишем единственный терм. Обратим внимание, что в этом терме индекс i повторяется дважды, а значит, в соответствие с первым правилом мы понимаем, что суммирование ведётся по индексу i, а точнее, по всем возможным значениям, которые принимает этот индекс.

Рассмотрим ещё один пример: пусть нам нужно умножить матрицу на вектор . Результатом будет являться вектор . По определению:

Соглашение Эйнштейна позволяет избавиться от знака суммы:

Заметим, что в терм индекс i входит один раз, а индекс j входит два раза, а значит, суммирование ведётся по индексу j.

Определение 1. Индекс, который входит в терм дважды, называется фиктивным индексом.

Определение 2. Свободным индексом назовём все индексы в терме, не являющие фиктивными.

Отметим, что каждый фиктивный индекс может быть заменён любым другим фиктивным индексом, при условии, что

1. Новый фиктивный индекс не входит в множество свободных индексов терма.

2. Новый фиктивный индекс принимает то же множество значений, что и старый фиктивный индекс.

Чтобы объяснить проще, рассмотрим следующий код на языке Python:

for i in range(M):    for j in range(N):        b[i, j] = A[i, j] * v[j]

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

Правило 2. В каждом терме не может встречаться более двух одинаковых индексов.

Второе правило говорит нам, что мы можем написать , но не можем написать или , несмотря на то, что на практике такие выражения всё же имеют смысл.
Больше примеров:

здесь i является фиктивным индексом, т.к. повторяется дважды;

здесь i является свободным индексом, а j фиктивным;

здесь и i, и j являются фиктивными индексами;

здесь и i, и j являются фиктивными индексами;

не правильно по второму правилу (индекс i входит в терм трижды);

Из примеров выше можно заключить, что когда мы считаем число вхождений индексов в терм, мы не делаем разницы между верхними и нижними индексами, и считаем их вместе. Ещё один важный пример: когда мы видим выражение следующего вида

Мы должны понимать, что это выражение записано верно, и не противоречит второму правилу. Действительно, если посчитать все вхождения индексов, то получится, что индекс i входит 3 раза, как и индекс j, но в выражении записано два терма, а не один, и если посчитать вхождение индексов в каждый терм отдельно (как того и требует второе правило), то мы увидим, что ничего не нарушается.

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

Рассмотрим несколько примеров для закрепления этого правила:

этот пример мы уже рассматривали выше, здесь i является свободным индексом левой части уравнения, и свободным индексом правой части уравнения;

пример посложнее. Посчитаем вхождения индексов для каждого терма: в первый терм правой части k и j входят дважды, значит, они являются фиктивными индексами, i входит один раз, значит, является свободным. Во второй терм правой части k входит два раза, i один, значит, k фиктивный, i свободный. В левой части индекс i входит один раз, а значит, является свободным. Итог: индекс i является свободным для обеих частей уравнения, а значит, правило 3 выполнено.

Рассмотрим так же несколько примеров, в которых третье правило не выполняется:

слева i является свободным индексом, но справа свободны индексы i и j;

слева свободен индекс j, но справа свободен индекс i. Свободные индексы не совпадают;

здесь слева свободен индекс i, а справа свободны индексы i, j;

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

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

Там сейчас не очень важно, что из себя представляется каждая , и что такое . Наша задача сейчас исключительно синтаксическая игра. Нужно упростить выражение, особо не вникая в смысл происходящего.
Прежде всего видно, что свободными индексами являются , а фиктивными, соответственно индексы . Расположим индексы в соседних множителях так, чтобы в первом множителе индекс, по которому идёт суммирование, стоял снизу, а во втором тот же самый индекс стоял сверху. Так же заметим, что множителями являются тензоры , и у них в верхнем регистре уже стоит . Чтобы повысить читаемость, будем оборачивать множители в скобки, и только потом ставить индексы. Само же упрощённое выражение переписывается из исходного почти дословно:

Ура, мы научились упрощать сложные выражения с помощью соглашения Эйнштейна!

Обсуждаем einsum

einsum это функция, присутствующая в нескольких популярных библиотеках для Python (NumPy, TensorFlow, PyTorch). Во всех библиотеках, в которых эта функция реализована, она работает одинаково (с точностью до функционала структур, определённых в конкретной библиотеке), поэтому нет смысла рассматривать один и тот же пример в разных библиотеках, достаточно рассказать про einsum в одной конкретной библиотеке. Далее в статье я буду использовать NumPy. einsum применяет соглашение Эйнштейна о суммировании к переданным массивам. Функция принимает множество опционально аргументов, про них лучше почитать в документации, мы же сейчас разберём, как передавать шаблон, по которому функция будет применять соглашение Эйнштейна.

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

Теперь пусть мы хотим перемножить их программно. Ну, это можно довольно просто сделать с помощью трёх вложенных циклов:

M = np.zeros((3, 2))for i in range(3):    for j in range(2):        for k in range(5):            M[i, j] = A[i, k] * B[k, j]

Либо, используя функцию einsum можно написать это произведение в одну строчку:

M = np.einsum("ik,kj->ij", A, B)

Разберёмся, что за магия происходит в этой строчке. einsum принимает один обязательный аргумент: шаблон, по которому будет применено соглашение Эйнштейна. Шаблон этот выглядит так:

"{индексы, определяющие размерность первого массива},{индексы, определяющие размерность второго массива}->{индексы, определяющие размерность результирующего массива}"

Поведение шаблона einsum определяется следующими правилами:

• Если один и тот же индекс встречается слева и справа от запятой (до стрелочки), то суммирование будет вестись по этому индексу;

• Если после стрелочки ничего не написано, то суммирование произойдёт по всем встреченным осям;

• Никакой индекс не должен встречаться 3 и более раз;

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

Рассмотрим теперь некоторое количество примеров разной степени сложности, чтобы закрепить понимание einsum:

Одна einsum, чтобы править всеми

Пример 1. Сумма всех значений вектора:

vector = np.array([1, 2, 3, 4, 5])result = np.einsum("i->", vector)print(result)

Пример 2. Сумма всех значений матрицы:

matrix = np.array([[1, 2], [3, 4], [5, 6]])result = np.einsum("ij->", matrix)print(result)

Пример 3. Сумма значений по столбцам:

matrix = np.array([[1, 2], [3, 4], [5, 6]])result = np.einsum("ij->j", matrix)print(result)

Пример 4. Сумма значений по строкам:

matrix = np.array([[1, 2], [3, 4], [5, 6]])result = np.einsum("ij->i", matrix)print(result)

Пример 5. Транспонирование (я об этом не написал, но оси, по которым суммирования не произошло, мы можем возвращать в любом порядке):

matrix = np.array([[1, 2], [3, 4], [5, 6]])result = np.einsum("ij->ji", matrix)print(result)

Пример 6. Умножение матрицы на вектор:

matrix = np.array([[1, 2], [3, 4], [5, 6]])vector = np.array([[1, 2]])result = np.einsum("ij,kj->ik", matrix, vector)print(result)

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

Пример 7. Умножение матрицы на матрицу:

matrix1 = np.array([[1, 2], [3, 4], [5, 6]])matrix2 = np.array([[1, 0], [0, 1]])result = np.einsum("ik,kj->ij", matrix1, matrix2)print(result)

Пример 8. Скалярное произведение векторов:

vector1 = np.array([[1, 2, 3]])vector2 = np.array([[1, 1, 1]])result = np.einsum("ik,jk->", vector1, vector2)print(result)

Пример 9. След матрицы:

matrix1 = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])result = np.einsum("ii->", matrix1)print(result)

Пример 10. Адамарово (покомпонентное) произведение:

matrix1 = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])matrix2 = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])result = np.einsum("ij,ij->ij", matrix1, matrix2)print(result)

Это может показаться контринтутивно, но, как написано выше: если не понятно, что делает einsum запиши через циклы:

result = np.zeros(matrix1.shape, dtype="int32")for i in range(result.shape[0]):    for j in range(result.shape[1]):        result[i, j] += matrix1[i, j] * matrix2[i, j]print(result)

Пример 11. Кронекерово (внешнее) произведение векторов:

vector1 = np.array([1, 2, 3])vector2 = np.array([1, 0, 0])result = np.einsum("i,j->ij", vector1, vector2)print(result)

Пример 12. Транспонирование тензора:

A = np.array([[[0, 1], [1, 2], [2, 3]], [[1, 2], [2, 3], [3, 4]], [[2, 3], [3, 4], [4, 5]]])result = np.einsum("ijk->jki", A)print(result)

Пример 13. Произведение тензора на матрицу по третьей моде:

A = np.array([[[0, 1], [1, 2], [2, 3]], [[1, 2], [2, 3], [3, 4]], [[2, 3], [3, 4], [4, 5]]])U = np.array([[1, 2], [2, 3]])result = np.einsum("ijk,nk->ijn", A, U)print(result)

Итоги

Конечно, einsum поставляет только дополнительный синтаксический сахар. Всегда можно использовать цепочки вложенных циклов, множество библиотечных функций (np.dot, np.outer, np.tensordot, np.transpose, np.cumsum и т.д.), и вообще не использовать einsum. Но если потратить время и понять, как она работает, то можно научиться писать гораздо более сжатый, и, не побоюсь этого слова, эффективный код.

Ссылки

Ролик с примерами einsum (ещё больше примеров).

Соглашение Эйнштейна (база)

Соглашение Эйнштейна (продвинутая часть)

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

Да, это кажется безумием, но более привычные методы преобразования картинок работают в 1,5-2,5 раза медленнее (если нужен не read-only объект). Сегодня я покопаюсь в кишках обеих библиотек, расскажу почему так получилось и кто виноват. А также покажу финальный результат, который работает так же, только быстрее. Никаких репозиториев или пакетов не будет, только рассказ и рабочий код в конце. Но давайте обо всём по порядку.

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

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

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

Для разнообразия сегодня я буду запускать бенчмарки на Raspberry Pi 4 1800 MHz под 64-разрядной Raspberry Pi OS. В конце концов, где ещё может понадобиться компьютерное зрение, как не на Малинке :-)

Как работает преобразование в NumPy

Существует два общепринятых способа конвертировать изображение Pillow в NumPy, с равной вероятностью вы нагуглите один из них:

1. numpy.array(im) делает копию из изображения в массив NumPy.

2. numpy.asarray(im) то же самое, что numpy.array(im, copy=False), то есть якобы не делает копию, а использует память оригинального объекта. На самом деле всё несколько сложнее.

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

In [1]: from PIL import ImageIn [2]: import numpyIn [3]: im = Image.open('./canyon.jpg').resize((4096, 4096))In [4]: n = numpy.asarray(im)In [5]: n[:, :, 0] = 255ValueError: assignment destination is read-onlyIn [6]: n.flagsOut[6]:   C_CONTIGUOUS : True  F_CONTIGUOUS : False  OWNDATA : False  WRITEABLE : False  ALIGNED : True  WRITEBACKIFCOPY : False  UPDATEIFCOPY : False

Это сильно отличается от того, что будет, если использовать функцию numpy.array():

In [7]: n = numpy.array(im)In [8]: n[:, :, 0] = 255In [9]: n.flagsOut[9]:   C_CONTIGUOUS : True  F_CONTIGUOUS : False  OWNDATA : True  WRITEABLE : True  ALIGNED : True  WRITEBACKIFCOPY : False  UPDATEIFCOPY : False

При этом, если провести измерение, функция asarray() действительноработает значительно быстрее:

In [10]: %timeit -n 10 n = numpy.array(im)257 ms  1.27 ms per loop (mean  std. dev. of 7 runs, 10 loops each)In [11]: %timeit -n 10 n = numpy.asarray(im)179 ms  786 s per loop (mean  std. dev. of 7 runs, 10 loops each)

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

Интерфейс массивов NumPy

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

    @property    def __array_interface__(self):        shape, typestr = _conv_type_shape(self)        return {            "shape": shape,            "typestr": typestr,            "version": 3,            "data": self.tobytes(),        }

_conv_type_shape() описывает тип и размер массива, который должен получиться. А всё самое интересное происходит в методе tobytes(). Если проверить, сколько этот метод выполняется, станет понятно, что в общем-то NumPy от себя ничего не добавляет:

In [12]: %timeit -n 10 n = im.tobytes()179 ms  1.27 ms per loop (mean  std. dev. of 7 runs, 10 loops each)

Время точно совпадает с временем функции asarray(). Кажется виновник найден, осталось заменить вызов этой функции или ускорить её, и дело в шляпе, верно? Ну, не всё так просто.

Внутреннее устройство памяти в Pillow и NumPy

Устройство массивов в NumPy описывается чрезвычайно просто это непрерывный кусок памяти, начинающийся с определенного указателя. Плюс есть смещения (strides), которые задаются отдельно по каждому измерению.

В Pillow всё устроено принципиально иначе. Изображение хранится чанками, в каждом чанке находится целое количество строк изображения. Каждый пиксель занимает 1 или 4 байта (не от 1 до 4, а ровно). Соответственно, для каких-то режимов изображения какие-то байты не используются. Например, для RGB не используется последний байт в каждом пикселе, а для черно-белых изображений с альфа-каналом (режим LA) не используются два средних байта для того, чтобы альфа-канал был в последнем байте пикселя.

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

Я думаю, теперь понятно, для чего нужен метод tobytes() он переводит внутреннее представление изображения Pillow в непрерывный поток байтов одним куском без пропусков: как раз такое, какое может использовать NumPy. NumPy уже получая на вход объект bytes, может либо сделать копию, либо использовать его в режиме read-only. Тут я не уверен, сделано ли это, чтобы нельзя было обойти неизменность объектов bytesв Python, или есть какие-то реальные ограничения на уровне C API. Но, например, если на вход вместо bytes податьbytearray, то массив не будет read-only.

Но давайте всё же посмотрим на упрощенную версию tobytes():

    def tobytes(self):        self.load()        # unpack data        e = Image._getencoder(self.mode, "raw", self.mode)        e.setimage(self.im)        data, bufsize, s = [], 65536, 0        while not s:            l, s, d = e.encode(bufsize)            data.append(d)        if s < 0:            raise RuntimeError(f"encoder error {s} in tobytes")        return b"".join(data)

Тут видно, что создается "raw"энкодер и из него получаются чанки изображения не менее 65 килобайт памяти. Это и есть первое копирование: к концу функций у нас всё изображение в виде небольших чанков лежит в массиве data. Последней строкой происходит второе копирование: все чанки собираются в одну большую байтовую строку.

Кто виноват и что делать

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

Первое, что хочется отметить: отказываться от энкодера не вариант. Кто знает, какие детали реализации он от нас срывает. Переносить это всё на уровень Python или переписывать часть на C последнее дело.

Кажется, намного разумнее было бы в tobytes()заранее выделить буфер нужного размера, и уже в него записывать чанки. Но очевидно, что интерфейс энкодера так не работает: он уже возвращает чанки упакованные в объекты bytes. Тем не менее, если эти чанки не складировать, а сразу копировать в буфер, эти данные не будут вымываться из L2 кэша и быстро попадут куда надо. Что-то вроде такого:

def to_mem(im):    im.load()    e = Image._getencoder(im.mode, "raw", im.mode)    e.setimage(im.im)    mem = ... # we don't know yet    bufsize, offset, s = 65536, 0, 0    while not s:        l, s, d = e.encode(bufsize)        mem[offset:offset + len(d)] = d        offset += len(d)    if s < 0:        raise RuntimeError(f"encoder error {s} in tobytes")    return mem

Что же будет вместо mem. В идеале это должен быть массив NumPy. Создать его не представляет проблем, мы уже видели какие у него будут параметры в __array_interface__:

In [13]: shape, typestr = Image._conv_type_shape(im)In [14]: data = numpy.empty(shape, dtype=numpy.dtype(typestr))

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

In [15]: mem = data.reshape((data.size,))In [16]: mem[0:4] = b'abcd'ValueError: invalid literal for int() with base 10: b'abcd'

В данном случае кажется странным, что нельзя в массив байтов по срезу поместить байты. Но не забывайте, что, во-первых, слева могут быть не только байты, а во-вторых, библиотека называется NumPy, то есть работает с числами. К счастью, NumPy дает доступ и к непосредственной памяти массива прямо из Python. Это свойство data:

In [17]: data.dataOut[17]: <memory at 0x7f78854d68>In [18]: data.data[0] = 255NotImplementedError: sub-views are not implementedIn [19]: data.data.shapeOut[19]: (4096, 4096, 3)In [20]: data.data[0, 0, 0] = 255

Там находится объект memoryview. Вот только этот memoryviewкакой-то странный: он тоже многомерный, как и сам массив NumPy, ещё у него такой же тип объектов, как у самого массива. К счастью, это легко исправляется методом cast:

In [21]: mem = data.data.cast('B', (data.data.nbytes,))In [22]: mem.nbytes == mem.shape[0]Out[22]: TrueIn [23]: mem[0], mem[1]Out[23]: (255, 0)In [24]: mem[0:4] = b'1234'In [25]: mem[0], mem[1]Out[25]: (49, 50)

Складываем пазл вместе:

def to_numpy(im):    im.load()    # unpack data    e = Image._getencoder(im.mode, 'raw', im.mode)    e.setimage(im.im)    # NumPy buffer for the result    shape, typestr = Image._conv_type_shape(im)    data = numpy.empty(shape, dtype=numpy.dtype(typestr))    mem = data.data.cast('B', (data.data.nbytes,))    bufsize, s, offset = 65536, 0, 0    while not s:        l, s, d = e.encode(bufsize)        mem[offset:offset + len(d)] = d        offset += len(d)    if s < 0:        raise RuntimeError("encoder error %d in tobytes" % s)    return data

Проверяем:

In [26]: n = to_numpy(im)In [27]: numpy.all(n == numpy.array(im))Out[27]: TrueIn [28]: n.flagsOut[28]:   C_CONTIGUOUS : True  F_CONTIGUOUS : False  OWNDATA : True  WRITEABLE : True  ALIGNED : True  WRITEBACKIFCOPY : False  UPDATEIFCOPY : FalseIn [29]: %timeit -n 10 n = to_numpy(im)101 ms  260 s per loop (mean  std. dev. of 7 runs, 10 loops each)

Круто! Имеем ускорение в 2,5 раза с тем же функционалом и меньшее количество аллокаций.

Бенчмарки

В статье я взял достаточно большую картинку для тестов. Нет, дело не в том, что to_numpy()не дает ускорения на меньших размерах (ещё как даёт!). Дело в том, что в общем случае очень сложно добиться какого-то постоянного времени работы, когда дело касается выделения памяти. Аллокатор может затребовать новую память у системы, а может и старую сохранить. Может решить заполнить её нулями, а может и так отдать. В этом смысле работа с большими массивами хотя бы дает стабильный результат: мы всегда получаем худший случай.

Код:

In [30]: for i in range(6, 0, -1):    ...:     i = 128 * 2 ** i    ...:     print(f'\n\nSize: {i}x{i}   \t{i*i // 1024} KPx')    ...:     im = Image.new('RGB', (i, i))    ...:     print('\tnumpy.array()')    ...:     %timeit n = numpy.array(im)    ...:     print('\tnumpy.asarray()')    ...:     %timeit n = numpy.asarray(im)    ...:     print('\tto_numpy()')    ...:     %timeit n = to_numpy(im)    ...:     im = None    ...: 

Результаты:

Итого, получилось избавиться от лишней аллокации памяти, ускорить работу от 1,5 до 2,5 раз, попутно немного разобраться как NumPy работает с памятью.

Привет, Хабр!

Сегодня хочу поделиться своим небольшим домашним проектом:

ModulationPy (GiHub)

- модуль для моделирования цифровых схем модуляции (это которые PSK, QAM и т.п.). Проект был вдохновлен другой питоновской библиотекой: CommPy; однако, в рассмотренном классе задач с ней удалось даже немного посоревноваться!

На данный момент доступны два класса схем модуляции:

• M-PSK: Phase Shift Keying (фазовая цифровая модуляция)

• M-QAM: Quadratured Amplitude Modulation (квадратурная амплитудная модуляция)

  где M - это порядок модуляции.

Интересен модуль может быть, скорее всего, в разрезе образовательных целей в сфере беспроводной связи (подбор модуляций исходил именно из нее), однако, вдруг кому-то пригодится и для научных изысканий. Не MatLab'ом насущным едины!

Примеры использования

Итак, для начала скачиваем библиотеку с PyPI:

$ pip install ModulationPy

Либо устанавливаем из исходников, но на PyPI на момент написания статьи все-таки актуальная версия.

Зависимости две:

• numpy>=1.7.1

• matplotlib>=2.2.2 (для построения сигнальных созвездий)

Продемонстрируем на деле.

Итак, например, мы хотим использовать QPSK с поворотом фазы pi/4, двоичным входом и наложением по Грею (Gray mapping). Импортируем и инициализируем:

import numpy as npfrom ModulationPy import PSKModemmodem = PSKModem(4, np.pi/4,                 gray_map=True,                 bin_input=True)

Проверяем то ли мы инициализировали, нарисовав сигнальное созвездие методом plot_const():

modem.plot_const()

То же самое сделаем для 16-QAM (но для десятичных чисел на входе; указывать фазовый сдвиг не нужно - подразумевается наиболее распространенная прямоугольная QAM):

from ModulationPy import QAMModemmodem = QAMModem(16,                 gray_map=True,                  bin_input=False)modem.plot_const()

На данный момент модуляция QAM реализована по примеру функции qammod в Octave [4]. И, да, реализованы только четные (в том смысле, что результат log2(M) - четное число) схемы модуляции (4-QAM, 16-QAM, 64-QAM). Пусть и не совсем полный набор, но как бы то ни было, в популярных стандартах беспроводной связи все равно нет "нечетных" схем модуляции (насколько я знаю).

Далее предлагаю перейти, собственно, к главному в модемах: к модуляции и демодуляции. Для этого нам понадобятся два метода:modulate() и demodulate() , доступные в обоих классах.

Метод modulate() принимает на вход всего один аргумент:

• вектор входных значений (1-D ndarray of ints) - либо единиц и нулей, если выбрана опция bin_input=True , либо целых десятичных чисел от 0 до M-1, если bin_input=False ;

Методdemodulate() ожидает максимум два аргумента:

• вектор, который должен быть демодулирован (1-D ndarray of complex symbols) ;

• значение дисперсии аддитивного шума (float, по умолчанию 1.0).

Например, вот как это будет выглядеть для QPSK (двоичный вход/выход):

import numpy as npfrom ModulationPy import PSKModemmodem = PSKModem(4, np.pi/4,                  bin_input=True,                 soft_decision=False,                 bin_output=True)msg = np.array([0, 0, 0, 1, 1, 0, 1, 1]) # input messagemodulated = modem.modulate(msg) # modulationdemodulated = modem.demodulate(modulated) # demodulationprint("Modulated message:\n"+str(modulated))print("Demodulated message:\n"+str(demodulated)) >>>  Modulated message:   [0.70710678+0.70710678j  0.70710678-0.70710678j    -0.70710678+0.70710678j  -0.70710678-0.70710678j]>>> Demodulated message:   [0. 0. 0. 1. 1. 0. 1. 1.]

Или тоже QPSK, но уже с недвоичным входом / выходом:

import numpy as npfrom ModulationPy import PSKModemmodem = PSKModem(4, np.pi/4,                  bin_input=False,                 soft_decision=False,                 bin_output=False)msg = np.array([0, 1, 2, 3]) # input messagemodulated = modem.modulate(msg) # modulationdemodulated = modem.demodulate(modulated) # demodulationprint("Modulated message:\n"+str(modulated))print("Demodulated message:\n"+str(demodulated))>>> Modulated message:[ 0.70710678+0.70710678j -0.70710678+0.70710678j  0.70710678-0.70710678j -0.70710678-0.70710678j] >>> Demodulated message:[0, 1, 2, 3]

Пример для 16-QAM (десятичный вход / выход):

import numpy as npfrom ModulationPy import QAMModemmodem = PSKModem(16,                  bin_input=False,                 soft_decision=False,                 bin_output=False)msg = np.array([i for i in range(16)]) # input messagemodulated = modem.modulate(msg) # modulationdemodulated = modem.demodulate(modulated) # demodulationprint("Demodulated message:\n"+str(demodulated))>>> Demodulated message:[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]

В общем и целом, я думаю, понятно. Доступные опции старался делать по примеру доступных в рамках матлабовского Communication Toolbox. Подробное описание приведено в README.md проекта.

BER performance

Продемонстрируем адекватно ли модемы работают в случае присутствия шума (возьмем классический АБГШ, он же AWGN), используя простейшую модель приема-передачи:

Сверяться будем с теоретическими кривыми [5] (если кому интересно, все формулы описаны тоже в README.md).

Исходные коды для моделирования представлены по ссылкам: M-PSK, M-QAM.

Результаты:

Да, с огрехами на малых значениях битовых ошибок (из-за относительно небольшого количества усреднений, я полагаю), однако.... it works!

Производительность

А теперь я хотел бы вернуться к вопросу "соревнования" с упомянутой выше CommPy.

Что нас отличает:

• организация кода, стилистические различия (побочное, но не пренебрежимое);

• я использовал более быстрый алгоритм демодуляции [6] (подробно описан в матлабовской документации [7], ну, и я добавил все в тот же README.md).

И вот, что получилось "забенчмаркать":

• Скрипт: CommPy_vs_ModulationPy.ipynb

• Платформа: https://jupyter.org/try (Classic Notebook)

• Длина фрейма: 10 000 (элементов)

Результаты:

Разработчикам CommPy результаты, вроде, понравились (см. данный issue) - поэтому возможно в обозримом будущем что-то из моего ModulationPy будет перекочевывать в CommPy (я не против, главное, чтобы пользу приносило). Но это, как говорится, поживем - увидим.

И, да, пусть результаты производительности и не дотянули до MatLab (по крайней мере исходя из данного примера: см. вкладку "Examples"), я все равно считаю достигнутое неплохим стартом!

Послесловие

Наверное, проекту не хватает еще некоторых видов модуляции (тех же 32-QAM и 128-QAM или же используемой в DVB-S2/S2X APSK), однако, честно скажу, что не могу обещать их скорого добавления.

Проект всегда был для меня в большей мере площадкой для изучения языка Python и библиотеки NumPy на практике (и сопутствующих инструментов: юнит-тесты (не успел правда в данном случае перейти на pytest - каюсь), CI (использую Travis), подготовка модуля для PyPi и т.д.), однако, теперь, слава богу, всему этому есть приложение и в рамках рабочих задач!

Однако, все же буду рад вашим issue и pull request'ам! И если возьметесь интегрировать наработки в CommPy, тоже будет очень круто!

В общем, не серчайте, если вдруг не отвечу достаточно быстро, и да пребудет с вами сила науки!

Литература и ссылки

1. Haykin S. Communication systems. John Wiley & Sons, 2008. p. 93

2. Goldsmith A. Wireless communications. Cambridge university press, 2005. p. 88-92

3. MathWorks: comm.PSKModulator (https://www.mathworks.com/help/comm/ref/comm.pskmodulator-system-object.html?s_tid=doc_ta)

4. Octave: qammod (https://octave.sourceforge.io/communications/function/qammod.html)

5. Link Budget Analysis: Digital Modulation, Part 3 (www.AtlantaRF.com)

6. Viterbi, A. J. (1998). An intuitive justification and a simplified implementation of the MAP decoder for convolutional codes. IEEE Journal on Selected Areas in Communications, 16(2), 260-264.

7. MathWorks: Approximate LLR Algorithm (https://www.mathworks.com/help/comm/ug/digital-modulation.html#brc6ymu)