Русский
Русский
English
Статистика
Реклама

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

Автоматизация и промышленная электроника когда одним Arduino сыт не будешь

19.05.2021 12:19:31 | Автор: admin
Если играться с контроллерами, то почему с маленькими?

Очень часто, когда речь заходит об автоматизации чего-либо, то в разговоре всплывает Arduino, его производные или же Raspberry PI и прочие одноплатники. Но есть отличие от домашних поделок, где можно пользоваться чем угодно ради экономии и потому, что это простое и доступное решение. В сфере автоматизации/модернизации объектов, связанных с промышленностью, речь идёт исключительно о специализированных промышленных контроллерах и системах визуализации, диспетчеризации/удалённого управления и все это исключительно с сертификатами соответствия и лицензиями.
Решений такого класса море и порой сложно в них разобраться. Разумеется, все возможные варианты разобрать невозможно, но мы с коллегами уже несколько лет работаем в этой сфере и потому какое-то количество опыта набралось. Мы поделимся своим и если вам есть, что сказать просим писать комментарии.

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

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

Самая маленькая команда должна включать в себя:

  • проектанта
  • опытного инженера в области промышленной электроники
  • электромонтеров/электромехаников
  • программиста

Если же компания постоянно развивается и работает в этой области, то она неуклонно растёт, так как несколько объектов делается одновременно и все они находятся на разной стадии.

Стадия подготовки к тендеру


О начале тендера, на выполнение чего-то, мы не всегда узнаем заранее. Иногда информация запаздывает, тендер ожидается на днях, или долго шли подготовительные работы, чтобы принять участие. Но это не меняет сути: нужно посчитать стоимость автоматизации объекта, к примеру. Сначала изучается планируемый объем и срок работ (если, конечно, сроки внятно прописаны) и принимается решение о выполнимости своими силами, либо с привлечением подрядчиков.
Часто в нагрузку к системам автоматики докидывают смежные (в понимании заказчика/инвестора) вещи. Нам постоянно навязывают обычную электрику, вплоть до внешнего освещения, локальные сети и телефонные линии в промышленных и административных помещениях, электрические сети низкого напряжения и понижающие трансформаторные подстанции. И вот хотелось бы иметь только автоматику в чистом виде, да нет её. Так как для инвесторов, что электроника, что электрика одного поля ягоды. Попутно нужно понимать, что компания, занимающаяся автоматизацией, заканчивает объект последней. Именно она будет связывать воедино разрозненные структуры пожарные системы, контроли доступа, видео наблюдение, связь с серверной и оптическими линиями, диспетчеризацию и удалённый доступ. Если же технологический процесс объекта содержит в себе независимые самоуправляемые системы или машины (автоматические линии, вентиляционные централи, сигнализацию, запасной генератор, солнечные панели и преобразователи к ним и т.д.), то задача усложняется ещё и необходимостью соединить между собой оборудование разных производителей, порой с совсем разными протоколами промышленной связи.
Вооружившись знаниями за все предыдущие годы и выполненные объекты, можно попробовать посчитать общую себестоимость работ (включая материалы, покупку оборудования, количество человеко-часов, стоимость всех выездов на объект для проведения монтажа, пуско-наладки и последующих сервисных выездов в течении гарантийных 5-7 лет).
Q: На основании чего производится расчёт?
Пожеланий инвестора и приблизительных догадок об общем бюджете объекта. Всегда можно воспользоваться открытой информацией. Оттуда можно выудить сколько процентов от общей суммы пойдёт на данную сферу. Грубо говоря, если объект будет стоить 10 млн евро, то системы автоматизации плюс электрика это около 1,5-2 млн.
Хорошо, когда есть возможность пообщаться с генеральным подрядчиком, который координирует весь ход работ, от выравнивания территории под фундаменты и дороги, до капитального строительства и озеленения территории. Всегда есть кто-то, кто пытается, не выходя за рамки первичного бюджета, построить объект и сдать его срок. Не завидуем им совсем, но пытаемся сотрудничать на всех этапах реализации. Это действительно упрощает жизнь в дальнейшем.
Итак, получаем информацию от генерального подрядчика, что общая квота от этого бюджета на автоматику и электрику не сможет превысить тех же 2 млн. Теперь начинается самое интересное в тендере участвует минимум 3 компании на каждую сферу работ. А значит мы встретимся ещё, как минимум, с двумя предложениями от других компаний, что претендуют выполнить этот же объем работ, что и мы. Нужно предложить на тендере сумму меньше них (а мы, разумеется не знаем сколько они заявят), ну и не остаться без прибыли, тем более не уйти в минус.
Q: Что мы используем и почему?
Зависит от возможностей оборудования и пожеланий заказчика. Из контроллеров это Siemens s1200-s1500 и его аналоги. Панели к ним Weintek и Astraada. Реле, пускатели, клеммные колодки Eaton и т.д. Частотных преобразователи Danfoss, Mitsubishi, LS IS. Коммуникация Modbus RTU, Modbus TCP, Ethernet, Profibus и оптика. Подробнее о подборе оборудования будет в следующей статье относительно проектирования.

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

Сначала изучается ТЗ.Как известно без внятного ТЗ получается непонятно что. Инвестор не будет пояснять расплывчатые формулировки и нужно готовиться ко всему. Подводные камни ждут в каждом предложении и нужно уметь их находить. К примеру, если пишется, что необходимо выполнить систему автоматизации и удалённого управления с возможностью работы при кратковременных пропаданиях питания от сети, то это означает, что всего лишь необходимо доложить по одному блоку бесперебойного питания либо аккумуляторы в каждых шкаф управления техпроцессом и обычный UPS, через который запитана серверная стойка или компьютер диспетчера с системой управления. Задача поддержать работу всех контролёров, во всех шкафах управления и всех роутеров и модемов, чтобы собрать все статусы ошибок различных устройств, автоматических линий и машин, записать все это в логи и передать на сервер/компьютер с системой удалённого управления, который через SMS-шлюз сообщит наружу подготовленному списку номеров, что пропало напряжение сети и все оборудование остановилось.
Как можно понять из требований, это достаточно недорогое удовольствие обойдётся в несколько сотен евро, и его легко учесть на стадии проектирования. Однако формулировка может быть несколько иной выполнить систему автоматизации с возможностью продолжения работы объекта при пропадании питания от сети, до момента её появления. И вот это уже совсем другое. Здесь хотят иметь либо вторую линию резервного питания со своим понижающим трансформатором и шкаф автоматического переключения нагрузки с одного источника питания на другой с анализаторами параметров сети. Либо, что встречается гораздо чаще заказчик хочет иметь резервный источник питания в виде дизельного генератора, который сможет поддержать работу объекта в течении нескольких часов. Подобные генераторы стоят десятки тысяч евро. И не забываем про шкаф переключения нагрузки, благо в большинстве случаев современные дизельные генераторы имеют свои неплохие контроллеры с внешним LCD экраном и возможностью отслеживания параметров сети сразу по нескольким параметрам и управлением силовыми контакторами, что переключают нагрузку от одного источника (сети) на другой (генератор) и наоборот при восстановлении нормальной работы первичного источника.
К примеру, генераторы фирмы FOGO, Elem Power или их аналоги, с контроллером, к примеру, datacom-d500-lite или inteli lite amf 25. Тут шкаф переключения нагрузки (AVR) будет состоять всего из 2х силовых контакторов с дополнительными контактами, чтобы снимать состояние, какой из них включён и сделать их одновременно включение взаимоисключающим.

Фото панели управления генератором

После прочтения ТЗ все обычные работы:прокладка кабельных трасс, монтаж кабельных лотков, монтаж наружного освещения, заземления, громоотвода и молниезащиты (если есть в требованиях), изготовление шкафов управления (их мы делаем сами), всё калькулируется по нормам. В нашем случае эти нормы европейские. Потом прикидывается общая стоимость всех возможных шкафов управления с их наполнением в виде контроллеров и их обвязки, частотных преобразователей или устройств плавного пуска, серверных стоек, компьютеров и телевизоров (это относительная новинка, так как стали отказываться от синоптических таблиц и предпочитают иметь выведенную систему визуализации техпроцесса на 4-6 телевизоров по 50-55", склеенных в мульти-экран), считается за сколько рабоче-часов будет изготовлено, смонтировано. Сюда же добавляется прибыль. Есть несколько быстрых шаблонов по подсчёту стоимости к примеру, в конечную стоимость одного шкафа управления, включена оплата работы проектанта, стоимость сборки шкафа и стоимость программирования контроллера и панели HMI. Чтобы быстро считать многие вещи упрощаются. Итоговая стоимость шкафа управления будет = себестоимость всех компонентов * 1,5(2). Разумеется, сам по себе шкаф ничего не делает. Его ещё нужно привезти на объект, смонтировать на своём месте, проложить к нему кабельные линии и трассы, уложить питающие и управляющие кабеля, подключить это все и запустить. Приблизительно можно подсчитать общую длину будущих кабельных линий, учесть проколы под дорогами (это дорого), посмотреть требования к силовым линиям (медные или алюминиевые), подумать на счёт трансформаторной подстанции и генератора. Многие вещи сразу умножаются на 2, чтобы избежать нехватки средств в будущем. К примеру, тот же генератор весит 3-5 тонн. Значит нужно организовать выгрузку с помощью крана и прочный фундамент, соответственно ещё и заземление. Внимательно изучаются расположения разных технологических зданий внутри объекта. По современным нормам уже никто не хочет связывать отдалённые здания витой парой хоть 5й, хоть 6й категории. Только оптика, только хардкор. Считаем и стоимость укладки, сварки оптики и оборудования для соединения между собой разных зданий. Внутри зданий и на малых расстояниях, для экономии, используем utp категории 5е/6e.

Готовый шкаф управления. Запуск ещё идёт, но большая часть техпроцесса уже налажена

После того, как мы все посчитали, оказывается что мы как обычно не вписываемся в ожидаемые рамки бюджета. В ход идут оптовые закупки оборудования напрямую у производителя и соответственно скидки за многолетнюю лояльность. Иногда скидки могут достигать десятков процентов. Это ничего не значит, если вы покупаете на объект 1 контроллер за тысячу евро и один частотный преобразователь за пару сотен. Но когда суммарно на объекте несколько десятков частотных преобразователей, несколько контроллеров и несколько шкафов, то экономия получается порядочной. Это наша подушка безопасности. Если что-то упущено на стадии подготовки к тендеру, то мы используем сэкономленные средства на это. Если же все хорошо, то это плюс к нашей прибыли. В конечном счёте, никто из нас не альтруист и все мы здесь за этим. Хотя есть ещё и моральное удовольствие от хорошо выполненной работы. Но это уже невозможно посчитать)

Фото запущенного объекта на память)

Вот, пожалуй, и все, что можно в нескольких словах сказать о подготовке к тендеру и подводных камнях на этом этапе. Большие объекты, вроде Amazon, Volkswagen, АЭС (да, в нашей практике и такое было) строятся не быстро, но и бюджеты отнюдь не маленькие.
В одной статье описать все стадии автоматизации невозможно, далее планирую затронуть стадии проектирования, программирования и запуска оборудования. Пожалуй, больше всего интересного происходит во время запуска. Именно там встречается масса проблем и находятся решения.


Подробнее..

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

14.05.2021 18:17:37 | Автор: admin

Эта статья является продолжением двух предыдущих: Ударим биспектром по бездорожью, или как найти золото в Сибири, в которой мы рассмотрели геологическую модель месторождения золота на территории Новосибирской области и Ищем рудное золото на острове Сумбава, Индонезия, в которых мы построили геологически обусловленную модель машинного обучения для поиска золота или других рудных минералов по всему Тихоокеанскому рудному поясу, используя для геологического моделирования открытые данные на платформе Google Earth Engine (GEE).


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



Сопоставление моделей на территории Индонезии и Сибири


У нас есть уже обученный на детальных открытых данных классификатор для Индонезии, который мы хотели бы использовать и в Сибири. Сравним детальные модели плотности без учета влияния рельефа, построив для Западной Сибири модель с теми же параметрами, что и для Западной Сумбавы и выбрав равные участки. На левой модели точками на поверхности показаны скважины с высокими (Au) и низкими пробами золота (No Au), а на правой контурами отображены известные региональные месторождения золота (Au) и молибдена (Mo):



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


Классификация золоторудных территорий


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



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


Аналогично применим построенный классификатор к участку Западной Сибири (очевидно, здесь мы не можем оценить точность классификации автоматически):



Визуально результат соответствует ожидаемому и классификатор выделяет все те признаки, которые мы перечислили выше. Проверим на 3D модели:



Полученный результат геологически обоснован, поскольку выделенные области находятся в зонах накопления на границе выходов твердых пород, как мы это подробно разбирали ранее. При этом в области молибденового месторождения тоже выделены потенциальные золотые депозиты, что не удивительно, поскольку на региональной карте месторождения указаны по преобладающему минералу. Кроме того, на онлайн карте N-45-XIII. Схема прогноза полезных ископаемых, м-б 1:500 000 эти два месторождения объединены как золото-молибденовое месторождение.


Заключение


Начиная с пары статей по использованию платформы Google Earth Engine (GEE) и получению серий спутниковых данных для любой территории планеты, мы продолжили построением геологических моделей в разных регионах с последующей оценкой их геологической достоверности и автоматизированной классификацией, и, наконец, совместили полученные результаты, создав автоматизированное решение построения детальных карт золоторудности разных территорий.


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


Ссылки


Подробнее..

CatBoost и ML-конкурсы

15.05.2021 14:20:40 | Автор: admin

Анализ данных и базоваямодель

Вступление

Эта статья основана на данных конкурса, который компания Driven Data опубликовала для решения проблем с источниками воды в Танзании.

Информация для конкурса была получена Министерством водных ресурсов Танзании с использованием платформы с открытым исходным кодом под названием Taarifa. Танзаниясамая большая страна в Восточной Африке с населением около 60 миллионов человек. Половина населения не имеет доступа к чистой воде, а 2/3 населения страдает от плохой санитарии. В бедных домах семьям часто приходится тратить несколько часов пешком, чтобы набрать воду из водяных насосов.

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

Данные

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

Пункты водоснабжения разделены на исправные, нефункциональные и исправные, но нуждающиеся в ремонте. Цель конкурсапостроить модель, прогнозирующую функциональность точек водоснабжения.
Данные содержат 59400 строк и 40 столбцов. Целевая метка содержится в отдельном файле.
Показатель, используемый для этого соревнования,это classification rate, который вычисляет процент строк, в которых прогнозируемый класс совпадает с фактическим классом в тестовом наборе. Максимальное значение равно 1, а минимальное0. Цель состоит в том, чтобы максимизировать classification rate.

Анализ данных

Описания полей в таблице с данными:

  • amount_tshобщий статический напор (количество воды, доступное для точки водоснабжения)

  • date_recordedдата сбора данных

  • funderкто спонсировал постройку колодца

  • gps_heightвысота на которой находится колодец

  • installerорганизация построившая колодец

  • longitudeGPS координаты (долгота)

  • latitudeGPS координаты (широта)

  • wpt_nameназвание, если оно есть

  • num_privateнет информации

  • basinгеографический водный бассейн

  • subvillageгеографическая локация

  • regionгеографическая локация

  • region_codeгеографическая локация (код)

  • district_code географическая локация (код)

  • lgaгеографическая локация

  • ward географическая локация

  • populationколичество населения около колодца

  • public_meetingда/нет

  • recorded_by кто собрал данные

  • scheme_managementкто управляет колодцем

  • scheme_nameкто управляет колодцем

  • permitтребуется ли разрешение на доступ

  • construction_yearгод постройки

  • extraction_typeтип колодца

  • extraction_type_groupгруппа типа колодца

  • extraction_type_classкласс типа колодца

  • managementкак управляется колодец

  • management_groupгруппа управления колодца

  • paymentстоимость воды

  • payment_typeтип стоимости воды

  • water_qualityкачество воды

  • quality_groupгруппа качества воды

  • quantityколичество воды

  • quantity_groupгруппа количества воды

  • sourceисточник воды

  • source_typeтип источника воды

  • source_classкласс источника воды

  • waterpoint_typeтип колодца

  • waterpoint_type_groupгруппа типа колодца

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

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

    • уменьшать количество примеров у преобладающих классов(under-sampling)

    • дублировать строки с метками, которых мало(over-sampling)

    • генерировать синтетические выборкиэто случайная выборка атрибутов из экземпляров в классе меньшинства (SMOTE)

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

    Больше про тактики работы с несбалансированными датасетами можно почитать здесь.

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

    Некоторые данные содержат пустые значения.

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

    Следующая тепловая карта представляет наличие/отсутствие связи между переменными. Стоит обратить внимание на соотношение между permit, installer и funder.

    Давайте посмотрим на общую картину взаимоотношений на дендрограмме.

    В характеристиках водяных насосов есть та, которая показывает количество воды. Мы можем проверить, как количество воды связано с состоянием насосов (quantity_group).

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

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

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

    Большинство колодцев с неизвестным качеством из quality_group не работают.

    Есть еще одна интересная характеристика водных точеких тип (waterpoint_type_group).

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

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

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

    Danidaсовместная организация Танзании и Дании по финансированию скважин, и хотя у них много работающих водозаборов, процент неисправных очень высок. Похожая ситуация с RWSSP (программа сельского водоснабжения и канализации), Dhv и некоторыми другими. Следует отметить, что большинство скважин, профинансированных Германией и частными лицами, находятся в рабочем состоянии. Напротив, большое количество скважин, которые финансируются государством, не функционируют. Большинство пунктов водоснабжения, установленных центральным правительством и районным советом, также не работают.

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

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

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

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

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

    Часть данных была заполнена значением 0 вместо реальных данных. Мы также можем видеть, что amount_tsh выше у исправных колодцев (label = 0). Также следует обратить внимание на выбросы в функции amount_tsh. В качестве особенности можно отметить перепад высот и тот факт, что значительная часть населения проживает на высоте 500 метров над средним уровнем моря.

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

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

  • Фича installer содержит множество повторов с разными регистрами, орфографическими ошибками и сокращениями. Давайте сначала переведем все в нижний регистр. Затем по простым правилам уменьшаем количество ошибок и делаем группировку.

  • После очистки мы заменяем все элементы, которые встречаются менее 71 раз (0,95 квантиля), на other.

  • Повторяем по аналогии с фичей funder. Порог отсечки98.

  • Данные содержат фичи с очень похожими категориями. Выберем только по одной из них. Поскольку в датасете не так много данных, мы оставляем функцию с наименьшим набором категорий. Удаляем следующие фичи: scheme_management, quantity_group, water_quality, payment_type, extraction_type, waterpoint_type_group, region_code.

  • Заменим latitude и longitude значения у выбросов медианным значением для соответсвующего региона из region_code.

  • Аналогично повторям для пропущенных значений для subvillage и scheme_name.

  • Пропущенные значения в public_meeting и permit заменяем медианным значением.

  • Для subvillage, public_meeting, scheme_name, permit, создаем дополнительные категориальные бинарные фичи, которы будут отмечать данные с пропущенными значениями. Так как мы заменяем их на медианные, то для модели оставим информацию, что пропуски были.

  • Фичи scheme_management, quantity_group, water_quality, region_code, payment_type, extraction_type, waterpoint_type_group, date_recorded, и recorded_by можно удалить, так как там повторяются данные из других фичей, для модели они будут бесполезны.

Модель

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

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

def fit_model(train_pool, test_pool, **kwargs):    model = CatBoostClassifier(        max_ctr_complexity=5,        task_type='CPU',        iterations=10000,        eval_metric='AUC',        od_type='Iter',        od_wait=500,        **kwargs    )return model.fit(        train_pool,        eval_set=test_pool,        verbose=1000,        plot=False,        use_best_model=True)

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

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

def classification_rate(y, y_pred):    return np.sum(y==y_pred)/len(y)

Поскольку данных мало, разбивать полную выборку на обучающую и проверочную частине лучший вариант. В этом случае лучше использовать методику OOF (Out-of-Fold). Мы не будем использовать сторонние библиотеки; давайте попробуем написать простую функцию. Обратите внимание, что разбиение набора данных на фолды необходимо стратифицировать.

def get_oof(n_folds, x_train, y, x_test, cat_features, seeds):    ntrain = x_train.shape[0]    ntest = x_test.shape[0]              oof_train = np.zeros((len(seeds), ntrain, 3))    oof_test = np.zeros((ntest, 3))    oof_test_skf = np.empty((len(seeds), n_folds, ntest, 3))    test_pool = Pool(data=x_test, cat_features=cat_features)     models = {}    for iseed, seed in enumerate(seeds):        kf = StratifiedKFold(            n_splits=n_folds,            shuffle=True,            random_state=seed)                  for i, (train_index, test_index) in enumerate(kf.split(x_train, y)):            print(f'\nSeed {seed}, Fold {i}')            x_tr = x_train.iloc[train_index, :]            y_tr = y[train_index]            x_te = x_train.iloc[test_index, :]            y_te = y[test_index]            train_pool = Pool(data=x_tr, label=y_tr, cat_features=cat_features)            valid_pool = Pool(data=x_te, label=y_te, cat_features=cat_features)model = fit_model(                train_pool, valid_pool,                loss_function='MultiClass',                random_seed=seed            )            oof_train[iseed, test_index, :] = model.predict_proba(x_te)            oof_test_skf[iseed, i, :, :] = model.predict_proba(x_test)            models[(seed, i)] = modeloof_test[:, :] = oof_test_skf.mean(axis=1).mean(axis=0)    oof_train = oof_train.mean(axis=0)    return oof_train, oof_test, models

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

Кривая обучения одного изфолдовКривая обучения одного изфолдов

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

Посмотрев на важность фичей в модели, мы можем убедиться в отсутствии явной утечки информации (лика).

После усреднения прогнозов:

balanced accuracy: 0.6703822994494413classification rate: 0.8198316498316498

Попробуем загрузить результаты на сайт с конкурсом и посмотрим на результат.

Учитывая, что на момент написания этой статьи результат топ-5 был только на 0,005 лучше, мы можем сказать, что базовая модель получилась хорошей.

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

balanced accuracy: 0.6549535670689709classification rate: 0.8108249158249158

Результат заметно хуже.

Итоги

В статье:

  • познакомились с данными и поискали идеи для модели;

  • очистили и подготовили данные для модели;

  • приняли решение использовать CatBoost, так как основная масса фичей категориальные;

  • написали функцию для OOF-предсказания;

  • получили отличный результат для базовой модели.

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

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

Код из статьи можно посмотреть здесь.

Подробнее..

Легенды и мифы геофизики

23.05.2021 12:20:57 | Автор: admin

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



Видите взаимосвязь ортофотоснимка и рельефа? Если да, то вы или геолог или можете им стать: корреляция компонентов (разложения в пространственный спектр) составляет 41% для длины волны 20 м, 58% для 50 м и 99% для 300 м (Jupyter Python ноутбук с вычислениями доступен по ссылкам ниже). Большинство геофизиков поклянутся, что это у вас спектры порченые (записано с натуры), игнорируя и геофизику и прилагаемые вычисления и ссылки на публикации.


Миф первый, рельеф и сила тяжести


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


Важнейшим доказательством изостазии является отсутствие связи между рельефом и силой тяжести.

Увы, именно в отсутствие такой связи верят многие из встреченных мной геофизиков (но не геологи), хотя это абсолютно неверно! Что интересно, в англоязычной версии этой же викистраницы написано прямо противоположное, смотрите пример с айсбергом в разделе Deposition and erosion:


An analogy may be made with an iceberg, which always floats with a certain proportion of its mass below the surface of the water. If snow falls to the top of the iceberg, the iceberg will sink lower in the water. If a layer of ice melts off the top of the iceberg, the remaining iceberg will rise. Similarly, Earth's lithosphere "floats" in the asthenosphere.

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


Посмотрим, что говорит оглавление русскоязычного учебника по геофизике (Викулин, 2004):


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

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


Более того, корреляция компонентов силы тяжести (или вертикальной компоненты гравитационного поля) и рельефа близка к 100%, как показано в работе "Radially symmetric coherence between satellite gravity and multibeam bathymetry grids" (Marks and Smith, 2012). На графике ниже график радиальной когеренции (справа) из указанной статьи дополнен мною двумерной коррелограммой (слева) для глобальной модели гравитационных аномалий Sandwell and Smith Gravity Anomaly и глобальной модели рельефа GEBCO 2019:



Если на графике радиальной когеренции (корреляции) не понятно, почему для первых километров масштаба связь компонентов мала, то на двумерной коррелограмме все становится ясно измеренная сила тяжести относится к поверхности океана, а топография, очевидно, ко дну. Таким образом, здесь мы видим, что до масштабов ~40 км корреляция компонент близка к 100%. Более того, для региона Индонезия высокая корреляция (>75%) наблюдается до масштабов в сотни километров и средняя (>50%) для масштабов в тысячи километров, как показано в статьях по ссылкам ниже.


Миф второй, рельеф и аномалия Буге


Что интересно, миф второй полностью противоречит первому, что не мешает многим геофизикам верить в оба разом. Встречайте Аномалию Буге (Бугера):


Аномалия Бугера расхождение между теоретически вычисленным и реальным значением гравитационного поля Земли в определённой точке референц-эллипсоида.

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


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

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


И следующий шаг, или редукция:


Уточнённая, или полная редукция Бугера позволяет полностью учесть влияние рельефа местности.

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


Посмотрим практический пример, как компоненты силы тяжести глобальной модели WGM2012 после редукции Буге (слева) и еще одной редукции (справа) коррелируют с компонентами рельефа глобальной модели GEBCO 2019:



Как видим, про "полностью учесть рельеф местности" речь тут явно не идет то есть результат редукции Буге противоречит его определению, а значит, его использование неправомерно. Обратим внимание на правый график тут мы ожидаем максимум корреляции 100%, а наблюдаем всего около 75%, что показывает неточность самой модели WGM2012 (по сравнению с использованной выше Sandwell and Smith Gravity Anomaly, авторы которой вовсе не вычисляют редукцию Буге по уже понятным причинам), таким образом, здесь реальная корреляция редукции Буге с рельефом получится около -100%.


Если обратиться к учебнику "Гравиразведка" (Утёмов, 2009), там в соответствующем разделе есть разъяснение полученного результата:


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

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


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


Миф третий, сферический конь эллипсоид в вакууме


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


В геофизике модель геопотенциала представляет собой теоретический анализ измерения и расчета эффекты Земли гравитационного поля

Да-а, расчета эффекты Земли это, конечно, не совсем русский язык, или даже совсем не русский. Равно как и теоретический анализ измерения тоже тот еще шедевр. Тем не менее, можно понять, что модель геопотенциала это модель (глобального) гравитационного поля Земли. Это определение уже ошибочное на самом деле, рассматривается вертикальная компонента гравитационного поля, то есть сила тяжести (гравитационное поле векторное, а сила тяжести скалярная, и связь между ними весьма сложная). Но, главное, геофизикам такая модель вообще не нужна! А кому нужна, описано дальше в этой же викистатье:


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

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


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


Миф четвертый, о пространственных спектрах


У меня в коллекции есть такие перлы, как спектры у вас порченые (освятить, что ли, хотят?) и другие. Попросту говоря, геофизики не имеют никакого понятия, что это такое.


Давайте посмотрим, как же практически работать с пространственными спектрами. Как всегда в математике, все эти преобразования по сути чрезвычайно просты. Например, результат применения гауссова фильтра масштаба N метров есть низкочастотная фильтрация (low-pass filter), разница между исходным изображением и фильтрованным есть высокочастотная фильтрация (high-pass filter) и разница между двумя фильтрами масштаба N и M метров есть полосовая фильтрация (band-pass filter). Для выделения компонента пространственного спектра шириной l метров для длины волны L метров достаточно выполнить полосовую фильтрацию с фильтрами масштаба L-l/2, L+l/2. Спектр мощности может быть представлен как стандартное отклонение полученных компонент, а через логарифм спектра мощности посчитана и фрактальная размерность. Для примера посмотрите график радиальной когерентности из статьи НАСА и двумерную спектрограмму, полученную указанным методом, в обсуждении первого мифа выше.


Низкочастотная фильтрация это улучшенная замена Буге преобразования:



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



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


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


Заключение


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


Бонус


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


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

Так что учтите если вы ссылаетесь на результаты из публикаций, скажем, НАСА (аналогично корреляции между пространственными спектрами рельефа и гравиметрии посчитав таковую между рельефом и космическими снимками) то будете обвинены в ереси ("математических манипуляциях") и все ваши "порченые спектры" придется немедленно освятить (боюсь, путем окунания компьютера в чан со святой водой и никак иначе). Я предупредил! :)


Ссылки


Marks, K.M., Smith, W.H.F., 2012. Radially symmetric coherence between satellite gravity and multibeam bathymetry grids. Mar Geophys Res 33, 223227. https://doi.org/10.1007/s11001-012-9157-1


Викулин А.В., 2004. ВВЕДЕНИЕ В ФИЗИКУ ЗЕМЛИ. УЧЕБНОЕ ПОСОБИЕ ДЛЯ ГЕОФИЗИЧЕСКИХ СПЕЦИАЛЬНОСТЕЙ ВУЗОВ


Утёмов Э.В., 2009. Гравиразведка


Spectral Coherence between Gravity and Bathymetry Grids


WGM2012 spatial components of free-air gravity and topography correlation


Bouguer and Free-Air Gravity Anomalies in terms of spatial spectrum


There is a high correlation between DEM and ortho photos


3D Density Inversion by Circular Hough Transform [Focal Average] and Fractality Index


[Gaussian Filtering on Spheroid [Sandwell & Smith]](

Подробнее..

Перевод 5 разных библиотек Python, которые сэкономят ваше время

12.06.2021 18:20:44 | Автор: admin

В этой подборке, переводом которой мы решили поделиться к старту курса о машинном и глубоком обучении, по мнению автора, каждая библиотека заслуживает отдельной статьи. Всё начинается с самого начала: предлагается библиотека, которая сокращает шаблонный код импортирования; заканчивается статья пакетом удобной визуализации данных для исследовательского анализа. Автор также касается работы с картами Google, ускорения и упрощения работы с моделями ML и библиотеки, которая может повысить качество вашего проекта в области обработки естественного языка. Посвящённый подборке блокнот Jupyter вы найдёте в конце.


PyForest

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

Вот почему PyForest это одна из самых удобных библиотек, которые я знаю. С её помощью в ваш блокнот Jupyter можно импортировать более 40 популярнейших библиотек (Pandas, Matplotlib, Seaborn, Tensorflow, Sklearn, NLTK, XGBoost, Plotly, Keras, Numpy и другие) при помощи всего одной строки кода.

Выполните pip install pyforest. Для импорта библиотек в ваш блокнот введите команду from pyforest import *, и можно начинать. Чтобы узнать, какие библиотеки импортированы, выполните lazy_imports().

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

Emot

Эта библиотека может повысить качество вашего проекта по обработке естественного языка. Она преобразует эмотиконы в их описание. Представьте, например, что кто-то оставил в Твиттере сообщение I [здесь в оригинале эмодзи "красное сердце", новый редактор Хабра вырезает его] Python. Человек не написал слово люблю, вместо него вставив эмодзи. Если твит задействован в проекте, придётся удалить эмодзи, а значит, потерять часть информации.

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

Чтобы установить Emot, выполните команду pip install emot, а затем командой import emot импортируйте её в свой блокнот. Нужно решить, с чем вы хотите работать, то есть с эмотиконами или с эмодзи. В случае эмодзи код будет таким: emot.emoji(your_text). Посмотрим на emot в деле.

Выше видно предложение I [эмодзи "красное сердце"] Python, обёрнутое в метод Emot, чтобы разобраться со значениями. Код выводит словарь со значением, описанием и расположением символов. Как всегда, из словаря можно получить слайс и сосредоточиться на необходимой информации, например, если я напишу ans['mean'], вернётся только описание эмодзи.

Geemap

Говоря коротко, с её помощью можно интерактивно отображать данные Google Earth Engine. Наверное, вы знакомы с Google Earth Engine и всей его мощью, так почему не задействовать его в вашем проекте? За следующие несколько недель я хочу создать проект, раскрывающий всю функциональность пакета geemap, а ниже расскажу, как можно начать с ним работать.

Установите geemap командой pip install geemap из терминала, затем импортируйте в блокнот командой import geemap. Для демонстрации я создам интерактивную карту на основе folium:

import geemap.eefolium as geemapMap = geemap.Map(center=[40,-100], zoom=4)Map

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

Dabl

Позвольте мне рассказать об основах. Dabl создан, чтобы упростить работу с моделями ML для новичков. Чтобы установить её, выполните pip install dabl, импортируйте пакет командой import dabl и можно начинать. Выполните также строчку dabl.clean(data), чтобы получить информацию о признаках, например о том, есть ли какие-то бесполезные признаки. Она также показывает непрерывные, категориальные признаки и признаки с высокой кардинальностью.

Чтобы визуализировать конкретный признак, можно выполнить dabl.plot(data).

Наконец, одной строчкой кода вы можете создать несколько моделей вот так: dabl.AnyClassifier, или так: dabl.Simplefier(), как это делается в scikit-learn. Но на этом шаге придётся предпринять некоторые обычные шаги, такие как создание тренировочного и тестового набора данных, вызов, обучение модели и вывод её прогноза.

# Setting X and y variablesX, y = load_digits(return_X_y=True)# Splitting the dataset into train and test setsX_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)# Calling the modelsc = dabl.SimpleClassifier().fit(X_train, y_train)# Evaluating accuracy scoreprint(Accuracy score, sc.score(X_test, y_test))

Как видите, Dabl итеративно проходит через множество моделей, включая Dummy Classifier (фиктивный классификатор), GaussianNB (гауссовский наивный Байес), деревья решений различной глубины и логистическую регрессию. В конце библиотека показывает лучшую модель. Все модели отрабатывают примерно за 10 секунд. Круто, правда? Я решил протестировать последнюю модель при помощи scikit-learn, чтобы больше доверять результату:

Я получил точность 0,968 с обычным подходом к прогнозированию и 0,971 с помощью Dabl. Для меня это достаточно близко! Обратите внимание, что я не импортировал модель логистической регрессии из scikit-learn, поскольку это уже сделано через PyForest. Должен признаться, что предпочитаю LazyPredict, но Dabl стоит попробовать.

SweetViz

Это low-code библиотека, которая генерирует прекрасные визуализации, чтобы вывести ваш исследовательский анализ данных на новый уровень при помощи всего двух строк кода. Вывод библиотеки интерактивный файл HTML. Давайте посмотрим на неё в общем и целом. Установить её можно так: pip install sweetviz, а импортировать в блокнот строкой import sweetviz as sv. И вот пример кода:

my_report = sv.analyze(dataframe)my_report.show_html()

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

Заключение

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

Этот материал не только даёт представление о полезных пакетах экосистемы Python, но и напоминает о широте и разнообразии проектов, в которых можно работать на этом языке. Python предельно лаконичен, он позволяет экономить время и в процессе написания кода, выражать идеи максимально быстро и эффективно, то есть беречь силы, чтобы придумывать новые подходы и решения задач, в том числе в области искусственного интеллекта, получить широкое и глубокое представление о котором вы можете на нашем курсе "Machine Learning и Deep Learning".

Узнайте, как прокачаться и в других специальностях или освоить их с нуля:

Другие профессии и курсы
Подробнее..

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

16.06.2021 22:19:34 | Автор: admin

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

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

В результате родился проект Объясняем код. Посмотреть, что это такое можно на code-explained.com. Код проекта выложен на Гитхаб.

Чем я вдохновлялся

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

Например, если вы редактируете программный код, то результат должен отображаться в реальном времени, без повторной компиляции. Иначе, вы не увидите промежуточные состояния, и у вас не родятся связанные с ними идеи. Я не буду подробно пересказывать эту мысль, тем более, что на Хабре уже подробно о ней рассказали. Почитайте, если еще не успели.

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

Как сегодня изучают алгоритмы

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

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

На ИТ-ресурсах создают анимации.

А кто-то даже пытается объяснить работу алгоритма через танец.

Почему эти подходы казались мне неэффективными

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

Как я перешел от технозависимости к человечности

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

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

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

Как это технически реализовано

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

Чтобы сэкономить на копировании данных, я применил Immutable.js. Immutable.js это библиотека неизменяемых коллекций. Модификации таких коллекций всегда возвращают новую коллекцию, которую легко сохранить. Так я избегаю глубокого копирования. Не буду вдаваться в подробности, на Хабре про immutable.js уже писали. Для примера кусочек кода с итерированием по хеш-таблице:

while (true) { // Итерируемся по списку    this.addBP('check-not-found'); // Метод сохраняем состояние    if (this.newList.get(this.newListIdx) === null) {        // this.newList -- это немутабельный список        break;    }    this.addBP('check-found'); // Выполнена очередная строчка, сохраняем состояние    if (EQ(this.newList.get(this.newListIdx), this.number)) {        this.addBP('found-key');        return true;    }    this.fmtCollisionCount += 1; // Для динамических комментариев иногда нужно сохранять статистикуу    this.newListIdx = (this.newListIdx + 1) % this.newList.size; // Переходим к следующему индекксу    this.addBP('next-idx');}

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

Много строк возникает из-за особенностей браузерных API и производительности в разных браузерах. Например, большой проблемой оказалось сделать так, чтобы браузеры не склеивали последовательные изменения друг с другом. Если добавить div с определённой начальной позицией, и потом сразу же поменять координаты на конечные, то браузер склеит эти два изменения в одно. Div сразу окажется в конечной позиции без анимации. Чтобы такое не происходило, приходится вставлять задержку в два фрейма анимации с помощью window.requestAnimationFrame().

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

Код проекта на гитхабе

Что дальше?

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

Подробнее..

Геопространственное моделирование с применением методов машинного обучения

18.06.2021 18:20:53 | Автор: admin


Всем привет! Меня зовут Константин Измайлов, я руководитель направления Data Science в Delivery Club. Мы работаем над многочисленными интересными и сложными задачами: от формирования классических аналитических отчетов до построения рекомендательных моделей в ленте приложения.

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

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

Статья написана по мотивам выступления с Евгением Макиным на конференции Highload++ Весна 2021. Для тех, кто любит видео, ищите его в конце статьи.

Бизнес-модель работы Delivery Club


Бизнес-модель Delivery Club состоит из двух частей:

  • ДДК (доставка Деливери Клаб): мы передаем заказ в ресторан и доставляем его клиенту, то есть ресторану остается только приготовить заказ к определенному времени, когда придет курьер.
  • МП (маркетплейс): мы передаем заказ в ресторан, а он своими силами доставляет заказ в пределах своей согласованной зоны.

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

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

Рисуем зону доставки ресторана


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



Как процесс выглядел раньше


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

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



Стоит упомянуть и про SLA (Service Level Agreement соглашение о максимальной длительности отрисовки зоны доставки для одного партнера): онбординг партнера или подготовка его зоны для внедрения в нашу платформу составляли порядка 40 минут для одного заведения. Представьте, что к вам подключилась городская сеть с сотней ресторанов, а если это ещё и жаркий сезон, например, после проведения рекламной акции Вот наглядное доказательство неэффективности ручной отрисовки:

$T = R * SLA = 100 * 40\ минут =\ \sim 67\ часов\ ручной\ работы$


где $T$ время, которое будет затрачено на отрисовку всех зон доставки партнера,
$R$ количество ресторанов,
$SLA$ время на отрисовку одной зоны.

Проблемы ручной отрисовки зон:


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

Baseline


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

При этом оставались недостатки:

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

Через Delivery Club ежедневно проходят сотни тысяч заказов, и артефакты зон доставки могли дать большую нагрузку на колл-центр.



Пробовали алгоритмы на основе Convex Hull и Alpha Shape, с их помощью можно было бы устранить артефакты. Например, используя границы водных массивов как опорные точки для построения форм удавалось обходить реки. Однако всё равно отсутствовал индивидуальный подход к партнерам.

Преимущества технологии H3


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


Источник: eng.uber.com/h3

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

Также стоит отметить, что:

  1. Существует хорошая библиотека для работы с H3, которую и выбрала наша команда в качестве основного инструмента. Библиотека поддерживает многие языки программирования (Python, Go и другие), в которых уже реализованы основные функции для работы с гексагонами.
  2. Наша реляционная аналитическая база Postgres поддерживает работу с нативными функциями H3.
  3. При использовании гексагональной сетки благодаря ряду алгоритмов, работающих с индексами, можно очень быстро получить точную информацию о признаках в соседних ячейках, например, определить вхождение точки в гексагон.
  4. Преимуществом гексагонов является возможность хранить информацию о признаках не только в конкретных точках, но и в целых областях.

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


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


  2. Убираем точки-выбросы. Анализируем все постройки и сразу отбрасываем те, которые нас не интересуют. Например, какие-то мелкие нежилые объекты. Далее с помощью DBSCAN формируем кластеры точек и отбрасываем те, которые для нас не являются важными: например, если кластер находится далеко от ресторана или нам экономически невыгодно доставлять туда.


  3. Далее на основе очищенного набора точек применяем триангуляцию Делоне.


  4. Создаем сетку гексагонов H3.


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


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



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

    • минимизацию времени доставки при фиксированном покрытии;
    • максимизацию охвата пользователей при фиксированном времени доставки.

    Пример функции ошибки для минимизации времени доставки:

    ${L_{min}}_{time}\;=\;min(\sum_{i=1}^n\;({t_{rest}}_i)/n),$


    где $L_{min_ {time}}$ функция ошибки минимизации времени доставки с фиксированным покрытием,
    $t_{rest_ {i}}$ время от ресторана i до клиента,
    $n$ количество клиентов в зоне доставки.

  7. Далее строим временной градиент в получившихся зонах (с очищенными выбросами) и с заранее определенными интервалами (например, по 10-минутным отрезкам пешего пути).



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

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

Внедрение


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

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

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

Но тут пришел COVID-19

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

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

Оценка


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



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

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

$T = 100 * 3,6\ секунды =\ \sim 6\ минут$


Ускорение в 670 раз!

Текущая ситуация и планы


Сервис работает в production. Зоны автоматически строятся по кнопке. Появился более гибкий инструмент для работы со стоимостью доставки для клиентов в зависимости от их удаленности от ресторана. 99,9% ресторанов (изредка ещё нужно ручное вмешательство) перешли на алгоритмические зоны, а наш алгоритм поспособствовал переходу бэкенда на H3.

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

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

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

Всем спасибо!

Подробнее..

Перевод Как культура жителей города влияет на дизайн карт метро Нью-Йорк

13.06.2021 14:11:44 | Автор: admin
Что бы вы порекомендовали тому, кто впервые приезжает в Нью-Йорк? Посетить Центральный парк? Посмотреть шоу на Бродвее? Увидеть Статую Свободы?

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

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

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


Самая новая карта Нью-Йоркского метрополитена

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


Карта метрополитена Стамбула


Карта метрополитена Афин


Карта лондонского метрополитена

Менялась ли карта одной из старейших транспортных систем мира?


Задавшись этим вопросом, я отправилась по линии F в Музей транспорта Нью-Йорка. В нём представлено множество карт метрополитена, начиная от самой первой, спроектированной в 1904 году, до самой последней, выпущенной в 2020 году. Как изменились эти карты за последние 116 лет?


Слева направо: IRT 1904 года, BMT 1924 года, IND 1939 года

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


Карта Нью-Йоркского метрополитена, Массимо Виньелли, 1972 год

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

Почему эта карта не нравилась нью-йоркцам?


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

  • Оттенки цветов линий были близки друг к другу, поэтому их трудно было различать;
  • Водоёмы не были окрашены в синий (они были бежевыми);
  • Парки не были окрашены в зелёный, как в реальной жизни;
  • Станции были расположены рядом друг с другом, даже если переход между ними был невозможен.

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

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

После того, как карту Виньелли перестали использовать, Нью-Йорк вернулся к географической карте, которая используется и сегодня. Майкл Херц спроектировал и разработал карту совместно с Комитетом карт метрополитена, которым руководил Джон Тауранак; Херц сыграл важнейшую роль в создании дизайна новой карты. Он добавил на неё названия улиц и достопримечательностей. Несмотря на критику, которой Херц подверг карту Виньелли, в конечном итоге они с Тауранаком воспользовались восьмицветной палитрой Виньелли.

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

Как культура влияет на ощущение от продукта?


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

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

Какие основные факторы влияют на карту?


  1. Город Нью-Йорк обладает мультикультурной структурой:

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

Также карта отражает связь, которую пассажиры имеют со своими сообществами. Эти сообщества представлены на карте косвенным образом. Это ощущение принадлежности отражено на карте Херца посредством названий районов: Бронкс, Чайна-таун, Мотт-Хейвен. Все они позволяют пассажирам не только понять то, куда они едут, но и увидеть на карте важную часть своей жизни.

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

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


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

2. Город имеет структуру сетки улиц:

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

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


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


New York Times, 2012: Туристы испытывали сложности с сопоставлением своих маршрутов с обозначениями на карте. В эпоху до Google Maps её художественное оформление делало поездки в чреве города более пугающими.

На карте Виньелли казалось, что Колумбус-Серкл на 59-й улице всего в паре кварталов от Центрального парка. На первый взгляд, линии 1, 2 и 3 представляли отдельные станции, однако это не так. Станция всех трёх линий расположена под входом в Центральный парк. Поезда линий 1, 2 и 3 останавливаются на той же станции, что и остальные. Это создавало разрыв между картой метро и реальным миром.

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

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

Как повысить удобство продукта для пользователей?


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

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

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

Карта Нью-Йоркского метрополитена так же знаменита, как и сам город. За последние годы дизайн карты сильно не менялся. Однако были разработаны новые карты, спроектированные в надежде на повышение удобства пользования. Существуют и альтернативные проекты дизайна, нацеленные на редизайн карты подземки.

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

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

Источники


Tylor, E. B. (1889). Primitive culture: Researches into the development of mythology, philosophy, religion, language, art and custom. New York: H. Holt.

Norman, D. A. (2004). Emotional design: Why we love (or hate) everyday things. New York: Basic Books.

Norman, D. A. (2013). The design of everyday things. MIT Press.

Kottak, Conrad Phillip. (2013). Cultural anthropology: appreciating human diversity. New York: McGraw-Hill.

Moalosi, R., Popovic, V. & Hickling-Hudson, A. Culture-orientated product design. Int J Technol Des Educ 20, 175190 (2010). https://doi.org/10.1007/s10798-008-9069-1

Tejada, Soledad O., The Public and the Personal: Mapping the NYC Subway System as an Urban Memoryscape (2020). Library Map Prize. 7. https://elischolar.library.yale.edu/library_map_prize/7

https://www.nngroup.com/articles/ten-usability-heuristics/

https://www.nytimes.com/interactive/2019/12/02/nyregion/nyc-subway-map.html

https://www.oreilly.com/content/redesigning-the-new-york-city/

http://www.culturelookingsideways.com/subway-civility-intro



На правах рекламы


Наша компания предлагает VDS в аренду с Windows или Linux. Не экономим на железе только современное оборудование и одни из лучших дата-центров в России и ЕС.

Подписывайтесь на наш чат в Telegram.

Подробнее..

Перевод 6 причин, по которым вам следовало бы отказаться от гистограмм

12.05.2021 14:11:23 | Автор: admin

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

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

И нетрудно догадаться почему. Гистограммы весьма интуитивно наглядны: любой поймет их с первого взгляда. Более того, они объективно представляют реальность, не так ли? А вот и нет.

Гистограмма может ввести в заблуждение и привести к ошибочным выводам даже на простейшем наборе данных!

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

  1. Они слишком сильно зависят от количества интервалов.

  2. Они слишком сильно зависят от максимума и минимума переменной.

  3. Они не дают возможности заметить значимые значения переменной.

  4. Они не позволяют отличить непрерывные переменные от дискретных.

  5. Они делают сравнение распределений сложным.

  6. Их построение затруднено, если в памяти находятся не все данные.

Ладно, я понял: гистограммы не идеальны. Но есть ли у меня выбор? Конечно есть!

В конце статьи я порекомендую другой график, называемый CDP, который минует эти недостатки.

Итак, что же не так с гистограммой?

1. Она слишком сильно зависит от количества интервалов.

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

Переменная представляет собой максимальную частоту сердечных сокращений (ударов в минуту), полученную у 303 людей во время некоторой физической активности (данные взяты из набора данных UCI по сердечным заболеваниям: источник).

Как изменяется гистограмма при изменении количества интервалов. [Рисунок автора]Как изменяется гистограмма при изменении количества интервалов. [Рисунок автора]

Глядя на верхний левый график (который мы получим по умолчанию в Python и R), у нас сложится впечатление хорошего распределения с одним пиком (модой). Однако если бы мы рассмотрели бы другие варианты гистограммы, мы получили бы совершенно другую картину. Разные гистограммы одних и тех же данных могут привести к противоречивым выводам.

2. Она слишком сильно зависит от максимума и минимума переменной.

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

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

Как меняется гистограмма при изменении максимального значения. [Рисунок автора]Как меняется гистограмма при изменении максимального значения. [Рисунок автора]

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

3. Не дает возможности заметить значимые значения переменной.

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

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

Те же данные, разная ширина интервала. На левом графике невозможно обнаружить высокую концентрацию нулей. [Рисунок автора]Те же данные, разная ширина интервала. На левом графике невозможно обнаружить высокую концентрацию нулей. [Рисунок автора]

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

График справа получен путем сужения интервалов и дает более четкое представление о реальности. Но дело в том, что как бы вы ни сужали интервалы, вы никогда не будете уверены, содержит ли первый интервал только 0 или какие-то другие значения.

4. Не позволяет отличить непрерывные переменные от дискретных.

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

Возьмем переменную Возраст (Age). Вы можете получить Возраст = 49 лет (когда возраст округлен) или Возраст = 49,828884325804246 лет (когда возраст рассчитывается как количество дней с момента рождения, деленное на 365,25). Первая дискретная переменная, вторая непрерывная.

Слева непрерывная переменная. Справа дискретная переменная. Однако на верхних графиках они выглядят одинаково. [Рисунок автора]Слева непрерывная переменная. Справа дискретная переменная. Однако на верхних графиках они выглядят одинаково. [Рисунок автора]

Тот, что слева, непрерывен, а тот, что справа, дискретен. Однако на верхних графиках (по умолчанию в Python) вы не увидите никакой разницы между ними: они выглядят совершенно одинаково.

5. Сложно сравнивать распределения.

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

  • все население (для справки)

  • люди моложе 50 страдающие сердечными заболеваниями

  • люди моложе 50 НЕ страдающие сердечными заболеваниями

  • люди старше 60 лет страдающие сердечными заболеваниями

  • люди старше 60 и НЕ страдающие сердечными заболеваниями.

Вот что мы получили бы в итоге:

Сравнение гистограмм. [Рисунок автора]Сравнение гистограмм. [Рисунок автора]

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

6. Сложно построить, если в памяти находятся не все данные.

Если все ваши данные находятся в Excel, R или Python, построить гистограмму легко: в Excel вам просто нужно кликнуть по иконке гистограммы, в R выполнить команду hist(x), а в Python plt.hist(х).

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

| INTERVAL_LEFT | INTERVAL_RIGHT | COUNT |

|---------------|----------------|---------------|

| 75.0 | 87.0 | 31 |

| 87.0 | 99.0 | 52 |

| 99.0 | 111.0 | 76 |

| ... | ... | ... |

Но получить ее с помощью SQL-запроса не так просто, как кажется. Например, в Google Big Query код будет выглядеть так:

WITHSTATS AS (  SELECT     COUNT(*) AS N,    APPROX_QUANTILES(VARIABLE_NAME, 4) AS QUARTILES  FROM    TABLE_NAME),BIN_WIDTH AS (  SELECT    -- freedman-diaconis formula for calculating the bin width    (QUARTILES[OFFSET(4)]  QUARTILES[OFFSET(0)]) / ROUND((QUARTILES[OFFSET(4)]  QUARTILES[OFFSET(0)]) / (2 * (QUARTILES[OFFSET(3)]  QUARTILES[OFFSET(1)]) / POW(N, 1/3)) + .5) AS FD  FROM     STATS),HIST AS (  SELECT     FLOOR((TABLE_NAME.VARIABLE_NAME  STATS.QUARTILES[OFFSET(0)]) / BIN_WIDTH.FD) AS INTERVAL_ID,    COUNT(*) AS COUNT  FROM     TABLE_NAME,    STATS,    BIN_WIDTH  GROUP BY     1)SELECT   STATS.QUARTILES[OFFSET(0)] + BIN_WIDTH.FD * HIST.INTERVAL_ID AS INTERVAL_LEFT,  STATS.QUARTILES[OFFSET(0)] + BIN_WIDTH.FD * (HIST.INTERVAL_ID + 1) AS INTERVAL_RIGHT,  HIST.COUNTFROM   HIST,   STATS,   BIN_WIDTH

Немного громоздко, не правда ли?

Альтернатива: график кумулятивного распределения.

Узнав 6 причин, по которым гистограмма не является идеальным выбором, возникает естественный вопрос: Есть ли у меня альтернатива? Хорошие новости: существует лучшая альтернатива, которая называется График кумулятивного распределения (Cumulative Distribution Plot - CDP). Я знаю, что это название не такое запоминающееся, но гарантирую, оно того стоит.

График кумулятивного распределения это график квантилей переменной. Другими словами, каждая точка CDP показывает:

  • по оси x: исходное значение переменной (как в гистограмме);

  • по оси y: сколько наблюдений имеют такое же или меньшее значение.

Давайте посмотрим на пример с переменной максимальной частотой пульса.

График кумулятивного распределения максимальной частоты сердечных сокращений. [Рисунок автора]График кумулятивного распределения максимальной частоты сердечных сокращений. [Рисунок автора]

Возьмем точку с координатами x = 140 и y = 90 (30%). По горизонтальной оси вы видите значение переменной: 140 ударов сердца в минуту. По вертикальной оси вы видите количество наблюдений, у которых частота сердцебиение равна или ниже 140 (в данном случае 90 человек, что означает 30% выборки). Следовательно, у 30% нашей выборки максимальная частота сердцебиения составляет 140 или менее ударов в минуту.

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

Вдобавок CDP намного полезнее. Если задуматься, вам часто приходится отвечать на такие вопросы, как у скольких из них от 140 до 160? Или у скольких из них больше 180?. Имея перед глазами CDP, вы можете дать немедленный ответ. С гистограммой это было бы невозможно.

CDP решает все проблемы, которые мы видели выше. Фактически, по сравнению с гистограммой:

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

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

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

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

5. Упрощает сравнение распределений. На одном графике легко сравнить два или более распределения, поскольку это просто кривые, а не области. Кроме того, ось y всегда находится в диапазоне от 0 до 100%, что делает сравнение еще более простым. Для сравнения, это пример, который мы видели выше:

Сравнение распределений в CDP. [Рисунок автора]Сравнение распределений в CDP. [Рисунок автора]

6. Его легко построить, даже если у вас нет всех данных в памяти. Все, что вам нужно, это квантили, которые можно легко получить с помощью SQL:

SELECT   COUNT(*) AS N,  APPROX_QUANTILES(VARIABLE_NAME, 100) AS PERCENTILESFROM  TABLE_NAME

Как построить график кумулятивного распределения в Excel, R, Python

В Excel вам нужно построить два столбца. Первый с 101 числом, равномерно распределенными от 0 до 1. Второй столбец должен содержать процентили, которые могут быть получены по формуле: =PERCENTILE(DATA, FRAC), где DATA - это вектор, содержащий данные, а FRAC - это первый столбец: 0,00, 0,01, 0,02, 0,03,, 0,98, 0,99, 1. Затем вам просто нужно построить график по этим двум столбцам, разместив значения переменной на оси x.

В R это делается в одну строчку:

plot(ecdf(data))

В Python:

from statsmodels.distributions.empirical_distribution import ECDFimport matplotlib.pyplot as pltecdf = ECDF(data)plt.plot(ecdf.x, ecdf.y)

Спасибо за внимание! Надеюсь, эта статья оказалась для вас полезной.

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


Перевод материала подготовлен в рамках онлайн-курса "Machine Learning. Basic". Всех заинтересованных приглашаем на день открытых дверей курса, где можно будет узнать все подробности об обучении и пообщаться с преподавателем.

- Узнать подробнее о курсе "Machine Learning. Basic"

- Смотреть онлайн-встречу "День открытых дверей"

Подробнее..

Погружаемся в статистику вместе с Python. Часть 1. Z-статистика и p-value

14.05.2021 16:12:54 | Автор: admin

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

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

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

  • "Теория вероятности и математическая статистика" Л. Н. Фадеева, А. В. Лебедев;

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

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

Ну и наконец книга "Статистика для всех" Сара Бослаф. Я бы порекомендовал для прочтения только эту книгу, но мне кажется, что если бы я ее читал, абсолютно не зная теорию вероятности и мат. статистику, то я бы очень много вообще не понял, ну или очень долго вникал в текст, формулы и картинки. Там конечно есть разделы с основами для конкретных новичков, но на мне такое не работает. Хотя, более чем вероятно, что этой книги для вас этого окажется более чем достаточно. В книге намеренно нет примеров с кодом, но это скорее плюс чем минус, потому что в книге есть примеры численных расчетов. Так что можно кодить и сравнивать свои результаты с результатами из книги.

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

Вместо того что бы развозить математику и доказательства (кто-то скажет, что это кощунство) я предлагаю немного покодить и пойти по принципу "если что-то непонятно то моделируй и экспериментируй до тех пор пока не станет понятно". Само программирование потребует от нас хотя бы небольших познаний библиотек NumPy, matplotlib, seaborn и SciPy - открытых вкладок с документацией будет более чем достаточно. Итак, начнем.

Пример не для тех кто проголодался

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

import numpy as npfrom scipy import statsimport matplotlib.pyplot as plt#import seaborn as snsfrom pylab import rcParams#sns.set()rcParams['figure.figsize'] = 10, 6%config InlineBackend.figure_format = 'svg'np.random.seed(42)

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

norm_rv = stats.norm(loc=30, scale=5)samples = np.trunc(norm_rv.rvs(365))samples[:30]
array([32., 29., 33., 37., 28., 28., 37., 33., 27., 32., 27., 27., 31.,       20., 21., 27., 24., 31., 25., 22., 37., 28., 30., 22., 27., 30.,       24., 31., 26., 28.])

Конечно же теперь нам интересно выяснить среднее время доставки пиццы и его среднеквадратическое отклонение:

samples.mean(), samples.std()
(29.52054794520548, 4.77410133275075)

Можно сказать, что время доставки пиццы занимает где-то 30\pm5 минут.

А еще, было бы интересно посмотреть на то, как распределены данные:

sns.histplot(x=samples, discrete=True);

Глядя на такой график, мы вполне можем допустить, что время доставки пиццы имеет нормальный закон распределения с параметрами \mu = 30 и \sigma = 5 . Кстати, а почему мы решили, что распределение нормальное? Потому что гистограмма хорошо смотрится на фоне функции распределения плотности вероятности нормального распределения? Если речь идет о визуальном предпочтении, то с таким же успехом мы можем подогнать и нарисовать функции распределений плотности гамма, бета и даже треугольного распределения:

norm_rv = stats.norm(loc=30, scale=5)beta_rv = stats.beta(a=5, b=5, loc=14, scale=32)gamma_rv = stats.gamma(a = 20, loc = 7, scale=1.2)tri_rv = stats.triang(c=0.5, loc=17, scale=26)x = np.linspace(10, 50, 300)sns.lineplot(x = x, y = norm_rv.pdf(x), color='r', label='norm')sns.lineplot(x = x, y = beta_rv.pdf(x), color='g', label='beta')sns.lineplot(x = x, y = gamma_rv.pdf(x), color='k', label='gamma')sns.lineplot(x = x, y = tri_rv.pdf(x), color='b', label='triang')sns.histplot(x=samples, discrete=True, stat='probability',             alpha=0.2);

Распределение точно нормальное?

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

  • на перекрестке пришлось ждать две минуты пока светофор загорится зеленым;

  • ударился ногой и из-за хромоты шел дольше обычного;

  • доставщик оказался скейтбордистом и передвигался быстрее обычного;

  • дорогу перебежала черная кошка и пришлось идти другим более долгим путем;

  • развязались шнурки и пришлось тратить время на их завязывание;

  • развязались шнурки и доставщик упал, поэтому пришлось тратить время на отряхивание грязи и завязывание шнурков.

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

unif_rv = stats.uniform(loc=0, scale=4)unif_rv.rvs()
0.8943833540778106

А завтра, после-завтра и после-после-завтра это время может быть равно:

unif_rv.rvs(size=3)
array([3.85289016, 0.0486179 , 3.87951531])

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

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

  • мог забыть пицу, вернуться за ней, а потом вспомнить, что взял не ту пицу - т.е. вернуться за пицей два раза;

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

Если мы придумали всего 15 случайных переменных: X_{1}, X_{2},.., X_{15} , то можно сказать, что общее время доставки является их суммой и тоже является случайной величиной, которую можно обозначить буквой Y :

Y = X_{1} + X_{2} + .. + X_{15} = \sum_{i = 1}^{15} X_{i}

Но если значения каждой из переменных X_{1}, X_{2},.., X_{15} распределены равномерно, то как будет распределена их сумма - переменная Y ? Давайте сгенерируем 10 тысяч таких сумм и посмотрим на гистограмму:

Y_samples = [unif_rv.rvs(size=15).sum() for i in range(10000)]sns.histplot(x=Y_samples);

Не напоминает ли вам эта картинка тот самый колокол, того самого нормального распределения? Если напоминает, то вы только, что поняли смысл центральной предельной теоремы: распределение суммы случайных переменных стремится к нормальному распределению при увеличении количества слагаемых в этой сумме.

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

  • на перекрестке пришлось ждать две минуты пока светофор загорится зеленым. Но ждать придется только если мы подошли к светофору, который уже горит красным;

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

  • доставщик оказался скейтбордистом и передвигался быстрее обычного. При условии что он не забыл скейтборд дома;

  • дорогу перебежала черная кошка и пришлось идти другим более долгим путем. Существует ненулевая вероятность того, что доставщик не является фанатом бабы Нины;

  • развязались шнурки и пришлось тратить время на их завязывание. Как часто развязываются шнурки?;

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

А то, что каждая из переменных X_{1}, X_{2},.., X_{15} носит условный характер означает, что они могут входить в сумму в самых разных комбинациях. Например, сегодня время доставки задавалось, как:

Y = X_{3} + X_{5}+ X_{7} + X_{11}

А завтра это время может задаваться как:

Y = X_{1} + X_{3}+ X_{9} + X_{10}+ X_{13} + X_{15}

Будет ли теперь Y распределена нормально? Учитывая что сумма нормально распределенных величин тоже имеет нормальное распределение, то можно дать утвердительный ответ. Именно поэтому, когда мы взглянули на распределение 365-и значений времени доставки, мы практически сразу решили, что перед нами нормальное распределение, даже несмотря на то что оно вовсе не похоже на идеальный колокол.

Z-значения

Допустим, по прошествии года у нас появился новый сосед и он так же как и мы решил ежедневно заказывать пицу. И вот по прошествии трех дней мы наблюдаем, как этот сосед обвиняет доставщика в слишком долгом ожидании заказа. Мы решаем поддержать доставщика и говорим, что в среднем время доставки занимает 30\pm5 минут, на что наш сосед отвечает, что все три раза он ждал больше 40 минут, а это всяко больше 35 минут.

Почему наш сосед так уверен в долгой доставке? И вообще, оправдана ли его уверенность? Очевидно он, как и некоторые люди, думает, что 30\pm5 минут означает, что доставка может длиться 27, 31, даже 35 минут, но никак не 23 или 38 минут. Однако, мы заказывали пицу 365 раз и знаем, что доставка может длиться и 20 и даже 45 минут. А фраза 30\pm5 минут, означает лишь то, что какая-то значительная часть доставок бодет занимать от 25 до 35 минут. Зная параметры распределения, мы даже можем смоделировать несколько тысяч доставок и прикинуть величину этой части:

N = 5000t_data = norm_rv.rvs(N)t_data[(25 < t_data) & (t_data < 35)].size/N
0.7008

Где-то две трети значений укладываются в интервал от 25 до 35 минут. А сколько значений будет превосходить 40 минут?

t_data[t_data > 40].size/N
0.0248

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

0.02**3
8.000000000000001e-06

Конечно, компьютерное моделирование - это очень хорошо. Но в данном случае лучше воспользоваться так называемыми Z-значениями.

Вычислить Z-значение можно по следующей формуле:

Z = \frac{y - \mu}{\sigma}

Где y - это время доставки, т.е. какое-то конкретное значение случайной переменной Y , которая имеет нормальное распределение, а \mu и \sigma это мат. ожидание и среднеквадратическое отклонение, т.е. параметры распределения, в нашем случае они равны 30 и 5 минут соответственно. Давайте рассчитаем Z-значение для сорока минут:

Z = \frac{40 - 30}{5} = 2

Что мы сейчас сделали? В числителе мы вычислили на какую величину наше время доставки отличается от среднего времени доставки, а далее мы просто поделили это значение на стандартное отклонение времени доставки. Но как интерпретировать данный результат и зачем вообще использовать Z-значение? Что бы понять это придется немного "порисовать":

fig, ax = plt.subplots()x = np.linspace(norm_rv.ppf(0.001), norm_rv.ppf(0.999), 200)ax.vlines(40.5, 0, 0.1, color='k', lw=2)sns.lineplot(x=x, y=norm_rv.pdf(x), color='r', lw=3)sns.histplot(x=samples, stat='probability', discrete=True);

На данном графике мы изобразили гистограмму наших данных samples, но теперь высота прямоугольников показывает не количество вхождений каждого значения в выборку, а вероятность их появления в выборке. Красной линией мы нарисовали функцию распределения плотности вероятности значений времени доставки. Выше мы пытались экспериментальным путем определить какую долю составляют значения времени больше 40 минут. И данный график должен натолкнуть нас на мысль о том, что есть два способа сделать это. Первый - экспериментальный, т.е. мы моделируем, скажем, 5000 доставок, строим гистограмму и находим сумму высот прямоугольников, расположенных правее черной линии:

np.random.seed(42)N = 10000values = np.trunc(norm_rv.rvs(N))fig, ax = plt.subplots()v_le_41 = np.histogram(values, np.arange(9.5, 41.5))v_ge_40 = np.histogram(values, np.arange(40.5, 51.5))ax.bar(np.arange(10, 41), v_le_41[0]/N, color='0.8')ax.bar(np.arange(41, 51), v_ge_40[0]/N)p = np.sum(v_ge_40[0]/N)ax.set_title('P(Y>40min) = {:.3}'.format(p))ax.vlines(40.5, 0, 0.08, color='k', lw=1);

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

fig, ax = plt.subplots()x = np.linspace(norm_rv.ppf(0.0001), norm_rv.ppf(0.9999), 300)ax.plot(x, norm_rv.pdf(x))ax.fill_between(x[x>41], norm_rv.pdf(x[x>41]), np.zeros(len(x[x>41])))p = 1 - norm_rv.cdf(41)ax.set_title('P(Y>40min) = {:.3}'.format(p))ax.hlines(0, 10, 50, lw=1, color='k')ax.vlines(41, 0, 0.08, color='k', lw=1);

Значения вероятности в первом и втором случае получились довольно близкими. Но надо оговориться, что в первом случае мы использовали дискретную случайную величину, а во втором - непрерывную, поэтому код выглядит несколько "читерским", как буд-то бы все намеренно подгонялось так, что бы два результата были похожи. Конечно, возможно я что-то перемудрил, но вроде бы все верно. В первом случае мы считаем сколько раз доставка составляла более 40 минут, т.е. считали сколько доставок длились 41, 42, 43 и т.д. минут. Во втором случае мы просто считаем сколько непрерывных величин оказалось правее значения 41.0. В принципе, модуль SciPy позволяет создавать собственные распределения случайных переменных, в руководстве даже имеется пример, как создать дискретно-нормальное распределение. Но если задуматься, то наличие дискретных значений и использование непрерывных функций распределения, это своего рода - неизбежность, ведь подавляющее большинство измерений носит дискретный характер: рост, вес, доход и т.д. Так что далее мы не будем больше оговариваться по этому поводу, а лучше вернемся к Z-значениям.

Итак, допустим мы оказались в средиземье и каким-то образом выяснили что рост хобитов и гномов в сантиметрах распределен как N(91;8^{2}) и N(134;6^{2}) . Если рост Фродо равен 99 сантиметрам а рост Гимли 143 сантиметра, то как понять чей рост более типичен среди своих народов? Что бы выяснить это мы можем изобразить функцию распределения плотности вероятности для каждого народа с отмеченными значениями, а заодно определить долю тех, кто превышает эти значения:

fig, ax = plt.subplots(nrows=1, ncols=2, figsize = (12, 5))nrv_hobbit = stats.norm(91, 8)nrv_gnome = stats.norm(134, 6)for i, (func, h) in enumerate(zip((nrv_hobbit, nrv_gnome), (99, 143))):    x = np.linspace(func.ppf(0.0001), func.ppf(0.9999), 300)    ax[i].plot(x, func.pdf(x))    ax[i].fill_between(x[x>h], func.pdf(x[x>h]), np.zeros(len(x[x>h])))    p = 1 - func.cdf(h)    ax[i].set_title('P(H > {} см) = {:.3}'.format(h, p))    ax[i].hlines(0, func.ppf(0.0001), func.ppf(0.9999), lw=1, color='k')    ax[i].vlines(h, 0, func.pdf(h), color='k', lw=2)    ax[i].set_ylim(0, 0.07)fig.suptitle('Хобиты слева, гномы справа', fontsize = 20);

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

Выполнить сравнение типичности гораздо проще и нагляднее если воспользоваться Z-значениями:

\begin{align*} & Z_{\textrm{Frodo}} = \frac{99 - 91}{8} = 1 \\ & \\ & Z_{\textrm{Gimli}} = \frac{143 - 134}{6} = 1.5 \end{align*}
fig, ax = plt.subplots()N_rv = stats.norm()x = np.linspace(N_rv.ppf(0.0001), N_rv.ppf(0.9999), 300)ax.plot(x, N_rv.pdf(x))ax.hlines(0, -4, 4, lw=1, color='k')ax.vlines(1, 0, 0.4, color='r', lw=2, label='Фродо')ax.vlines(1.5, 0, 0.4, color='g', lw=2, label='Гимли')ax.legend();

Огромным преимуществом Z-значений является то, что они "стандартизированы", т.е. преобразованы так словно они взяты из стандартного нормального распределения N(0;1) , именно поэтому два Z-значения нарисованы на фоне единственной кривой. Однако, в общем случае, даже само рисование графиков вовсе не обязательно, потому что меньшие по модулю Z-значения обладают большей частотой появления. Одной записи:

\left | Z_{\textrm{Frodo}} \right | < \left | Z_{\textrm{Gimli}} \right |

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

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

Z = \frac{y - \mu}{\sigma}

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

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

Z-статистика

Теперь давайте снова вернемся к нашему соседу, который возмущен слишком долго доставкой пицы. Выше мы вычислили Z-значение для 40 минут:

Z = \frac{40 - 30}{5} = 2

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

fig, ax = plt.subplots()N_rv = stats.norm()x = np.linspace(N_rv.ppf(0.0001), N_rv.ppf(0.9999), 300)ax.plot(x, N_rv.pdf(x))ax.hlines(0, -4, 4, lw=1, color='k')ax.vlines(2, 0, 0.4, color='g', lw=2);

Это значение находится справа от вершины на расстоянии 2\sigma , что довольно много. Однако, сосед утверждает, что ждал пицу больше сорока минут три дня подряд. Возможно он и не вычислял среднее время ожидания своей пицы, но три значения - это уже статистика!

Характеристики генеральной совокупности называют параметрами, а характеристики выборки - статистиками. Замеряя время доставки на протяжении 365-и дней, мы сделали вывод о параметрах генеральной совокупности, т.е. времени всех возможных доставок пицы, решив что эти значения берутся из N(30;5^{2}) . А зная распределение мы можем поэкспериментировать. Например, наш сосед сделал всего три заказа, и по его ощущениям среднее время доставки, было больше сорока минут. А что если мы повторим этот эксперимент 5000 раз и посмотрим на то как распределено среднее трех заказов:

sns.histplot(np.trunc(norm_rv.rvs(size=(N, 3))).mean(axis=1),             bins=np.arange(19, 42));

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

  • либо сосед врет;

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

Вычисленное выше Z-значение для сорока минут (Z = 2), позволяет оценить долю (вероятность появления) значений больших сорока:

1 - N_rv.cdf(2)
0.02275013194817921

Поэтому не удивительно, что вероятность получить среднее время трех доставок большее 40 минут исчезающе мала:

(1 - N_rv.cdf(2))**3
1.1774751750476762e-05

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

Z = \frac{\bar{x} - \mu}{\frac{\sigma}{\sqrt{n}}}

где \bar{x} - это среднее значение для нашей выборки, \mu и \sigma среднее значение и стандартное отклонение для генеральной совокупности, а n - размер выборки.

Давайте предположим что мы сделали три заказа и среднее значение оказалось равным 35 минутам, тогда Z-статистика будет вычисляться так:

Z = \frac{35 - 30}{\frac{5}{\sqrt{3}}} \approx 1.73

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

\begin{bmatrix}\mu - \left | \mu - \bar{x} \right |; \mu + \left | \mu - \bar{x} \right |\end{bmatrix}

который в нашем случае выглядит как [25; 35] минут. Как и ранее мы можем найти данную вероятность с помощью моделирования:

N = 10000means = np.trunc(norm_rv.rvs(size=(N, 3))).mean(axis=1)means[(means>=25)&(means<=35)].size/N
0.9241
N = 10000fig, ax = plt.subplots()means = np.trunc(norm_rv.rvs(size=(N, 3))).mean(axis=1)h = np.histogram(means, np.arange(19, 41))ax.bar(np.arange(20, 25), h[0][0:5]/N, color='0.5')ax.bar(np.arange(25, 36), h[0][5:16]/N)ax.bar(np.arange(36, 41), h[0][16:22]/N, color='0.5')p = np.sum(h[0][6:16]/N)ax.set_title('P(25min < Y < 35min) = {:.3}'.format(p))ax.vlines([24.5 ,35.5], 0, 0.08, color='k', lw=1);

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

x, n, mu, sigma = 35, 3, 30, 5z = abs((x - mu)/(sigma/n**0.5))N_rv = stats.norm()fig, ax = plt.subplots()x = np.linspace(N_rv.ppf(0.0001), N_rv.ppf(0.9999), 300)ax.plot(x, N_rv.pdf(x))ax.hlines(0, x.min(), x.max(), lw=1, color='k')ax.vlines([-z, z], 0, 0.4, color='g', lw=2)x_z = x[(x>-z) & (x<z)] # & (x<z)ax.fill_between(x_z, N_rv.pdf(x_z), np.zeros(len(x_z)), alpha=0.3)p = N_rv.cdf(z) - N_rv.cdf(-z)ax.set_title('P({:.3} < z < {:.3}) = {:.3}'.format(-z, z, p));

Обратите внимание, на то что Z-статистика зависит как от среднего выборки - \bar{x} , так и от величины выборки - n . Если мы закажем пицу 5, 30 или 100 раз, то какова вероятность того, что среднее значение времени доставок окажется в интервале [29; 31]? Давайте взглянем:

fig, ax = plt.subplots(nrows=1, ncols=3, figsize = (12, 4))for i, n in enumerate([5, 30, 100]):    x, mu, sigma = 31, 30, 5    z = abs((x - mu)/(sigma/n**0.5))    N_rv = stats.norm()    x = np.linspace(N_rv.ppf(0.0001), N_rv.ppf(0.9999), 300)    ax[i].plot(x, N_rv.pdf(x))    ax[i].hlines(0, x.min(), x.max(), lw=1, color='k')    ax[i].vlines([-z, z], 0, 0.4, color='g', lw=2)    x_z = x[(x>-z) & (x<z)] # & (x<z)    ax[i].fill_between(x_z, N_rv.pdf(x_z), np.zeros(len(x_z)), alpha=0.3)    p = N_rv.cdf(z) - N_rv.cdf(-z)    ax[i].set_title('n = {}, z = {:.3}, p = {:.3}'.format(n, z, p));fig.suptitle(r'Z-статистики для $\bar{x} = 31$', fontsize = 20, y=1.02);

При 5 заказах среднее выборки попадает в интервал [29;31] скорее случайно, чем систематически. При 30 заказх около четверти средних значений так и не войдут в заданный интервал. И только при сотне заказов мы можем быть более-менее уверены в том что отклонение среднего выборки от среднего генеральной совокупности не будет больше 1 минуты.

С другой стороны, мы можем рассуждать и по другому: если среднее генеральной совокупности равно 30 минутам, то какова вероятность получить среднее выборки равное 31 минуте если мы сделаем 5, 30 или 100 заказов? Очевидно, что при n=5 среднее выборки, может отклоняться очень сильно, следовательно, вероятность получить \bar{x} = 31 минуте очень высока. Но при n=100 среднее выборки практически не отклоняется от среднего генеральной совокупности, поэтому получить случайным образом \bar{x} = 31 при n=100 практически невозможно. Что это значит? А это значит, что если мы сделали 100 заказов и получили среднее время доставки равное 31 минуте, то скорее всего мы ошибаемся насчет того, что среднее генеральной совокупности равно 30 минутам.

Но как быть с нашим соседом, ведь он сделал всего 3 заказа? Неужели он действительно прав? Судя по всему - да. Даже если среднее время ожидания составило 40 минут, то Z-статистика будет равна 3.81, а площадь под кривой будет практически равна 1:

x, n, mu, sigma = 41, 3, 30, 5z = abs((x - mu)/(sigma/3**0.5))N_rv = stats.norm()fig, ax = plt.subplots()x = np.linspace(N_rv.ppf(1e-5), N_rv.ppf(1-1e-5), 300)ax.plot(x, N_rv.pdf(x))ax.hlines(0, x.min(), x.max(), lw=1, color='k')ax.vlines([-z, z], 0, 0.4, color='g', lw=2)x_z = x[(x>-z) & (x<z)] # & (x<z)ax.fill_between(x_z, N_rv.pdf(x_z), np.zeros(len(x_z)), alpha=0.3)p = N_rv.cdf(z) - N_rv.cdf(-z)ax.set_title('P({:.3} < z < {:.3}) = {:.3}'.format(-z, z, p));

В свою очередь, это значит, что вероятность случайного отклонения среднего выборки из трех значений, взятых из генеральной совокупности с \mu=30 и \sigma=5 более чем на 10 минут исчезающе мала. В этом случае мы можем заключить только следующее:

  • либо наш сосед просто чрезвычайно "везучий" человек;

  • либо доставка пицы и вправду теперь выполняется дольше обычного.

Что из этих пунктов является наиболее вероятным? Скорее всего сосед прав насчет долгой доставки.

p-value

Мы смогли убедиться в том что Z-статистика позволяет оценить вероятность того, что среднее выборки \bar{x} размером n , взятой из генеральной совокупности попадет в заданный интервал значений. Это удобно тем, что позволяет сделать вывод о случайности полученного \bar{x} . Чем меньше модуль значения Z-статистики, тем меньше достоверность среднего. Например выше мы видели, что вероятность попадания \bar{x} при n=5 в интервал [29;31] составляет всего около 0.35. В то время как вероятность непопадания в заданный интервал равна 10.35=0.65. Поэтому мы и сделали вывод о том, что значение \bar{x}=31 при n=5 скорее обусловлено случайностью, чем какими-то объективными причинами.

Фактически, значение вероятности 0.65, полученное выше - это и есть то самое p-value которое как раз и показывает вероятность случайного выхода за границы интервала, задаваемого значением Z-статистики, что можно изобразить следующим образом:

x, n, mu, sigma = 31, 5, 30, 5z = abs((x - mu)/(sigma/n**0.5))N_rv = stats.norm()fig, ax = plt.subplots()x = np.linspace(N_rv.ppf(0.0001), N_rv.ppf(0.9999), 300)ax.plot(x, N_rv.pdf(x))ax.hlines(0, x.min(), x.max(), lw=1, color='k')ax.vlines([-z, z], 0, 0.4, color='g', lw=2)x_le_z, x_ge_z = x[x<-z], x[x>z]ax.fill_between(x_le_z, N_rv.pdf(x_le_z), np.zeros(len(x_le_z)), alpha=0.3, color='b')ax.fill_between(x_ge_z, N_rv.pdf(x_ge_z), np.zeros(len(x_ge_z)), alpha=0.3, color='b')p = N_rv.cdf(z) - N_rv.cdf(-z)ax.set_title('P({:.3} < z < {:.3}) = {:.3}'.format(-z, z, p));

Чем меньше p-value тем меньше вероятность того, что среднее выборки получено случайно. При этом p-value напрямую связано с двусторонними гипотезами, т.е. гипотезами о попадании величины в заданный интервал. Если мы получили какие-то результаты, но p-value оказалось довольно большим, то вряд ли эти результаты могут считаться значимыми. Причем, традиционно, уровень значимости, обозначаемый буквой \alpha равен 0.05, а это означает, что для подтверждения значимости результатов p-value должно быть меньше этого уровня. Однако, стоит обязательно отметить, что традиционный уровень значимости \alpha=0.05 может быть непригоден в некоторых областях исследований. Например, в сфере образования наверняка можно обойтись \alpha = 0.1 , а вот в квантовой физике, запросто придется снизить этот уровень до 5 сигм, т.е. \alpha будет равна:

1 - (N_rv.cdf(5) - N_rv.cdf(-5))
5.733031438470704e-07

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

С другой стороны, что бы подтвердить небольшие отклонения от среднего генеральной совокупности, придется очень сильно увеличивать размер выборки. Так, например, если мы хотим заявить с уровнем значимости \alpha=0.05 , что среднее время доставки пицы равно 31 минуте, а не 30 как считалось ранее, то придется сделать не менее 100 заказов.

На последок

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

Это кажется не очень правдоподобным, но давайте взглянем. Сгенерируем 1000 значений из равномерного, экспоненциального и Лапласова распределения, а затем, последовательно, для каждого распределения построим kde-графики распределений среднего значения выборок разного размера:

Код для гифки
import matplotlib.animation as animationfig, axes = plt.subplots(nrows=2, ncols=3, figsize = (12, 7))unif_rv = stats.uniform(loc=10, scale=10)exp_rv = stats.expon(loc=10, scale=1.5)lapl_rv = stats.laplace(loc=15)np.random.seed(42)unif_data = unif_rv.rvs(size=1000)exp_data = exp_rv.rvs(size = 1000)lapl_data = lapl_rv.rvs(size=1000)title = ['Uniform', 'Exponential', 'Laplace']data = [unif_data, exp_data, lapl_data]y_max = [60, 310, 310]n = [3, 10, 30, 50]for i, ax in enumerate(axes[0]):    sns.histplot(data[i], bins=20, alpha=0.4, ax=ax)    ax.vlines(data[i].mean(), 0, y_max[i], color='r')    ax.set_xlim(10, 20)    ax.set_title(title[i])def animate(i):    for ax in axes[1]:        ax.clear()    for j in range(3):        rand_idx = np.random.randint(0, 1000, size=(1000, n[i]))        means = data[j][rand_idx].mean(axis=1)        sns.kdeplot(means, ax=axes[1][j])        axes[1][j].vlines(means.mean(), 0, 2, color='r')        axes[1][j].set_xlim(10, 20)        axes[1][j].set_ylim(0, 2.1)        axes[1][j].set_title('n = ' + str(n[i]))    fig.set_figwidth(15)    fig.set_figheight(8)    return axes[1][0], axes[1][1],axes[1][2]    dist_animation = animation.FuncAnimation(fig,                                       animate,                                       frames=np.arange(4),                                      interval = 200,                                      repeat = False)#  Сохраняем анимацию в виде gif файла:dist_animation.save('dist_means.gif',                 writer='imagemagick',                  fps=1)

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

В общем, спасибо за внимание. Жму F5 и жду ваших и комментариев.

Подробнее..

Оценка кредитного портфеля на R

19.05.2021 18:19:31 | Автор: admin

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


Ниже приведен код на R с прототипом подхода. Не более одного экрана кода на прототип и никаких циклов (закладные для производительности и читаемости).


Декомпозиция


Поскольку мы делаем все с чистого листа, то задачу разбиваем на три шага:


  1. Формирование тестовых данных.
  2. Расчет даты погашения каждого займа.
  3. Расчет и визуализация динамики для заданного временнОго окна.

Допущения и положения для прототипа:


  1. Гранулярность до даты. В одну дату только одна транзакция. Если в один день будет несколько транзакций, то надо будет их порядок устанавливать (для соблюдения принципа FIFO). Можно использовать доп. индексы, можно использовать unixtimestamp, можно еще что-либо придумывать. Для прототипа это несущественно.
  2. Явных циклов for быть не должно. Лишних копирований быть не должно. Фокус на минимальное потребление памяти и максимальную производительность.
  3. Будем рассматривать следующие группы задержек: "< 0", "0-30", "31-60", "61-90", "90+".

Шаг 1.


Генерируем датасет. Просто тестовый датасет. Для каждого пользователя сформируем ~ по 10 записей. Для расчетов полагаем, что займ положительное значение, погашение отрицательное. И весь жизненный цикл для каждого пользователя должен начинаться с займа.


Генерация датасета
library(tidyverse)library(lubridate)library(magrittr)library(tictoc)library(data.table)total_users <- 100events_dt <- tibble(  date = sample(    seq.Date(as.Date("2021-01-01"), as.Date("2021-04-30"), by = "1 day"),    total_users * 10,    replace = TRUE)  ) %>%  # сделаем суммы кратными 50 р.  mutate(amount = (runif(n(), -2000, 1000)) %/% 50 * 50) %>%  # нашпигуем идентификаторами пользователей  mutate(user_id = sample(!!total_users, n(), replace = TRUE)) %>%  setDT(key = "date") %>%  # первая запись должна быть займом  .[.[, .I[1L], by = user_id]$V1, amount := abs(amount)] %>%  # для простоты оставим только одну операцию в день,   # иначе нельзя порядок определить и гранулярность до секунд надо спускать  # либо вводить порядковый номер займа и погашения  unique(by = c("user_id", "date"))

Шаг 2. Расчитываем даты погашения каждого займа


data.table позволяет изменять объекты по ссылке даже внутри функций, будем этим активно пользоваться.


Расчет даты погашения
# инициализируем аккумуляторaccu_dt <- events_dt[amount < 0, .(accu = cumsum(amount), date), by = user_id]ff <- function(dt){  # на вход получаем матрицу пользователей и их платежей на заданную дату  # затягиваем суммы займов  accu_dt[dt, amount := i.amount, on = "user_id"]  accu_dt[is.na(amount) == FALSE, accu := accu + amount][accu > 0, accu := NA, by = user_id]  calc_dt <- accu_dt[!is.na(accu), head(date, 1), by = user_id]  # нанизываем обратно на входной data.frame, сохраняя порядок следования  calc_dt[dt, on = "user_id"]$V1}repay_dt <- events_dt[amount > 0] %>%  .[, repayment_date := ff(.SD), by = date] %>%  .[order(user_id, date)]

Шаг 3. Считаем динамику задолженности за период


Расчет динамики
calcDebt <- function(report_date){  as_tibble(repay_dt) %>%    # выкидываем все, что уже погашено на дату отчета    filter(is.na(repayment_date) | repayment_date > !! report_date) %>%    mutate(delay = as.numeric(!!report_date - date)) %>%    # размечаем просрочки    mutate(tag = santoku::chop(delay, breaks = c(0, 31, 61, 90),                               labels = c("< 0", "0-30", "31-60", "61-90", "90+"),                               extend = TRUE, drop = FALSE)) %>%    # делаем сводку    group_by(tag) %>%    summarise(amount = sum(amount)) %>%    mutate_at("tag", as.character)}# Устанавливаем окно наблюденияdf <- seq.Date(as.Date("2021-04-01"), as.Date("2021-04-30"), by = "1 day") %>%  tibble(date = ., tbl = purrr::map(., calcDebt)) %>%  unnest(tbl)# строим графикggplot(df, aes(date, amount, colour = tag)) +  geom_point(alpha = 0.5, size = 3) +  geom_line() +  ggthemes::scale_colour_tableau("Tableau 10") +  theme_minimal()

Можем получить примерно такую картинку.


Один экран кода, как и требовалось.


Предыдущая публикация Storytelling R отчет против BI, прагматичный подход.

Подробнее..

Перевод Как создавать интерактивные линейные графики на Pandas и Altair

21.05.2021 20:16:58 | Автор: admin

Линейный график является неотъемлемой частью анализа данных. Он даёт нам представление о том, как величина изменяется при последовательных измерениях. В случае работы с временными рядами важность линейных графиков становится решающей. Тренд [направление], сезонность и корреляция вот некоторые характеристики, которые можно наблюдать на аккуратно сгенерированных линейных графиках. В этой статье мы будем создавать интерактивные линейные графики с помощью двух библиотек Python Pandas и Altair.

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


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

Давайте начнём с генерации данных. Типичный пример использования линейных графиков это анализ цен на акции. Один из самых простых способов получения данных о ценах на акции предоставляет библиотека pandas-datareader. Сначала нам нужно импортировать её вместе с уже установленной в Google Colab Pandas.

import pandas as pdfrom pandas_datareader import data

Мы получим цены на акции трёх компаний за один год. Необходимо указать начальную дату, конечную дату и источник данных:

start = '2020-1-1'end = '2020-12-31'source = 'yahoo'

Необходимо знать ещё одну деталь название акции:

apple = data.DataReader("AAPL", start=start ,end=end, data_source=source).reset_index()[["Date", "Close"]]ibm = data.DataReader("IBM", start=start ,end=end, data_source=source).reset_index()[["Date", "Close"]]microsoft = data.DataReader("MSFT", start=start ,end=end, data_source=source).reset_index()[["Date", "Close"]]

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

apple["Stock"] = "apple"ibm["Stock"] = "ibm"microsoft["Stock"] = "msft"stocks["Month"] = stocks.Date.dt.monthstocks = pd.concat([apple, ibm, microsoft])

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

Altair

Altair это библиотека статистической визуализации на Python. Как мы увидим в примерах, её синтаксис чист и прост для понимания. Также с помощью Altair очень просто создавать интерактивные визуализации. Я кратко объясню структуру Altair, а затем сосредоточусь на создании интерактивных линейных графиков. Если вы новичок в Altair, вот туториал по Altair в виде серии из 4 частей:

Список частей

Вот простой линейный график без какой-либо интерактивности:

alt.Chart(stocks).mark_line().encode(   x="Date",   y="Close",   color="Stock").properties(   height=300, width=500)

Основная структура начинается с объекта диаграммы верхнего уровня. Данные могут иметь формат фрейма данных Pandas, или можно написать строку с URL, указывающим на файл JSON или CSV. Затем указывается тип визуализации (mark_circle, mark_line и так далее).

Функция encode указывает Altair, что нужно построить в заданном фрейме данных. Таким образом, всё, что мы пишем в функции encode, должно быть связано с данными. Различать акции мы будем при помощи параметра color, который аналогичен параметру hue в Seaborn. Наконец, с помощью функции properties задаются определённые свойства графика.

Выделение в Altair один из способов интерактивного взаимодействия с пользователем, оно перехватывает действия пользователя.

selection = alt.selection_multi(fields=["Stock"], bind="legend")alt.Chart(stocks).mark_line().encode(   x="Date",   y="Close",   color="Stock",   opacity=alt.condition(selection, alt.value(1), alt.value(0.1))).properties(   height=300, width=500).add_selection(   selection)

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

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

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

hover = alt.selection(   type="single", on="mouseover", fields=["Stock"], nearest=True)

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

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

# line plotlineplot = alt.Chart(stocks).mark_line().encode(   x="Date:T",   y="Close:Q",   color="Stock:N",)# nearest pointpoint = lineplot.mark_circle().encode(   opacity=alt.value(0)).add_selection(hover)# highlightsingleline = lineplot.mark_line().encode(   size=alt.condition(~hover, alt.value(0.5), alt.value(3)))

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

point + singleline

На первом изображении показан исходный, или необработанный, график. На втором обновлённая версия в момент наведения курсора.

Заключение

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

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

Узнайте, как прокачаться и в других специальностях или освоить их с нуля:

Другие профессии и курсы
Подробнее..

Погружаемся в статистику вместе с Python. Часть 2. Распределение Стьюдента

25.05.2021 06:09:43 | Автор: admin

Доброго времени суток, хабраледи и хабраджентельмены! В этой статье мы продолжим погружение в статистику вместе с Python. Если кто пропустил начало погружения, то вот ссылка на первую часть. Ну, а если нет, то я по-прежнему рекомендую держать под рукой открытую книгу Сары Бослаф "Статистика для всех". Так же рекомендую запустить блокнот, чтобы поэкспериментировать с кодом и графиками.

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

Зачем выдумывать новые распределения?

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

import numpy as npfrom scipy import statsimport matplotlib.pyplot as pltimport seaborn as snsfrom pylab import rcParamssns.set()rcParams['figure.figsize'] = 10, 6%config InlineBackend.figure_format = 'svg'np.random.seed(42)

А вот теперь давайте представим, что мы пивовары, придумавшие рецепт и технологию варки классного пива. Продажи идут хорошо, да и в целом все неплохо, но все-таки очень интересно, как потребители оценивают качество нашего пива. Мы решили опросить 1000 человек и оценить качество напитка по 100-бальной шкале в сравнении с другими сортами пива, которые им доводилось пробовать ранее. Данные опроса могли бы выглядеть так:

gen_pop = np.trunc(stats.norm.rvs(loc=80, scale=5, size=1000))gen_pop[gen_pop>100]=100print(f'mean = {gen_pop.mean():.3}')print(f'std = {gen_pop.std():.3}')
mean = 79.5std = 4.95

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

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

[89, 99, 93, 84, 79, 61, 82, 81, 87, 82]

Теперь мы можем вычислить Z-статистику по уже известной нам формуле:

Z = \frac{\bar{x} - \mu}{\frac{\sigma}{\sqrt{n}}}

где \bar{x} - это среднее значение для нашей выборки, \mu и \sigma среднее значение и стандартное отклонение для генеральной совокупности, а n - размер выборки. Вот среднее нашей выборки:

sample = np.array([89,99,93,84,79,61,82,81,87,82])sample.mean()
83.7

Вот значение Z-статистики:

z = 10**0.5*(sample.mean()-80)/5z
2.340085468524603

И еще обязательно вычисляем p-value:

1 - (stats.norm.cdf(z) - stats.norm.cdf(-z))
0.019279327322753836

Для нас, как для пивоваров, это крайне интересные цифры и вот почему: Z-значение отстоит от 0 более чем на 2 сигмы, т.е. отклонение среднего балла 10 респондентов находится очень далеко от среднего бала генеральной совокупности, а вероятность того, что это отклонение произошло случайно, всего около 0.02. Другими словами, если взять 10 человек из совокупности людей, которые оценивают наше "старое" пиво как N(80, 5^{2}) , то вероятность того, что эти 10 человек могут случайно дать среднюю оценку "новому" пиву равную 83.7 баллам составляет всего около 2%. Скорее всего, если мы проведем полномасштабный опрос, то увидим, что потребители не просто не заметят изменения качества пива, а даже отметят его небольшое улучшение. Круто.

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

sample.std(ddof=1)
10.055954565441423
Врезка по поводу ddof в std

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

\begin{align*} & \sigma={\sqrt {{\frac {1}{n}}\sum _{i=1}^{n}\left(x_{i}-{\mu}\right)^{2}}} \\ & \\ & s={\sqrt {{\frac {1}{n - 1}}\sum _{i=1}^{n}\left(x_{i}-{\bar {x}}\right)^{2}}} \end{align*}

Как видите, отличия не очень сильные и заключаются только в том, что при расчете \sigma мы используем мат. ожидание генеральной совокупности - \mu и делим на n , а при вычислении s используем среднее значение выборки и делим на n - 1 . Зачем делить на n - 1 вместо n ? Все дело в том, что при расчете отклонения выборки вместо \mu используется \bar{x} из-за чего значения s получаются несколько ниже, чем реальное значение \sigma для генеральной совокупности. Поэтому знаменатель n - 1 нужен для того, чтобы несколько подправить результат - пододвинуть его чуть поближе к истинному значению.

Для того, чтобы учитывать поправку в методе std() пакета NumPy есть параметр ddof, который по умолчанию равен 0, но если мы применяем std() к небольшой выборке, то необходимо явно указать, что ddof=1. Влияние данного параметра на расчеты может быть очень значимым. Например, если мы возьмем 10000 выборок по 10 элементов из N(80, 5^{2})и построим распределение стандартных отклонений, то увидим, что при ddof=0 пик распределения находится левее реального значения стандартного отклонения. При ddof=1 пик распределения так же расположен левее реального значения, но все-таки ближе к нему, чем при ddof=0:

fig, ax = plt.subplots(nrows=1, ncols=2, figsize = (12, 5))for i in [0, 1]:    deviations = np.std(stats.norm.rvs(80, 5, (10000, 10)), axis=1, ddof=i)    sns.histplot(x=deviations ,stat='probability', ax=ax[i])    ax[i].vlines(5, 0, 0.06, color='r', lw=2)    ax[i].set_title('ddof = ' + str(i), fontsize = 15)    ax[i].set_ylim(0, 0.07)    ax[i].set_xlim(0, 11)fig.suptitle(r'$s={\sqrt {{\frac {1}{10 - \mathrm{ddof}}}\sum _{i=1}^{10}\left(x_{i}-{\bar {x}}\right)^{2}}}$',             fontsize = 20, y=1.15);

Является ли этот факт важным и можно ли теперь быть уверенными в выводе, который мы сделали на основе Z-статистики? Первое, что приходит в голову - посмотреть, как может быть распределено стандартное отклонение большого количества выборок из генеральной совокупности. Давайте сделаем 5000 тысяч выборок по 10 человек, взятых из распределения N(80, 5) и посмотрим на распределение стандартного отклонения баллов этих выборок:

deviations  = np.std(stats.norm.rvs(80, 5, (5000, 10)), axis=1, ddof=1)sns.histplot(x=deviations ,stat='probability');

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

В конце концов, добавление нового сорта хмеля могло привести к большей специфичности вкуса, следовательно, вызвать большую разнонаправленность в предпочтениях: кому-то этот вкус пива стал нравиться больше, кому-то меньше, но средняя оценка могла остаться прежней. Например, если генеральная совокупность оценок "нового" пива теперь распределена как N(80, 10^{2}) , то Z-статистика и p-value для имеющихся 10-и оценок будет выглядеть не так обнадеживающе:

z = 10**0.5*(sample.mean()-80)/10p = 1 - (stats.norm.cdf(z) - stats.norm.cdf(-z))print(f'z = {z:.3}')print(f'p-value = {p:.4}')
z = 1.17p-value = 0.242

Оказывается, что для распределения N(80, 5^{2}) резульататы оказались значимыми, но как только мы оценили стандартное отклонение генеральной совокупности тем же значением, что и у выборки, т.е. представили, что генеральная совокупность распределена как N(80, 10^{2}) , то это очень сильно отразилось на значимости имеющихся оценок. Ведь вероятность получить такое же или даже еще более отклоненное среднее теперь равна не 2%, а почти 25%. Хотя, всего-то заменили \sigma на s .

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

Z = \frac{\bar{x} - \mu}{\frac{\sigma}{\sqrt{n}}} \; ; \;\;\; T = \frac{\bar{x} - \mu}{\frac{s}{\sqrt{n}}}.

По сути мы просто придумали новую T-статистику, которая отличается от Z-статистики только тем, что в знаменателе стоит не \sigma генеральной совокупности, а s выборки. А теперь сделаем 10000 выборок из N(80, {5}^2) , вычислим для каждой выборки Z- и T-статистику, а распределение значений изобразим в виде гистограмм:

fig, ax = plt.subplots(nrows=1, ncols=2, figsize = (12, 5))N = 10000samples = stats.norm.rvs(80, 5, (N, 10))statistics = [lambda x: 10**0.5*(np.mean(x, axis=1) - 80)/5,              lambda x: 10**0.5*(np.mean(x, axis=1) - 80)/np.std(x, axis=1, ddof=1)]title = 'ZT'bins = np.linspace(-6, 6, 80, endpoint=True)for i in range(2):    values = statistics[i](samples)    sns.histplot(x=values ,stat='probability', bins=bins, ax=ax[i])    p = values[(values > -2)&(values < 2)].size/N    ax[i].set_title('P(-2 < {} < 2) = {:.3}'.format(title[i], p))    ax[i].set_xlim(-6, 6)    ax[i].vlines([-2, 2], 0, 0.06, color='r');

Я все-таки не удержался и сделал гифку:

Код для гифки
import matplotlib.animation as animationfig, axes = plt.subplots(nrows=1, ncols=2, figsize = (18, 8))def animate(i):    for ax in axes:        ax.clear()    N = 10000    samples = stats.norm.rvs(80, 5, (N, 10))    statistics = [lambda x: 10**0.5*(np.mean(x, axis=1) - 80)/5,                  lambda x: 10**0.5*(np.mean(x, axis=1) - 80)/np.std(x, axis=1, ddof=1)]    title = 'ZT'    bins = np.linspace(-6, 6, 80, endpoint=True)    for j in range(2):        values = statistics[j](samples)        sns.histplot(x=values ,stat='probability', bins=bins, ax=axes[j])        p = values[(values > -2)&(values < 2)].size/N        axes[j].set_title(r'$P(-2\sigma < {} < 2\sigma) = {:.3}$'.format(title[j], p))        axes[j].set_xlim(-6, 6)        axes[j].set_ylim(0, 0.07)        axes[j].vlines([-2, 2], 0, 0.06, color='r')    return axesdist_animation = animation.FuncAnimation(fig,                                       animate,                                       frames=np.arange(7),                                      interval = 200,                                      repeat = False)dist_animation.save('statistics_dist.gif',                 writer='imagemagick',                  fps=1)

Просто мне кажется, что так можно продемонстрировать, как едва заметные черты могут оказаться чрезвычайно важными для исследователя. На что я хочу обратить ваше внимание? Во-первых, и слева, и справа мы видим свиду два одинаковых колоколообразных распределения. Вполне уместно предположить, что это два стандартных нормальных распределения, верно? Но, вот что любопытно, в N(0, 1) интервал [-2\sigma; 2\sigma] должен содержать около 95.5% всех значений. Для Z-статистики это требование выполняется, а для придуманной нами T-статистики не очень, потому что только 92-93% ее значений укладываются в заданный интервал. Казалось бы, что отличие не столь велико, чтобы заподозрить какую-то закономерность, но мы можем провести больше экспериментов:

statistics = [lambda x: 10**0.5*(np.mean(x, axis=1) - 80)/5,              lambda x: 10**0.5*(np.mean(x, axis=1) - 80)/np.std(x, axis=1, ddof=1)]quantity = 50N=10000result = []for i in range(quantity):    samples = stats.norm.rvs(80, 5, (N, 10))    Z = statistics[0](samples)    p_z = Z[(Z > -2)&((Z < 2))].size/N    T = statistics[1](samples)    p_t = T[(T > -2)&((T < 2))].size/N    result.append([p_z, p_t])result = np.array(result)fig, ax = plt.subplots()line1, line2 = ax.plot(np.arange(quantity), result)ax.legend([line1, line2],           [r'$P(-2\sigma < {} < 2\sigma)$'.format(i) for i in 'ZT'])ax.hlines(result.mean(axis=0), 0, 50, color='0.6');

В каждом из 50 экспериментов мы видим одно и тоже. Можно было бы заподозрить, что есть ошибка в коде, но его не так много, и легко убедиться в том, что вычисления выполняются верно. Так в чем же дело? А дело в том, что мы получили совершенно новый тип распределения! Снова взгляните на гифку с гистограммами распределений значений Z- и T- статистик, присмотритесь к основанию колокола каждой из них. Вам не кажется, что у распределения T-статистик основание чуть шире? Это хорошо видно по выпирающим за красные линии, так называемым - хвостам. А то, что эти хвосты несколько больше, или как еще говорят - тяжелее, чем у нормального распределения, так же будет означать, что мы будем наблюдать несколько больше сильных отклонений от вершины распределения. Проще говоря, мы теперь можем учитывать дисперсию выборки при оценке параметров генеральной совокупности. Однако, мы так и не ответили на вопрос - хорошо ли, можно ли, да и вообще зачем оценивать \sigma генеральной совокупности значением стандартного отклонения s выборки.

Задумаемся о дисперсии

Когда мы рассматривали распределение значений Z-статистик, мы обращали внимание только на распределение выборочного среднего \bar{x} и даже не задавались вопросом о том, что выборочное стандартное отклонение s тоже может как-то распределяться. Давайте возьмем 10000 выборок из N(80, 5) по 10 элементов в каждой и посмотрим, как будет выглядеть зависимость стандартного отклонения выборок от их среднего значения:

# если график строиться слишком долго,# то смените формат svg на png:#%config InlineBackend.figure_format = 'png'N = 10000samples = stats.norm.rvs(80, 5, (N, 10))means = samples.mean(axis=1)deviations = samples.std(ddof=1, axis=1)T = statistics[1](samples)P = (T > -2)&((T < 2))fig, ax = plt.subplots()ax.scatter(means[P], deviations[P], c='b', alpha=0.7,           label=r'$\left | T \right | < 2\sigma$')ax.scatter(means[~P], deviations[~P], c='r', alpha=0.7,           label=r'$\left | T \right | > 2\sigma$')mean_x = np.linspace(75, 85, 300)s = np.abs(10**0.5*(mean_x - 80)/2)ax.plot(mean_x, s, color='k',        label=r'$\frac{\sqrt{n}(\bar{x}-\mu)}{2}$')ax.legend(loc = 'upper right', fontsize = 15)ax.set_title('Зависимость выборочного стандартного отклонения\nот выборочного среднего',             fontsize=15)ax.set_xlabel(r'Среднее значение выборки ($\bar{x}$)',              fontsize=15)ax.set_ylabel(r'Стандартное отклонение выборки ($s$)',              fontsize=15);

Данный график интересен тем, что показывает зависимость среднего выборки от ее стандартного отклонения. На нем видно, что выборки, у которых \bar{x} и s отклоняется очень сильно и одновременно практически не встречаются, т.е. точки стараются избегать углов графика. Это значит, что, если мы берем выборку из нормального распределения N(\mu, \sigma^{2}) , то мы вряд ли увидим слишком большое значение \left | \bar{x} - \mu \right | при большом значении \left | s - \sigma \right | . А еще из данного графика следует, что экстремальные значения (красные точки) могут быть получены, только если выполняется следующее соотношение:

\left | \bar{x} - \mu \right | >\frac{2\sigma s}{\sqrt{n}}

Учитывая, что \sigma=1 , т.е. то, что мы снова измеряем расстояние от вершины распределения в сигмах, и то, что для нашего конкретного примера \mu = 80 , а n = 10 мы можем записать это условие следующим образом:

\left | \bar{x} - 80 \right | >\frac{2s}{\sqrt{10}}

Границы данного условия обозначены на графике черной линией, а поскольку мы раньше уже вычисляли доли значений которые попадали в интервал [-2\sigma; 2\sigma] , то мы можем сказать, что над черной линией находится около 92,5% всех точек.

Как интерпретировать это расположение и цвет точек? Давайте вспомним пример, от которого мы отталкивались. Мы сварили пиво с новым сортом хмеля, нашли десять человек (наша выборка) и попросили их оценить качество нового пива по 100-больной шкале. На основании их оценок мы пытаемся понять, как могут оценить качество пива все, кто его купит в дальнейшем (генеральная совокупность). Допустим среднее полученных 10-и оценок равно 82-м баллам, а стандартное отклонение равно всего 2-м баллам. Могут ли такие оценки быть получены случайно, если предположить, что генеральная совокупность имеет мат.ожидание \mu=80 , а стандартное отклонение такое же как у выборки, т.е. \sigma = s = 2 ? Чтобы выяснить это нужно вычислить Z-статистику:

Z = \frac{\bar{x} - \mu}{\frac{\sigma}{\sqrt{n}}} = \frac{82 - 80}{\frac{2}{\sqrt{10}}} \approx 3.16

Чтобы оценить вероятность случайного получения такого значения из N(80, 2^{2}) нужно вычислить p-value:

z = 10**0.5*(82-80)/2p = 1 - (stats.norm.cdf(z) - stats.norm.cdf(-z))print(f'p-value = {p:.2}')
p-value = 0.0016

Вероятность того что 10 человек случайно дадут средню оценку равную 82-м баллам очень мала и составляет менее 2%. Но данное утверждение является верным только если мы предполагаем что эти десять человек взяты из генеральной совокупности всех потребителей в которой их оценки распределены как N(80, 2^{2}) . И если мы действительно предполагаем, что \sigma = s = 2 , то мы можем быть вполне уверены, что оценки обусловлены новым сортом хмеля в составе пива, а не случайностью.

Такой результат может наблюдаться, например, в том случае, если мы пригласили экспертов в качестве респондентов. Только эксперты могли бы отметить небольшое улучшение качества пива (малое \left | \bar{x} - \mu \right | ) и при этом быть единодушны в оценках (малое s ).

С другой стороны мы могли бы пригласить крайне разношерстную публику в качестве 10 респондентов. Они могли бы оценить качество пива так же в 82 балла, но со стандартным отклонением равным, допустим, 9-и баллам. Может ли такая оценка быть получена случайно? Проверим:

z = 10**0.5*(82-80)/9p = 1 - (stats.norm.cdf(z) - stats.norm.cdf(-z))print(f'p-value = {p:.2}')
p-value = 0.48

Значит получить случайным образом 10 оценок с результатом \bar{x} = 82 и s = 9 из распределения N(80, 9^{2}) более чем вероятно. Следовательно, мы не можем утверждать, что новый сорт хмеля в составе пива как-то повлиял на изменение мнения генеральной совокупности о его качестве.

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

Код для гифки
import matplotlib.animation as animationfig, ax = plt.subplots(figsize = (15, 9))def animate(i):    ax.clear()    N = 10000        samples = stats.norm.rvs(80, 5, (N, i))    means = samples.mean(axis=1)    deviations = samples.std(ddof=1, axis=1)    T = i**0.5*(np.mean(samples, axis=1) - 80)/np.std(samples, axis=1, ddof=1)    P = (T > -2)&((T < 2))        prob = T[P].size/N    ax.set_title(r'зависимость $s$ от $\bar{x}$ при $n = $' + r'${}$'.format(i),                 fontsize = 20)    ax.scatter(means[P], deviations[P], c='b', alpha=0.7,               label=r'$\left | T \right | < 2\sigma$')    ax.scatter(means[~P], deviations[~P], c='r', alpha=0.7,               label=r'$\left | T \right | > 2\sigma$')    mean_x = np.linspace(75, 85, 300)    s = np.abs(i**0.5*(mean_x - 80)/2)    ax.plot(mean_x, s, color='k',            label=r'$\frac{\sqrt{n}(\bar{x}-\mu)}{2}$')    ax.legend(loc = 'upper right', fontsize = 15)    ax.set_xlim(70, 90)    ax.set_ylim(0, 10)    ax.set_xlabel(r'Среднее значение выборки ($\bar{x}$)',              fontsize='20')    ax.set_ylabel(r'Стандартное отклонение выборки ($s$)',              fontsize='20')    return axdist_animation = animation.FuncAnimation(fig,                                       animate,                                       frames=np.arange(5, 21),                                      interval = 200,                                      repeat = False)dist_animation.save('sigma_rel.gif',                 writer='imagemagick',                  fps=3)

С увеличением элементов в выборке мы наблюдаем, что выборочные \bar{x} и s стремятся к \mu и \sigma распределения N(\mu, \sigma^{2}) ,из которого эти выборки были взяты. Так что при больших значениях n мы можем спокойно пользоваться Z-статистиками, зная что при увеличении n стандартные отклонения выборки и генеральной совокупности будут все меньше отличаться друг от друга.

Но постойте! Разве мы не вернулись к тому, с чего начали? Ведь мы снова так и не ответили на вопрос - почему при малых размерах выборки мы предполагаем, что стандартное отклонение генеральной совокупности должно быть таким же, как у выборки. В нашем конкретном примере, с которого все началось, оценки 10-и респондентов выглядели так:

[89,99,93,84,79,61,82,81,87,82]

Среднее этой выборки \bar{x} = 83.7 , а стандартное оклонение s = 10.06 , причем мы приняли во внимание, что полномасштабный опрос, до введения в состав нового сорта хмеля, показал, что оценки распределены как N(80, 5^{2}) . Но вместо Z-статистики, мы решили использовать придуманную T-статистику, единственное отличие которой от Z-статистики состоит в том, что мы просто заменили в формуле \sigma генеральной совокупности на s выборки. По сути мы просто решили проверить гипотезу относительно генеральной совокупности, которая теперь по каким-то причинам распределена не как N(80, 5^{2}) , а как N(80, 10.06^{2}) - это вообще правильно? Разве не нужно было бы проверить это предположение для самых разных значений s ?: Почему мы не предполагаем, что N(80, 1^{2}) , N(80, 5^{2}) , N(80, 7^{2}) или любое другое значение \sigma ?

Чтобы ответить на этот вопрос, давайте снова вернемся к изначальному смыслу доказательства гипотез. Суть состоит в том, чтобы на основании выборки из генеральной совокупности сделать какие-то выводы о параметрах распределения этой генеральной совокупности. В нашем примере мы были уверены, что оценки его качества распределены как N(80, 5^{2}) , однако, когда мы сварили пиво с новым сортом хмеля, то мы получили 10 оценок, стандартное отклонение которых составило 10 баллов. В самом начале мы видели, что получить выборку с таким отклонением из N(80, 5^{2}) просто невозможно. А это значит, что мы просто вынуждены как-то по новому оценить стандартное отклонение генеральной совокупности, иначе все наши выводы не будут иметь никакого отношения к реальности.

Давайте сделаем вот что: среднее нашей выборки \bar{x} = 83.7 , а стандартное оклонение s = 10.06 , попробуем оценить вероятность получить такую выборку из N(80, \sigma^{2}) при разных значениях \sigma . Но поскольку вероятность получить конкретное значение из непрерывного распределения стремится к нулю, то мы оценим вероятность получения выборки, у которой 83<\bar{x}<84 и 9.5<s<10.5 :

N = 10000sigma = np.linspace(5, 20, 151)prob = []for i in sigma:    p = []    for j in range(10):        samples = stats.norm.rvs(80, i, (N, 10))        means = samples.mean(axis=1)        deviations = samples.std(ddof=1, axis=1)        p_m = means[(means >= 83) & (means <= 84)].size/N        p_d = deviations[(deviations >= 9.5) & (deviations <= 10.5)].size/N        p.append(p_m*p_d)    prob.append(sum(p)/len(p))prob = np.array(prob)fig, ax = plt.subplots()ax.plot(sigma, prob)ax.set_xlabel(r'Стандартное отклонение генеральной совокупности ($\sigma$)',              fontsize=20)ax.set_ylabel('Вероятность',              fontsize=20);

Как видите, максимум вероятности достигается при \sigma \approx 10 . Именно по этой причине, когда мы имеем дело с небольшими выборками, мы проверяем гипотезу о распределении генеральной совокупности не относительно известного нам значения ее стандартного отклонения \sigma , а относительно стандартного отклонения выборки s . Если параметры генеральной совокупности как-то изменились, то мы заранее предполагаем, что ее стандартное отклонение теперь равно стандартному отклонению выборки, просто потому, что вероятность такого изменения максимальна. Вот и все.

А это точно мы изобрели T-статистику?

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

Несомненно - абсурд! Однако, что-то подобное чуть не произошло с Уильямом Госсетом, когда он работал на пивоварне "Гиннес" и открыл t-распределение. Оценить практическую значимость некоторых статистических открытий очень сложно, но не в тех случаях, когда на кону астрономические потенциальные прибыли и стратегическое развитие компании, или даже государства. Например, последовательный критерий отношений правдоподобия Вальда, открытый в 1943 году во время войны вообще засекретили, так как этот критерий позволял на 50% уменьшить среднее число наблюдений в выборке. Короче, статистика - это крайне полезная наука.

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

t = \frac{\bar{x}-\mu}{\frac{s}{\sqrt{n}}}

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

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

t={\frac {Y_{0}}{{\sqrt {{\frac {1}{n}}\sum \limits _{{i=1}}^{n}Y_{i}^{2}}}}},

которая позволяет прояснить понятие "степени свободы" выборки. Пусть случайные переменные Y_{i} берутся из стандартного нормального распределения, т.е. {\displaystyle Y_{i}\sim {\mathcal {N}}(0,1),;i=0,\ldots ,n} , тогда значение n , т.е. объем выборки, как раз и определяет степень свободы распределения Стьюдента. А принадлежность какой-нибудь случайной величины к распределению Стьюдента с определенным числом степеней свободы обозначается, как:

t\sim {\mathrm {t}}(n)

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

Код для гифки
import matplotlib.animation as animationfig, ax = plt.subplots(figsize = (15, 9))def animate(i):    ax.clear()    N = 15000        x = np.linspace(-5, 5, 100)    ax.plot(x, stats.norm.pdf(x, 0, 1), color='r')        samples = stats.norm.rvs(0, 1, (N, i))        t = samples[:, 0]/np.sqrt(np.mean(samples[:, 1:]**2, axis=1))    t = t[(t>-5)&(t<5)]    sns.histplot(x=t, bins=np.linspace(-5, 5, 100), stat='density', ax=ax)        ax.set_title(r'Форма распределения $t(n)$ при n = ' + str(i), fontsize = 20)    ax.set_xlim(-5, 5)    ax.set_ylim(0, 0.5)    return axdist_animation = animation.FuncAnimation(fig,                                       animate,                                       frames=np.arange(2, 21),                                      interval = 200,                                      repeat = False)dist_animation.save('t_rel_of_df.gif',                 writer='imagemagick',                  fps=3)

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

{\displaystyle f_{t}(y)={\frac {\Gamma \left({\frac {n+1}{2}}\right)}{{\sqrt {n\pi }}\,\Gamma \left({\frac {n}{2}}\right)}}\,\left(1+{\frac {y^{2}}{n}}\right)^{-{\frac {n+1}{2}}}}

Но лучше воспользуемся готовой реализацией в SciPy:

Код для гифки
import matplotlib.animation as animationfig, ax = plt.subplots(figsize = (15, 9))def animate(i):    ax.clear()    N = 15000        x = np.linspace(-5, 5, 100)    ax.plot(x, stats.norm.pdf(x, 0, 1), color='r')    ax.plot(x, stats.t.pdf(x, df=i))        ax.set_title(r'Форма распределения $t(n)$ при n = ' + str(i), fontsize = 20)    ax.set_xlim(-5, 5)    ax.set_ylim(0, 0.45)    return axdist_animation = animation.FuncAnimation(fig,                                       animate,                                       frames=np.arange(2, 21),                                      interval = 200,                                      repeat = False)dist_animation.save('t_pdf_rel_of_df.gif',                 writer='imagemagick',                  fps=3)

Здесь гораздо лучше видно, что при увеличении n (параметр df в коде) распределение Стьюдента стремится к нормальному распределению. Что в общем-то и неудивительно, так как ранее мы уже видели, что при увеличении n выборки, ее стандартное отклонение все меньше колеблется вокруг стандартного отклонения генеральной совокупности.

t-распределение на практике

Давайте сразу выполним t-тест из модуля SciPy для выборки из нашего примера с пивом:

sample = np.array([89,99,93,84,79,61,82,81,87,82])stats.ttest_1samp(sample, 80)
Ttest_1sampResult(statistic=1.163532240174695, pvalue=0.2745321678073461)

А теперь попробуем получить тот же самый результат самостоятельно:

T = 9**0.5*(sample.mean() -80)/sample.std()T
1.163532240174695

Обратите внимание, что при вычислении мы использовали не размер выборки n = 10 , а количество степеней свободы df , которое вычисляется как n - 1 . Да, количество степеней свободы распределения Стьюдента на 1 меньше объема выборки, а связано это, как вы уже догадались, с разницей в вычислениях стандартного отклонения для выборки и генеральной совокупности. Добиться такого же результата мы можем и так:

T = 10**0.5*(sample.mean() -80)/sample.std(ddof=1)T
1.1635322401746953

Хорошо, значение t-статистики мы вычислили, как вычислить значение p-value? Думаю, вы и до этого уже догадались - нужно проделать все то же самое, что мы делали при вычислении p-value для Z-статистики, но только с использованием t-распределения:

t = stats.t(df=9)fig, ax = plt.subplots()x = np.linspace(t.ppf(0.001), t.ppf(0.999), 300)ax.plot(x, t.pdf(x))ax.hlines(0, x.min(), x.max(), lw=1, color='k')ax.vlines([-T, T], 0, 0.4, color='g', lw=2)x_le_T, x_ge_T = x[x<-T], x[x>T]ax.fill_between(x_le_T, t.pdf(x_le_T), np.zeros(len(x_le_T)), alpha=0.3, color='b')ax.fill_between(x_ge_T, t.pdf(x_ge_T), np.zeros(len(x_ge_T)), alpha=0.3, color='b')p = 1 - (t.cdf(T) - t.cdf(-T))ax.set_title(r'$P(\left | T \right | \geqslant {:.3}) = {:.3}$'.format(T, p));

Мы видим, что p-value чуть больше 27%, т.е. вероятность того, что результат получен случайно - довольно велика, особенно в соответствии с самым распространенным уровнем значимости \alpha = 0.05 , который меньше p-value более, чем в 5 раз. Что ж, наши результаты не являются значимыми, но они пригодны для демонстрации вычисления доверительного интервала - границы значений математического ожидания генеральной совокупности, в пределах которых оно находится с заданной вероятностью \alpha , обычно равной 0.95:

\textrm{CI}_{1 - \alpha} = \bar{x} \pm \left ( t_{\frac{\alpha}{2}, \textrm{df}} \right )\left ( \frac{s}{\sqrt{n}} \right )

Чтобы вычислить данный интервал в SciPy, достаточно воспользоваться методом interval с заданными параметрами loc (смещение) и scale (масштаб) для нашей выборки:

sample_loc = sample.mean()sample_scale = sample.std(ddof=1)/10**0.5ci = stats.t.interval(0.95, df=9, loc=sample_loc, scale=sample_scale)ci
(76.50640345566619, 90.89359654433382)

Учитывая, что \bar{x} = 83.7 , то мы можем утверждать, что с вероятностью \alpha = 0.95 математическое ожидание генеральной совокупности находится в пределах интервала [76.5; 90.9] . Чем шире доверительный интервал, тем больше вероятность того, что мы получим совсем другое выборочное среднее \bar{x} , если возьмем другую выборку такого же объема из генеральной совокупности.

Напоследок

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

Конечно, хотелось бы в конце вставить какую-нибудь гифку, но хочу закончить фразой Герберта Спенсера: "Величайшая цель образования - не знание, а действие", так что запускайте свои анаконды и действуйте! Особенно это касается самоучек, вроде меня.

Жму F5 и жду ваших комментариев!

Подробнее..

Ещё один модуль рисования графиков

25.05.2021 20:20:29 | Автор: admin
Лет пятнадцать назад потребовалось мне в программе для диплома отобразить график. Была бы программа на Builder или Delphi, всё было бы ничего, но только писал я для Windows на MFC, а там с классами графиков как-то не очень. И написал я тогда собственный модуль построения графиков. Три пятилетки прошло, а модуль остался, был переработан и я его иногда использую в своих поделках в QNX, Linux и Windows. Быть может, он пригодится чем-либо и вам.

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

Функции растеризации вынесены в подключаемые классы. Всего возможно на текущий момент три варианта: рисование штатными функциями Windows-GUI из MFC (класс CVideo_Windows), рисование штатными функциями Qt (класс CVideo_Qt) и программная растеризация (класс CVideo_Software с доработкой этот модуль можно использовать на микроконтроллерах). Перекодировку символов в требуемый для классов растеризации формат осуществляет класс CTranslator.

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

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

Выглядят нарисованные графики вот так:



Создание класса графика, например, для Windows в MFC выполняется следующим образом:

 CWnd_Graphics cWnd_Graphics;//класс графиков CGrData *cGrData_SinPtr;//указатель на данные графика синуса CGrData *cGrData_CosPtr;//указатель на данные графика косинуса//создаём графикCRect cRect_Client;((CStatic*)GetDlgItem(IDC_STATIC_MAIN_MAP))->GetClientRect(&cRect_Client);cWnd_Graphics.Create(WS_VISIBLE,cRect_Client,(CStatic*)GetDlgItem(IDC_STATIC_MAIN_MAP));  //настраиваем график cWnd_Graphics.GetIGraphicsPtr()->SetBackGroundColor(CGrColor(192,192,192)); cWnd_Graphics.GetIGraphicsPtr()->SetLegendBackGroundColor(CGrColor(230,230,230)); cWnd_Graphics.GetIGraphicsPtr()->SetAxisColor(CGrColor(0,0,0),CGrColor(0,0,0)); cWnd_Graphics.GetIGraphicsPtr()->SetTextValueColor(CGrColor(0,0,0),CGrColor(0,0,0),CGrColor(0,0,0)); cWnd_Graphics.GetIGraphicsPtr()->SetSelectRectangleColor(CGrColor(0,0,255));  cWnd_Graphics.GetIGraphicsPtr()->SetName("Графики функций"); //создаём графики cGrData_SinPtr=cWnd_Graphics.GetIGraphicsPtr()->AddNewGraphic(); cGrData_SinPtr->SetEnable(true); cGrData_SinPtr->SetGrColor(CGrColor(255,0,0)); cGrData_SinPtr->SetGrLineStyle(CGrLineStyle(IVideo::LINE_TYPE_SOLID,1,false,false)); cGrData_SinPtr->SetName("График синуса");  for(size_t n=0;n<1024;n++) {  double x=n;  double y=sin(x*0.01);  cGrData_SinPtr->AddPoint(x,y); } cGrData_CosPtr=cWnd_Graphics.GetIGraphicsPtr()->AddNewGraphic(); cGrData_CosPtr->SetEnable(true); cGrData_CosPtr->SetGrColor(CGrColor(0,0,255)); cGrData_CosPtr->SetGrLineStyle(CGrLineStyle(IVideo::LINE_TYPE_SOLID,3,false,false)); cGrData_CosPtr->SetName("График косинуса"); for(size_t n=0;n<1024;n++) {  double x=n;  double y=cos(x*0.01);  cGrData_CosPtr->AddPoint(x,y); } //приводим масштаб, чтобы отображались все графики CGrRect cGrRect; cWnd_Graphics.GetIGraphicsPtr()->FindViewRectangle(cGrRect); cWnd_Graphics.GetIGraphicsPtr()->SetRectangle(cGrRect); cWnd_Graphics.GetIGraphicsPtr()->OnMagnify(); cWnd_Graphics.GetIGraphicsPtr()->GetRectangle(cGrRect); cWnd_Graphics.GetIGraphicsPtr()->SetOriginalScale(cGrRect);

Здесь класс cWnd_Graphics обеспечивает связку класса графиков CGraphics с Windows, перенаправляя в класс CGraphics события, возникающие в Windows и обеспечивая отображение графика в событии перерисовки ON_WM_PAINT. Для других ОС эту связку потребуется переписать, с учётом используемой ОС. В данном примере через cWnd_Graphics.GetIGraphicsPtr() можно непосредственно обратиться к классу графиков CGraphics и настроить параметры отображения графиков, а так же попросить класс графиков создать новый график и вернуть на него указатель AddNewGraphic (будет получен указатель на класс CGrData). Удалять этот указатель самостоятельно нельзя график можно удалить только через функцию DeleteGraphic. В дальнейшем работа с графиком выполняется через полученный указатель.

Всего доступны следующие функции управления графиками:

CGrData* AddNewGraphic(void);//добавить график и получить указатель на созданный график  void DeleteAllGraphics(void);//удалить все графики из памяти  void DeleteGraphic(CGrData *cGrDataPtr);//удалить график из памяти  void FindRectangle(CGrRect &cGrRect) const;//найти описывающий прямоугольник всех активных графиков  void FindRectangleOfEndPoints(CGrRect &cGrRect,size_t points) const;//найти описывающий прямоугольник всех активных графиков за последние points точек  void FindRectangleOfEndTime(CGrRect &cGrRect,long double time) const;//найти описывающий прямоугольник всех активных графиков за последние time время  void FindViewRectangle(CGrRect &cGrRect) const;//найти описывающий прямоугольник всех активных графиков с запасом по высоте  void FindViewRectangleOfEndPoints(CGrRect &cGrRect,size_t points) const;//найти описывающий прямоугольник всех активных графиков с запасом по высоте за последние points точек  void FindViewRectangleOfEndTime(CGrRect &cGrRect,long double time) const;//найти описывающий прямоугольник всех активных графиков с запасом по высоте за последнее time время  void SetTimeDivider(double value);//установить временной делитель  double GetTimeDivider(void) const;//получить временной делитель  //функции отображения  void CancelSelect(void);//убрать выделение  void Redraw(void);//перерисовать изображение  void RedrawAll(void);//перерисовать всё  void OnOriginalScale(void);//перейти в режим оригинального масштаба  //функции цветовых настроек  void SetBackGroundColor(const CGrColor &cGrColor);//задать цвет фона  void SetLegendBackGroundColor(const CGrColor &cGrColor);//задать цвет фона легенды  void SetAxisColor(const CGrColor &cGrColor_AxisX,const CGrColor &cGrColor_AxisY);//задать цвет осей  void SetGridColor(const CGrColor &cGrColor_GridX,const CGrColor &cGrColor_GridY);//задать цвет сетки  void SetSelectRectangleColor(const CGrColor &cGrColor);//задать цвет прямоугольника выделения  void SetTextValueColor(const CGrColor &cGrColor_TextX,const CGrColor &cGrColor_TextY,const CGrColor &cGrColor_TextLegend);//задать цвет текста  //функции настройки стиля линий  void SetAxisLineStyle(const CGrLineStyle &cGrLineStyle_AxisX,const CGrLineStyle &cGrLineStyle_AxisY);//задать стиль осей  void SetGridLineStyle(const CGrLineStyle &cGrLineStyle_GridX,const CGrLineStyle &cGrLineStyle_GridY);//задать стиль сетки  //функции настройки системы координат  void SetRectangle(const CGrRect &cGrRect);//задать прямоугольник системы координат  void SetGridStep(long double step_x,long double step_y);//задать шаг сетки  void GetRectangle(CGrRect &cGrRect) const;//получить прямоугольник системы координат  void GetGridSize(long double &step_x,long double &step_y) const;//получить оптимальный шаг сетки  //функции дополнительных настроек  void SetEnableMagnify(bool enable);//разрешить увеличение  void SetEnableValue(bool x_value,bool y_value);//разрешить отображение значений  void SetOriginalScale(const CGrRect &cGrRect);//задать оригинальный масштаб  void SetMoveMode(bool inversion);//задать режим инверсии мышки  bool GetSelectedRectangle(CGrRect &cGrRect) const;//вернуть прямоугольник выделения  void GetClientRectangle(CGrRect &cGrRect) const;//вернуть прямоугольник графика  void SetName(const std::string &name);//установить название графика  bool GetUserMode(void) const;//получить режим управления пользователем  void SetUserMode(bool state);//установить режим управления пользователем

В принципе, настроить отображение можно достаточно гибко.

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

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

Ну вот, собственно, и всё, что можно сказать об этой примитивной поделке.
P.S. Мне тут подсказали, что я неудачно назвал функции. Извиняюсь, я немецкий язык учил и счёл, что по-английски будет вот так. :)
Подробнее..

Перевод Clustergram визуализация кластерного анализа на Python

28.05.2021 14:21:34 | Автор: admin

В этой статье, переводом которой мы решили поделиться специально к старту курса о Data Science, автор представляет новый пакет Python для генерации кластерограмм из решений кластеризации. Библиотека была разработана в рамках исследовательского проекта Urban Grammar и совместима со scikit-learn и библиотеками с поддержкой GPU, такими как cuML или cuDF в рамках RAPIDS.AI.


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

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

Маттиас Шонлау предложил другой подход кластерограмму. Кластерограмма это двухмерный график, отражающий потоки наблюдений между классами по мере добавления кластеров. Это говорит вам о том, как перетасовываются ваши данные и насколько хороши ваши сплиты. Тал Галили позже реализовал кластерограмму для k-средних в R. Я использовал реализацию Таля, перенёс ее на Python и создал clustergram пакет Python для создания кластерограмм.

clustergram в настоящее время поддерживает метод k-средних, использование scikit-learn (включая реализацию Mini-Batch) и RAPIDS.AI cuML (если у вас есть GPU с поддержкой CUDA), Gaussian Mixture Model (только scikit-learn) и иерархическую кластеризацию на основе scipy.hierarchy. В качестве альтернативы мы можем создать кластерограмму на основе меток и данных, полученных с помощью альтернативных пользовательских алгоритмов кластеризации. Пакет предоставляет API, подобный sklearn, и строит кластерные диаграммы с помощью matplotlib, что даёт ему широкий выбор вариантов оформления в соответствии со стилем вашей публикации.

Установка

Установить clustergram можно при помощи conda или pip:

conda install clustergram -c conda-forge

или

pip install clustergram

В любом случае вам нужно установить выбранный бэкенд (scikit-learn и scipy или cuML).

from clustergram import Clustergramimport urbangrammar_graphics as uggimport seaborn as snsimport matplotlib.pyplot as pltfrom sklearn.preprocessing import scalesns.set(style='whitegrid')

Давайте рассмотрим несколько примеров, чтобы понять, как выглядит кластерограмма и что с ней делать.

Набор данных о цветке ириса

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

iris = sns.load_dataset("iris")g = sns.pairplot(iris, hue="species", palette=ugg.COLORS[1:4])g.fig.suptitle("Iris flowers", y=1.01)

Похоже, что setosa относительно чётко определённая группа, тогда как разница между versicolor и virginica меньше, поскольку они частично перекрываются (или, в случае ширины чашелистика, полностью).

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

Давайте начнём с кластеризации методом k-средних. Чтобы получить стабильный результат, мы можем запустить кластерную программу с 1000 инициализаций.

data = scale(iris.drop(columns=['species']))cgram = Clustergram(range(1, 10), n_init=1000)cgram.fit(data)ax = cgram.plot(    figsize=(10, 8),    line_style=dict(color=ugg.COLORS[1]),    cluster_style={"color": ugg.COLORS[2]},)ax.yaxis.grid(False)sns.despine(offset=10)ax.set_title('K-Means (scikit-learn)')

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

Мы ищем разделение, т. е. отвечаем на вопрос, принёс ли дополнительный кластер какое-либо значимое разделение? Шаг от одного кластера к двум большой хорошее и чёткое разделение. От двух до трёх свидетельство довольно хорошего раскола в верхней ветви. Но с 3 по 4 видимой разницы нет, потому что новый четвёртый кластер почти не отличается от существующей нижней ветви. Хотя сейчас она разделена на две части, это разделение не даёт нам много информации. Таким образом, можно сделать вывод, что идеальное количество кластеров для данных Iris три.

Мы также можем проверить некоторую дополнительную информацию, например оценку силуэта или оценку Калинского Харабазса.

fig, axs = plt.subplots(2, figsize=(10, 10), sharex=True)cgram.silhouette_score().plot(    xlabel="Number of clusters (k)",    ylabel="Silhouette score",    color=ugg.COLORS[1],    ax=axs[0],)cgram.calinski_harabasz_score().plot(    xlabel="Number of clusters (k)",    ylabel="Calinski-Harabasz score",    color=ugg.COLORS[1],    ax=axs[1],)sns.despine(offset=10)

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

Набор данных о пингвинах со станции Палмера

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

penguins = sns.load_dataset("penguins")g = sns.pairplot(penguins, hue="species", palette=ugg.COLORS[3:])g.fig.suptitle("Palmer penguins", y=1.01)

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

data = scale(penguins.drop(columns=['species', 'island', 'sex']).dropna())cgram = Clustergram(range(1, 10), n_init=1000)cgram.fit(data)ax = cgram.plot(    figsize=(10, 8),    line_style=dict(color=ugg.COLORS[1]),    cluster_style={"color": ugg.COLORS[2]},)ax.yaxis.grid(False)sns.despine(offset=10)ax.set_title("K-Means (scikit-learn)")

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

Можно ли сказать, что их три? Поскольку мы знаем, что их должно быть три... Ну, не совсем. Разница между разделениями 23 и 34 незначительна. Однако здесь виновником является метод K ближайших соседей, а не кластерограмма. Он просто не может правильно кластеризовать эти данные из-за наложений и общей структуры. Давайте посмотрим, как работает смешанная Гауссова модель (Gaussian Mixture).

cgram = Clustergram(range(1, 10), n_init=100, method="gmm")cgram.fit(data)ax = cgram.plot(    figsize=(10, 8),    line_style=dict(color=ugg.COLORS[1]),    cluster_style={"color": ugg.COLORS[2]},)ax.yaxis.grid(False)sns.despine(offset=10)ax.set_title("Gaussian Mixture Model (scikit-learn)")

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

Подобная ситуация случается очень часто. Идеального случая не существует. В конечном счёте нам необходимо принять решение об оптимальном количестве кластеров. Clustergam даёт нам дополнительные сведения о том, что происходит между различными вариантами, как они расходятся. Можно сказать, что вариант с четырьмя кластерами в данных Iris не помогает. Также можно сказать, что пингвины Палмера могут быть сложными для кластеризации с помощью k-средних, что нет решающего правильного решения. Кластерограмма не даёт простого ответа, но она даёт нам лучшее понимание, и только от нас зависит, как мы её [кластерограмму] интерпретируем.

Установить clustergram можно с помощью conda install clustergram -c conda-forge или pip install clustergram. В любом случае вам всё равно придётся установить бэкенд кластеризации: либо scikit-learn, либо cuML. Документация доступна здесь, а исходный код здесь, он выпущен под лицензией MIT.

Если вы хотите поиграть с примерами из этой статьи, блокнот Jupyter находится на GitHub. Вы также можете запустить его в среде interactive binder в браузере. Более подробную информацию можно найти в блоге Тала Галили и оригинальных статьях Матиаса Шонлау.

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

Узнайте, как прокачаться и в других специальностях или освоить их с нуля:

Другие профессии и курсы
Подробнее..

Визуализация аналитики APIM Gravitee в Grafana

30.05.2021 10:22:52 | Автор: admin

Бесспорно, интерфейс Gravitee представляет достаточно наглядные и удобные средства визуализации работы шлюзов Gravitee. Но в любом случае, возникает потребность предоставить доступ к этим инструментам службе мониторинга, владельцам или потребителям API и при этом они могут находится вне закрытого контура, в котором расположен менеджер API. Да и иметь всю доступную информацию по различным API на одном экране всегда удобнее.
Видеть происходящее на шлюзах, при этом не вдаваясь в особенности пользовательского интерфейса Gravitee, а администраторам - не тратить время на создание пользователей и разделение ролей и привилегий внутри Gravitee.
На Хабре уже была пара статей посвященных APIM Gravitee, тут и тут. По этому, в своей заметке, буду подразумевать, что читатель уже знаком с процессом установки/настройки APIM Gravitee и Grafana, рассмотрю только процесс настройки их интеграции.

Почему нельзя пойти простым путём?

По умолчанию, хранилищем для аналитики Gravitee является ElasticSearch. Информация накапливается в четырёх различных индексах, с посуточной разбивкой:

  • gravitee-request-YYYY.MM.DD - здесь хранится информация по каждому запросу (аналог access.log в nginx). Это наша основная цель;

  • gravitee-log-YYYY.MM.DD - здесь уже хранится более подробная информация о запросе (при условии, что включена отладка, см. рисунок ниже). А именно полные заголовки запросов и ответов, а также полезная нагрузка. В зависимости от настроек, логироваться может как обмен между потребителем и шлюзом, так и/или шлюзом и поставщиком API;

    Экран включения/отключения расширенного логированияЭкран включения/отключения расширенного логирования
  • gravitee-monitor-YYYY.MM.DD - этот нас не интересует;

  • gravitee-health-YYYY.MM.DD - этот нас не интересует.

И казалось бы, что может быть проще: подключай ElasticSearch в качестве источника данных в Grafana и визуализируй, но не всё так просто.
Во первых, в индексе хранятся только идентификаторы объектов, т.е. человеко-читаемых имён поставщиков и потребителей, вы там не увидите. Во вторых, получить полную информацию соединив данные из двух источников непосредственно в интерфейсе Grafana, крайне проблематично. Gravitee хранит информацию о настройках и статистику своей работы в разных местах. Настройки, в MongoDB или PostgreSQL, по сути статическая информация. Таким образом в одном месте у нас (в терминах Grafana) - таблица, в другом - временной ряд.

B как же быть?

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

Схема взаимодействия модулей GraviteeСхема взаимодействия модулей Gravitee

Ну что же, за дело!

Все ниже описанные действия актуальны для следующей конфигурации: CentOS 7, APIM Gravitee 3.6, СУБД PostgreSQL 11, ElasticSearch 7.+

Начнём с интеграции PostgreSQL и ElasticSearch. Сам процесс интеграции достаточно прост и делится на следующие шаги:

  1. Устанавливаем расширение multicorn11 и если не установлен pip, то ставим и его:

    yum install multicorn11 python3-pip
    
  2. Далее из pip-репозитория, устанавливаем библиотеку python3 для работы с ElasticSearch:

    pip3 install pg_es_fdw
    
  3. Далее, переходим к настройке PostgreSQL. Подключаемся целевой БД и добавляем расширение multicorn и подключаем необходимую библиотеку:

    GRANT USAGE on FOREIGN DATA WRAPPER multicorn TO gatewaytest;GRANT USAGE ON FOREIGN SERVER multicorn_es TO gatewaytest;
    
     CREATE EXTENSION multicorn; CREATE SERVER multicorn_es FOREIGN DATA WRAPPER multicorn  OPTIONS (wrapper 'pg_es_fdw.ElasticsearchFDW');
    
  4. Выдаём права, непривилегированному пользователю. В нашем случае это logreader:

    GRANT USAGE on FOREIGN DATA WRAPPER multicorn TO logreader;GRANT USAGE ON FOREIGN SERVER multicorn_es TO logreader;
    
  5. Для удобства, создадим отдельную схему logging, владельцем которой будет наш пользователь logreader:

    CREATE SCHEMA logging AUTHORIZATION logreader;
    
  6. Создадим родительскую таблицу, к которой мы будем подключать новые индексы и удалять не актуальные:

    CREATE TABLE logging.requests (  id varchar(36),  "@timestamp" timestamp with time zone,  api varchar(36),  "api-response-time" int,  application varchar(36),  custom json,  endpoint text,  gateway varchar(36),  "local-address" varchar(16),  method int,  path text,  plan varchar(36),  "proxy-latency" int,  "remote-address" varchar(16),  "request-content-length" int,  "response-content-length" int,  "response-time" int,  sort text,  status int,  subscription varchar(36),  uri text,  query TEXT,  score NUMERIC) PARTITION BY RANGE("@timestamp");
    

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

  7. Для подключения и отключения индексов, создадим небольшой shell-скрипт и будем запускать его раз в сутки через cron:

    #!/bin/shNEWPART=${1:-$(date +'%Y.%m.%d')}OLDPART=$(date --date='14 days ago' +'%Y.%m.%d')curl http://gateway.corp/testpsql gateway -U logreader -c "CREATE FOREIGN TABLE logging.\"requests_${NEWPART}\"PARTITION OF logging.requests   FOR VALUES FROM ('${NEWPART} 00:00:00') TO ('${NEWPART} 23:59:59')SERVER multicorn_esOPTIONS (host 'els-host',  port '9200',  index 'gravitee-request-${NEWPART}',  rowid_column 'id',  query_column 'query',  query_dsl 'false',    score_column 'score',  sort_column 'sort',  refresh 'false',  complete_returning 'false',  timeout '20',  username 'elastic-ro',  password 'Sup3rS3cr3tP@ssw0rd');"    psql gateway -U gatewaydev -c "drop foreign table logging.\"requests_${OLDPART}\""
    

    Немного пояснений:

    • NEWPART - текущая дата, для формирования имени партиции , при подключении нового индекса из ElasticSearch;

    • OLDPART - дата истекшего, неактуально индекса, здесь это 14 дней (определяется исходя из настроек ES Curator). Удалять партиции, ссылающиеся на несуществующие - обязательно. В противном случае запросы, к родительской таблице, будут прерываться с ошибками;

    • Вызов 'curl http://gateway.corp/test', необходим для того, что бы создавался индекс текущего дня, так как он создаётся в момент первого обращения к любому поставщику API. Если его не создать, то это будет приводить к ошибке, описанной выше. Такая проблема больше актуальна для тестовых стендов и стендов разработки;

    • Затем, создаём партицию на индекс текущего дня;

    • И на последнем шаге - удаляем неактуальный индекс.

    • Проверяем что всё работает

      TABLE logging.requests LIMIT 1;
      

      Если всё правильно, то должны получить похожий результат

    -[ RECORD 1 ]-----------+-------------------------------------id                      | 55efea8a-9c91-4a61-afea-8a9c917a6133@timestamp              | 2021-05-16 00:00:02.025+03api                     | 9db39338-1019-453c-b393-381019f53c72api-response-time       | 0application             | 1custom                  | {}endpoint                | gateway                 | 7804bc6c-2b72-497f-84bc-6c2b72897fa9local-address           | 10.15.79.29method                  | 3path                    | plan                    | proxy-latency           | 2remote-address          | 10.15.79.27request-content-length  | 0response-content-length | 49response-time           | 2sort                    | status                  | 401subscription            | uri                     | /testquery                   | score                   | 1.0
    

Рисуем графики

И вот мы подошли к тому, ради чего всё и делалось - визуализируем статистику Gravitee. Благодаря тому, что для доступа к аналитике используется единая точка входа, а именно СУБД PostgreSQL, это даёт дополнительные возможности. Например, выводить статическую информацию: количество поставщиков, количество потребителей и их статусы; количество и состояние подписок; параметры конфигурации для поставщика и многое другое, наряду с динамическими данными.
В том числе хотелось бы отметить, что у поставщиков и потребителей имеется раздел Metadata, которые можно заполнять кастомными данными и так же выводить в дашборды Grafana.

Вот тут:

Раздел Metadata в GraviteeРаздел Metadata в Gravitee

А вот так это можно отобразить в Grafana:

Вариант отображения Metadata в GrafanaВариант отображения Metadata в Grafana
SELECT  name "Наименование",  value "Значение"FROM  metadataWHERE  reference_id='${apis}'

Пример комплексного экрана

Вариант комплексного экранаВариант комплексного экрана

APIs (статика) - общее количество поставщиков и количество активных.

SELECT COUNT(*) AS "Всего" FROM apis;SELECT COUNT(*) AS "Активных" FROM apis WHERE lifecycle_state='STARTED';

Для Applications, запросы составляются по аналогично, только из таблицы applications

API Hits - количество вызовов по каждому поставщику. Тут уже немного по сложнее

SELECT  date_trunc('minute',"@timestamp") AS time,  apis.name,ee с Grafana  COUNT(*)FROM  logging.requests alJOIN  apis ON al.api = apis.idWHERE  query='@timestamp:[$__from TO $__to]'GROUP BY 1,2

Average response time by API - среднее время ответа, по каждому поставщику считается аналогичным способом.

SELECT  date_trunc('minute',"@timestamp") AS time,  apis.name,  AVG(al."api-response-time")FROM  logging.requests alJOIN  apis ON al.api = apis.idWHERE  query='@timestamp:[$__from TO $__to]'GROUP BY 1,2

Еще один интересный показатель Hits, by gateways, это равномерность распределения запросов по шлюзам. Считается так:

SELECT  date_trunc('minute',"@timestamp") as time,  al."local-address",  COUNT(*)FROM  logging.requests alWHERE  query='@timestamp:[$__from TO $__to]'GROUP BY 1,2
График распределения запросов по шлюзамГрафик распределения запросов по шлюзам

Заключение

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

P.S.

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

Конструктивная критика, пожелания и предложения приветствуются!

Подробнее..

Топ-5 софт-навыков дизайнера в банке

04.06.2021 14:18:19 | Автор: admin

Соавтор:Кузнецова Юлия Андреевна - UX-писатель Экосистемы РСХБ

Каким должен быть дизайнер в банке, чтобы и продукт хороший создавал, и коллеги не жаловались. Смотрим через призму софт-навыков вместе с UX-дизайнерами РСХБ.

#1 Коммуникабельность: не просто коллега человек

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

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

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

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

#2 Презентация: готов объяснить свою работу

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

Крутой дизайнер умеет переключить фокус с творческого и ранимого художника на хладнокровного адвоката. Если кто-то говорит: Это плохо сделано, это не значит, что ты плохой. Это значит, что коллега видит в данном решении слабые стороны. Продвинутый дизайнер отключает эмоции и включает голову. Объясняет и обосновывает свою позицию. Задаёт вопросы, внимательно слушает, делает все, чтобы добраться до истины и хорошего продукта.

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

#3 Аналитика: думает не только о красоте

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

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

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

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

#4 Понимание бизнес-процесса: переводит с банковского языка на человеческий

Юный дизайнер: мастерски работает в Фигме, любит компоненты, боготворит филигранную верстку.

Продвинутый дизайнер: понимает, как устроен бизнес и финансы.

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

Знает о банке

Знает о пользователе

как устроены процессы внутри банка

над чем работает банк

какие боли клиентов решает

как себя позиционирует и каких принципов придерживается

как работают конкуренты

кто и зачем пользуется банковскими продуктами

с какими проблемами сталкивается

что делает, чтобы решить эти проблемы

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

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

#5 Любознательность: постоянно учится новому

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

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

Подробнее..

Телеграмм-бот для анализа опционов

07.06.2021 10:04:02 | Автор: admin

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

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

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

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

Запоминание состояния бота между вебхуками

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

$id_init = file_get_contents('php://input');$id=sbs($id_init, '"from":{"id":',',"is_bot":'); //в эту переменную записываем уникальный номер пользователяfunction sbs ($str,$m1,$m2){ //из строки str возвращает подстроку между двумя метками-словами m1 и m2$p1=strpos($str,$m1)+strlen($m1); //длина слова-метки слева$p2=strpos($str,$m2);return substr($str,$p1,$p2-$p1);}

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

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

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

  3. файлы для построения графика

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

Построение графика для анализа портфеля

График строится по состоянию на момент экспирации опционов, что дает возможность упрощенно представить график каждого инструмента в виде y=kx+b,

где y это размер прибыли/убытка на момент экспирации

х стоимость базового актива

После записи данных в портфель создаются файлы с расширением png с помощью библиотеки GD в несколько этапов:

1) определение Х-координат всех точек перелома на графике (очевидно, что сумма всех линий даст ломанную кривую)

2) определение масштаба изображения (средняя всех координат точек перелома по оси Х, максимальное отклонение от среднего, и максимальный размер по оси Y)

3) создание ассоциативного массива точек, в котором координата X это ключ, координата Y величина, для всей цифровой полуплоскости:

$typ опцион колл, пут или фьючерс

$q количество ( отрицательное продажа)

$cena цена приобретения ценной бумаги

$strike страйк для опционов

$x0 начальная координата по оси Х

$sx масштаб по оси Х

function pparr($typ, $q, $cena, $strike,$x0,$sx){ //функция выдает одномерный массив - координаты x=>y точек по //типу цб, направлению (покупка продажа), цене приобретения и страйку (для опционов)  if ($q<0) { $q=-$q; $drct='-'; }  else $drct='+'; $a=array(); $b=array(); $delta=$sx; //расстояние между точками равно масштабу $scalx for ($i=0;$i<740;$i++){ //кол во точек 740 определено заранее $xkk=$x0+$delta*$i; //значение по оси X if ($typ=='fut') { if ($drct=='+') $a[$xkk]=($xkk-$cena)*$q; else $a[$xkk]=(-$xkk+$cena)*$q; } if ($typ=='call'){ if ($drct=='+') { if ($xkk<=$strike) $a[$xkk]=-$cena*$q; else $a[$xkk]=$q*($xkk-$strike-$cena);} else { if ($xkk<=$strike) $a[$xkk]=$q*$cena; else $a[$xkk]=(-$xkk+$strike+$cena)*$q;}  } if ($typ=='put'){ if ($drct=='+') { if ($xkk<=$strike) $a[$xkk]=(-$xkk+$strike-$cena)*$q; else $a[$xkk]=-$q*$cena;} else { if ($xkk<=$strike) $a[$xkk]=($xkk-$strike+$cena)*$q; else $a[$xkk]=$cena*$q;}  } $b[(string)$xkk]=(string)$a[$xkk]; }return $b; };

4) создание файлов изображений, при этом одновременно строится график для каждой строки из портфеля ( зеленый цвет) и результирующий для портфеля ( красный цвет). Кроме того, know-how заключается в том, что одновременно строится еще четыре изображения для увеличения/уменьшения изображения по оси Х и по оси Y. За счет этого достигается эффект работы он-лайн с клавишами X+,X-,Y+,Y- под графиком. Таким образом, для каждого пользователя в каждый момент времени существует пять файлов изображения.

Анализ рисков опционного портфеля

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

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

Гарантийное обеспечение для фьючерсов считается по формуле, которая предложена самой Московской биржей:

ГО=БГО+(Цена-Расчетная_Цена)*БП;

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

БП стоимость пункта в рублях, зависит, в частности, от курса доллара,

БГО базовое гарантийное обеспечение, определяется биржей, может произвольно увеличиваться, например, перед праздниками

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

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

Вот как выглядит покупка колла в нашем боте:

Запись в портфеле покупки опциона "колл" страйк 75000 дата экспирации 03.06.2021 по цене 25 Запись в портфеле покупки опциона "колл" страйк 75000 дата экспирации 03.06.2021 по цене 25 График зависимости прибыли/убытка по купленному опциону "колл" в зависимости от стоимости базового актива на дату экспирации График зависимости прибыли/убытка по купленному опциону "колл" в зависимости от стоимости базового актива на дату экспирации

Гарантийное обеспечение=23.

Что нам показывает график: если стоимость базового актива ( в данном случае стоимость фьючерса на курс рубля к доллару) [вечером] 03.06.2021 будет 75000 и ниже, то наш убыток составит 23 . При повышении этой стоимости до 75023 мы выйдем в безубыток, при дальнейшем росте получим прибыль.

Что мы имеем с точки зрения риска: не при каких обстоятельствах наш убыток не превысит сумму 23. Следовательно, наш опцион совершенно не похож на фьючерс, и в расчете ГО мы можем записать просто сумму 23.

Покупка пута примерно та же картина.

Продажа пута.

Запись в портфеле продажи опциона "пут" страйк 72750 по цене 44 с датой экспирации 03.06.2021Запись в портфеле продажи опциона "пут" страйк 72750 по цене 44 с датой экспирации 03.06.2021График зависимости прибыли/убытка по проданному опциону "пут" от стоимости базового актива на дату экспирацииГрафик зависимости прибыли/убытка по проданному опциону "пут" от стоимости базового актива на дату экспирации

Гарантийное обеспечение= 5436.

Можно убедиться, что при стоимости базового актива выше 72750 мы имеем прибыль 44. При снижении стоимости БА до 72706 мы выходим в ноль. При дальнейшем падении стоимости БА наш убыток НИЧЕМ НЕ ОГРАНИЧИВАЕТСЯ.

С точки зрения рисков это фьючерс, купленный по цене 72706. Подставляем это число в формулу ГО для фьючерса и получаем ГО для опциона! Этот ГО достаточно велик (5436), но может превратится в прибыль в течение нескольких дней.

Стоит ли овчинка выделки? Если при условных затратах сегодня 5436 мы получим через несколько дней 44? Мне кажется, риск очень велик. С другой стороны, курс рубля к доллару был таким примерно 3 месяца назад.

С продажей колла будет аналогичная ситуация.

А если одновременно продать пут и кол?

Запись в портфеле продажи опциона "пут" страйк 72750 и продажи опциона "колл" страйк 75000Запись в портфеле продажи опциона "пут" страйк 72750 и продажи опциона "колл" страйк 75000График зависимости прибыли/убытка по портфелю от стоимости базового актива на дату экспирации График зависимости прибыли/убытка по портфелю от стоимости базового актива на дату экспирации

Гарантийное обеспечение не изменилось!

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

с наибольшим ГО , оно и будет мерилом риска.

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

Интерфейс телеграмм-бота

На экране отображаются следующие группы данных:

  • подробная инструкция по работе с ботом

  • текущие котировки ближайших фьючерсов по трем базовым активам

  • таблица Портфель, отображающая состав портфеля, сделанная с помощью интерфейса телеграмм ( можно редактировать)

  • значение гарантийного обеспечения

  • таблица, дублирующая состав портфеля, построенная как изображение формата png, которую можно копировать и сохранять

  • нижняя клавиатура, которую можно скрывать, с кнопками: Добавить позицию в портфель, Анализ портфеля с помощью графика, Обновить котировки

Таблица Портфель создана с помощью InlineKeyboard.

При нажатии на клавиши в этой таблице происходят следующие действия:

  • клавиша в столбце К-во редактирует количество ценных бумаг, соответствующее строке, в которой нажата клавиша

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

  • клавишей Х можно удалить строку.

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

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

Заключение

В телеграмм-боте реализовано запоминание портфеля ценных бумаг для каждого пользователя. Под ценными бумагами понимаются опционы и фьючерсы, базовым активом для которых являются: курс рубля к доллару (Si), стоимость нефти брент (BR), а также индекс РТС (RI). Это самые высоколиквидные деривативы на московской бирже.

В телеграмм-бот заложен алгоритм подсчета гарантийного обеспечения, которое является мерой риска.

С помощью телеграмм-бота можно анализировать опционный портфель на графике прибылей/убытков (P/L график).

Протестировать телеграмм-бота можно по ссылке t.me/@test09062020bot. Или попробовать найти в телеграмме по названию опционный портфель.

Подробнее..

Анализируем слона вместе с коллегами

15.06.2021 14:10:21 | Автор: admin

Если ваша жизнь DBA, сопровождающего PostgreSQL, наполнена вопросами "а почему так медленно?" и "как сделать, чтобы запрос не тормозил?", наш сервис анализа и визуализации планов запросовexplain.tensor.ru сделает ее немного легче за счет привлечения коллег и обновленных подсказок.

м/ф "Следствие ведут Колобки"м/ф "Следствие ведут Колобки"

"Ландон из зе кепитал оф Грейт Британ"

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

Обсуждайте проблемный план там, где вам удобноОбсуждайте проблемный план там, где вам удобно

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

Подсказки к плану

Про базовый набор подсказок и способов ими воспользоваться я уже рассказывал в статье "Рецепты для хворающих SQL-запросов" - теперь мы сделали их еще больше и удобнее!

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

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

  • несколько подсказок у одного узла

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

  • на самом видном месте

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

Все подсказки - вместе, клик - и вы на местеВсе подсказки - вместе, клик - и вы на месте

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

Масштабируемая диаграмма

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

Карта будущего сражения за производительностьКарта будущего сражения за производительность

Пользуйтесь! Возникнут идеи или замечания - прошу в комментарии.

Подробнее..

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

17.06.2021 18:06:38 | Автор: admin

В школе все мы решали задачки вида едет из пункта А в пункт Б. Речь преимущественно шла о скорости и времени как быстро доберётся транспортное средство? Реальность, однако, подбрасывает задачки значительно интереснее: Существует масштабная ритейл-сеть по продаже товаров, которой необходимо, чтобы огромное количество номенклатурных позиций доезжало в каждый из 17000 магазинов, расположенных на половине площади самой большой страны в мире, вовремя и в нужном количестве. Для решения такой задачи в X5 Group существует ряд реализованных решений, и одним из самых важных является процесс автозаказа товаров.

Техническую поддержку этого направления в X5 Group обеспечивает команда 2-SAP Логистики.

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

Автозаказ это комплекс процессов управления запасами и заказами

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

Планирование заказа.

Формирование заказа.

Отправка заказа.

Экономическое обоснование поддерживаемого уровня.

Контроль за состоянием запасов.

Управление ассортиментом.

Для этого используются разные механизмы прогнозирования в зависимости от формата магазина и гибких настроек в системе.

Прогноз продаж, построенный SAP это прогноз спроса на период 42 дня, который строится на исторических данных по продажам. Период анализа продаж и модель прогноза определяется из настройки Профиля прогноза.

В периоде анализа продажи могут быть Плановые регулярные продажи и Внеплановые промо продажи.

Прогноз запускается ночью (по календарю) в фоновом режиме массово для всех комбинаций товар-магазин, в рамках ассортимента и является частью Автозаказа.

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

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

Около 04:30 эти данные по потребности из JDA поступают в SAP, где обрабатываются фоновым заданием с интервалом запуска каждые 15 минут, в результате чего создаются Автозаказы.

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

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

Через автозаказ пополняется до 80% основного ассортимента магазинов, а ручное пополнение выполняется для товаров in-out.

Схема реализации товара:

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

Так выглядит децентрализованная цепочка:

- в 00:00 стартует задание по АЗ, в период его работы по нашему товару рассчиталась потребность, сформировался автозаказ, и он бы отправлен из SAP в программу магазина GK(Пятёрочка), там у сотрудников магазина есть время отредактировать рассчитанное кол-во (в большую или меньшую сторону) до наступления времени автосогласования (например, до 11 утра, время зависит от тайминга поставщика, т.е. до какого момента он должен получить заказ, чтобы собрать и вовремя привезти заказ в магазин);

- после внесения корректировок эти изменения отправляются в SAP, и дальше происходит процесс автосогласования заказ согласовывается и автоматически отправляется поставщику по почте или его провайдеру (в зависимости от того, по какой схеме поставщик работает). Например, если поставщик работает по EDI (через электронный документооборот), то после доставки заказа в его систему, нам провайдер передает информацию о размещении заказа в своем ПО и в системе ER2 заказ переходит в статус Размещен;

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

Автозаказ (Ночной профиль) фоновое задание, имеющее внутри себя несколько шагов, выполняющих определённые функции. Результатом задания является создание заказов на закупку товаров. Запускается параллельно по магазинам.

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

  • Продажи за период в прошлом

  • Прогноз на период в будущее (на основании продаж строится прогноз SAP)

  • Остатки на объекте получателе (остаток в SAP на момент расчёта)

  • Открытые заказы (заказы на поставку и возвраты)

Продажи приходят из магазинов каждый вечер до расчёта фонового задания Автозаказа.

Факт получения продаж запускает расчет Автозаказа. Если до 03:00 продажи не получены SAP ERP, то происходит безусловный запуск задания Автозаказа.

Автозаказ формируется строго по графику заказа/поставки согласованному с поставщиком, так категории FRESH и ULTRA FRESH, как правило, заказываются и поставляются в магазины ежедневно.

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

Для того, чтобы заказать оптимальное кол-во товара, учитывается прогноз, который строится в SAP в момент расчета Автозаказа.

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

Автозаказ определяет цепочку и поставщика по каждому товару для каждого магазина, согласно записей книги источников поставок (смежный функционал, автоматизирован, и ведется также в SAP) и формирует заказ.

Запускается задание ежедневно в 00:00. Первыми выполняются расчеты по регионам Сибирь и Урал, далее по регионам Волга, ЮГ, Москва, Центральные регионы и Северо-запад. Последовательность выполнения крайне важна, т.к. наши магазины расположены в разных часовых поясах, и пока в Москве все сладко спят, в Сибири уже в разгаре рабочий день.

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

Для этого мы реализовали различные системы мониторинга, инструменты для анализа, а в случае необходимости и возможность отправки в магазины резервных шаблонов заказа.<o:p>

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

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

Подробнее..

Категории

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

  • Имя: Murshin
    13.06.2024 | 14:01
    Нейросеть-это мозг вселенной.Если к ней подключиться,то можно получить все знания,накопленные Вселенной,но этому препятствуют аннуннаки.Аннуннаки нас от неё отгородили,установив в головах барьер. Подр Подробнее..
  • Имя: Макс
    24.08.2022 | 11:28
    Я разраб в IT компании, работаю на арбитражную команду. Мы работаем с приламы и сайтами, при работе замечаются постоянные баны и лаги. Пацаны посоветовали сервис по анализу исходного кода,https://app Подробнее..
  • Имя: 9055410337
    20.08.2022 | 17:41
    поможем пишите в телеграм Подробнее..
  • Имя: sabbat
    17.08.2022 | 20:42
    Охренеть.. это просто шикарная статья, феноменально круто. Большое спасибо за разбор! Надеюсь как-нибудь с тобой связаться для обсуждений чего-либо) Подробнее..
  • Имя: Мария
    09.08.2022 | 14:44
    Добрый день. Если обладаете такой информацией, то подскажите, пожалуйста, где можно найти много-много материала по Yggdrasil и его уязвимостях для написания диплома? Благодарю. Подробнее..
© 2006-2024, personeltest.ru