В этой статье, переводом которой мы решили поделиться специально к старту курса о 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, где через подкреплённую теорией практику студенты учатся ориентироваться в данных, делать обоснованные выводы и действовать с открытыми глазами.

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

• Профессия Data Analyst

• Курс по Data Engineering

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

В статье с единой точки зрения обсуждаются три часто встречающихся на практике одновыборочных теста: тест знаков, t-тест и тест Уилкоксона (Signed-Rank Wilcoxon test) непараметрической процедуры, мощность которой сравнима с мощностью t-теста в случае нормально распределенной выборки, и превышает мощность t-теста в случае, если распределение выборки имеет более тяжелые хвосты по сравнению с нормальным распределением.

1. Определим модель для параметра положения (location model) следующим образом. Пусть обозначает случайную выборку, полученную по следующему закону

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

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

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

3.1. Классический тест знаков (sign test) основан на статистике

где для соответственно. Пусть

Тогда . Здесь предполагается, что ни одно из значений не равно нулю (на практике, равные нулю значения из выборки исключают, а объем выборки корректируют). При условии , статистика имеет биномиальное распределение с числом испытаний и вероятностью успеха . Пусть наблюдаемая величина тогда p-value для теста знаков равно , где функция биномиального распределения с параметрами и (R функция pbinom возвращает значения cdf для биномиального распределения).

Заметим, что в тесте знаков распределение статистики при нулевой гипотезе не зависит (свободно) от вида распределения .

3.2. Следующий традиционный t-тест (t-test) основан на сумме наблюдений. По аналогии можно записать

Заметим, что распределение статистики зависит от плотности распределения . Обычно t-тест записывают в форме t-отношения

где и соответственно, выборочное среднее и стандартное отклонение. Если выборка получена из нормального распределения, то статистика имеет t-распределение Стьюдента с степенью свободы. Пусть наблюдаемое по выборке значение . Тогда p-value для t-теста равно , где функция t-распределения Стьюдента c степенью свободы (R функция pt возвращает значения cdf для t-распределения). Это точное значение p-value в случае нормального распределения, в противном случае это аппроксимация.

3.3. Отличие t-теста от теста знаков состоит в том, что статистика t-теста является функцией расстояний элементов выборки относительно нуля в дополнение к их знакам.

Выбранная нами статистика теста Уилкоксона (signed-rank Wilcoxon test) хороша тем, что использует лишь ранги этих расстояний. Обозначим ранг среди всех , упорядоченных от меньшего значения к большему. Тогда статистика Уилкоксона имеет вид

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

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

Тогда p-value для теста Уилкоксона равно , где функция распределения статистики Уилкоксона для выборки размера (R функция psignrank возвращает значения cdf распределения ).

4. Техника построения доверительных интервалов широко используется при решении практических задач. Каждый из рассмотренных выше тестов: тест знаков, t-тест и тест Уилкоксона имеет соответствующую оценку и доверительный интервал для параметра положения . Рассмотрим далее имеющиеся результаты.

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

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

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

4.3. Оценкой параметра положения , связанной с тестом Уилкоксона является оценка Ходжеса-Лемана (Hodges-Lehmann)

Парные средние , называются средними Уолша (Walsh averages) выборки. Пусть упорядоченный набор средних Уолша. Тогда доверительный интервал для имеет вид , где квантиль signed-rank Wilcoxon распределения. Этот доверительный интервал не зависит от вида распределения ошибок при условии их симметрии относительно нуля. Отметим, что размах значений множество имеет порядок . Поэтому, для умеренных по размеру выборок, тест Уилкоксона менее зависим от дискретного характера распределения статистики критерия, то есть выбранный уровень значимости в этом случае ближе к найденному.

5. В качестве практического примера рассмотрим данные об объеме продаж (в штуках) для восьми товарных позиций в двух магазинах A и B за неделю. Ответим на вопрос, в каком магазине спрос на товары выше?

Составим выборку, каждый элемент которой представляет собой разницу в продажах соответствующей товарной позиции в магазинах A и B. Пусть характеризует центральное значение выборки. Следующая R сессия показывает результат применения теста Уилкоксона и t-теста для проверки правосторонней гипотезы

> Store_A <- c(82, 69, 73, 43, 58, 56, 76, 65)> Store_B <- c(63, 42, 74, 37, 51, 43, 80, 62)> response <- Store_A - Store_B> wilcox.test(response, alternative = "greater", conf.int = TRUE)Wilcoxon signed rank exact testdata:  responseV = 32, p-value = 0.02734alternative hypothesis: true location is greater than 095 percent confidence interval:   1 Infsample estimates:(pseudo)median           7.75 > t.test(response, alternative = "greater", conf.int = TRUE)One Sample t-testdata:  responset = 2.3791, df = 7, p-value = 0.02447alternative hypothesis: true mean is greater than 095 percent confidence interval: 1.781971      Infsample estimates:mean of x      8.75 

Тест Уилкоксона wilcox.test() возвращает статистику , p-value теста, оценку Ходжеса-Лемана для и доверительный интервал для . Т-тест t.test() имеет аналогичный синтаксис и результаты. Как видно, обе процедуры отвергают нулевую гипотезу на уровне , то есть можно сказать, что спрос на продукцию в магазине A выше.

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

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

Анализ модели о сдвиге параметров положения двух генеральных совокупностей начинается с описания свободной от распределения ранговой процедуры Манна-Уитни-Уилкоксона (Mann-Whitney-Wilcoxon, MWW), здесь строятся точечные и интервальные оценки для величины сдвига. Далее кратко описывается метод анализа, основанный на применении score функций и, с его помощью, также проверяется нулевая гипотеза о величине параметре сдвига. В заключение, модель для параметра положения формулируется в виде регрессионной задачи, решение которой также позволяет построить точечную и интервальную оценки для параметра сдвига.

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

1.Пустьи две непрерывные случайные величины:иобозначают функцию (cdf) и плотность (pdf) распределения случайной величины, аиобозначают соответственно функцию (cdf) и плотность (pdf) случайной величины. Будем говорить, чтои следуют модели для параметра положения (location model), если для некоторого параметра, имеем

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

2.Рассмотрим две независимые друг от друга выборки, извлеченные в соответствии с выбранной моделью. Пусть случайная выборка из генеральной совокупности с законом распределения(с функциями cdf и pdfисоответственно), и случайная выборка из генеральной совокупности с законом распределения(с функциями cdf и pdfисоответственно). Пусть также размер объединённой выборки. Рассмотрим гипотезу

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

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

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

> z <- c(12, 18, 11, 5, 11, 5, 11, 11)> rank(z)[1] 7.0 8.0 4.5 1.5 4.5 1.5 4.5 4.5

Пустьозначает рангсреди элементов объединенной выборки, то есть среди. Тогда статистика теста Уилкоксона имеет вид

На практике часто используется эквивалентнаяследующая статистика Манна-Уитни-Уилкоксона (Mann-Whitney-Wilcoxon, MWW). Рассмотрим набор из всехразностей видаи пусть число положительных таких разностей, то есть

В этом случае, верно равенство

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

4.Оценкой параметра сдвигав тесте Манна-Уитни-Уилкоксона является медиана выборки, составленной из всехпарных разностей вида (оценка Ходжеса-Лемана (Hodges-Lehmann))

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

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

5.Для примера применим тест Манна-Уитни-Уилкоксона к двум выборкам из t-распределения cстепенями свободы и значением параметра сдвига .

> x <- round(rt(11, 5) * 10 + 42, 1)> y <- round(rt(9, 5) * 10 + 50, 1)> x [1] 76.6 41.0 59.3 34.9 29.1 45.0 42.6 31.1 32.4 52.5 47.9> y [1] 58.3 47.2 40.1 45.8 62.0 58.7 64.8 48.1 49.5> wilcox.test(y, x, exact = TRUE, conf.int = TRUE, conf.level = 0.95)Wilcoxon rank sum exact testdata:  y and xW = 72, p-value = 0.09518alternative hypothesis: true location shift is not equal to 095 percent confidence interval: -1.0 18.4sample estimates:difference in location                   10.4

Здесь:с p-valueдоверительный интервалвключает истинное значение, точечная оценка равна, доверительная вероятность. По умолчанию p-value основано на точном распределении статистикидля малых выборок без хвостов. В случае выбора аргумента exact = FALSE и correct = FALSE (без корректировки на непрерывность) используется асимптотический метод, изложенный далее. В этом случае значение p-value равно.

> wilcox.test(y, x, exact = FALSE, correct = FALSE)Wilcoxon rank sum testdata:  y and xW = 72, p-value = 0.08738alternative hypothesis: true location shift is not equal to 0

6.Рассмотрим функцию, где функция(score функция) определена на интервалеи

Например,, где функция, обратная к функции cdf стандартной нормальной случайной величины. В этом случае преобразование посредством функции(Normal score function) ставит в соответствие рангам элементов выборки значения, которые ожидаемы тем, как если бы данные были получены из нормального распределения. Заметим, что в данном определении термин normal score употребляется в смысле rankit, а не в смысле standard score или z-score. Наряду с функцией normal score, возможны другие виды score функций, напримерзадаёт score функцию Уилкоксона.

Для выбранной модели для параметра положения определим следующую функцию переменной:

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

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

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

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

7.В следующей R сессии для данных из рассматриваемого в статье примера представлен расчёт асимптотической статистикии соответствующее ей значение p-value с использованием score функции Уилкоксона (пакет Rfit).

> x <- c(76.6, 41.0, 59.3, 34.9, 29.1, 45.0, 42.6, 31.1, 32.4, 52.5, 47.9)> y <- c(58.3, 47.2, 40.1, 45.8, 62.0, 58.7, 64.8, 48.1, 49.5)> # Объединим выборки x и y>   z = c(x, y)>   n1 = length(x)>   n2 = length(y)>   n = n1 + n2> # Выберем в качестве score функции функцию Уилкоксона>   scores = Rfit::wscores > # Найдём соответствующие score значения для рангов выборки z >   rs = rank(z)/(n + 1)>   asg = Rfit::getScores(scores, rs)> # Найдём значение статистики Sphi при нулевой гипотезе  >   Sphi = sum(asg[(n1 + 1):n])> # Найдём дисперсию Sphi>   asc = Rfit::getScores(scores, 1:n/(n + 1))>   varphi = ((n1 * n2)/(n * (n - 1))) * sum(asc^2)> # Найдём значения zphi и p-value >   zphi = Sphi/sqrt(varphi)>   alternative = "two.sided">   pvalue <-+     switch(+     alternative,+     two.sided = 2 * (1 - pnorm(abs(zphi))),+     less = pnorm(zphi),+     greater = 1 - pnorm(zphi)+   )> # Выведем результаты>   res <- list(Sphi = Sphi, statistic = zphi, p.value = pvalue)>   with(res, cat("statistic = ", statistic, ", p-value = ", p.value, "\n"))statistic =  1.709409 , p-value =  0.08737528

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

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

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

В следующей R сессии представлены результаты ранговой регрессии для имеющихся данных.

> z = c(x, y)> ci <- c(rep(0, n1), rep(1, n2))> fit <- Rfit::rfit(z ~ ci, scores = Rfit::wscores)> coef(summary(fit))            Estimate Std. Error  t.value      p.value(Intercept)     41.8   4.400895 9.498068 1.960951e-08ci              10.4   5.720594 1.817993 8.574801e-02

Таблица с результатами, наряду с самой оценкойсодержит её стандартную ошибку. Используя это, можно, например, построить приблизительныйдоверительный интервал для сдвига, основанный на квантиле уровня t-распределения сстепенями свободы:

> conf.level <- 0.95> estse <- coef(summary(fit))[2, 1:2]> alpha <- 1 - conf.level> alternative = "two.sided"> tcvs <- switch(+   alternative,+   two.sided = qt(1 - alpha / 2, n - 2) * c(-1, 1),+   less = c(-Inf, qt(1 - alpha, n - 2)),+   greater = c(qt(alpha, n - 2), Inf)+ )> conf.int <- estse[1] + tcvs * estse[2]> cat(100 * conf.level, " percent confidence interval:\n", conf.int)95  percent confidence interval: -1.618522 22.41852 

Построенные доверительный интервалнесколько шире своего дискретного аналога , найденного ранее.

Мы рассмотрели два подхода (процедура Манна-Уитни-Уилкоксона и метод score функций) для проверки гипотезы о наличии сдвига у параметров положения двух генеральных совокупностей. Для величины сдвига построены точечная оценка и два типа доверительных интервала.

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

Анализ модели о сдвиге параметров положения двух генеральных совокупностей начинается с описания свободной от распределения ранговой процедуры Манна-Уитни-Уилкоксона (Mann-Whitney-Wilcoxon, MWW), здесь строятся точечные и интервальные оценки для величины сдвига. Далее кратко описывается метод анализа, основанный на применении score функций и, с его помощью, также проверяется нулевая гипотеза о величине параметре сдвига. В заключение, модель для параметра положения формулируется в виде регрессионной задачи, решение которой также позволяет построить точечную и интервальную оценки для параметра сдвига.

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

1.Пустьи две непрерывные случайные величины:иобозначают функцию (cdf) и плотность (pdf) распределения случайной величины, аиобозначают соответственно функцию (cdf) и плотность (pdf) случайной величины. Будем говорить, чтои следуют модели для параметра положения (location model), если для некоторого параметра, имеем

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

2.Рассмотрим две независимые друг от друга выборки, извлеченные в соответствии с выбранной моделью. Пусть случайная выборка из генеральной совокупности с законом распределения(с функциями cdf и pdfисоответственно), и случайная выборка из генеральной совокупности с законом распределения(с функциями cdf и pdfисоответственно). Пусть также размер объединённой выборки. Рассмотрим гипотезу

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

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

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

> z <- c(12, 18, 11, 5, 11, 5, 11, 11)> rank(z)[1] 7.0 8.0 4.5 1.5 4.5 1.5 4.5 4.5

Пустьозначает рангсреди элементов объединенной выборки, то есть среди. Тогда статистика теста Уилкоксона имеет вид

На практике часто используется эквивалентнаяследующая статистика Манна-Уитни-Уилкоксона (Mann-Whitney-Wilcoxon, MWW). Рассмотрим набор из всехразностей видаи пусть число положительных таких разностей, то есть

В этом случае, верно равенство

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

4.Оценкой параметра сдвигав тесте Манна-Уитни-Уилкоксона является медиана выборки, составленной из всехпарных разностей вида (оценка Ходжеса-Лемана (Hodges-Lehmann))

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

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

5.Для примера применим тест Манна-Уитни-Уилкоксона к двум выборкам из t-распределения cстепенями свободы и значением параметра сдвига .

> x <- round(rt(11, 5) * 10 + 42, 1)> y <- round(rt(9, 5) * 10 + 50, 1)> x [1] 76.6 41.0 59.3 34.9 29.1 45.0 42.6 31.1 32.4 52.5 47.9> y [1] 58.3 47.2 40.1 45.8 62.0 58.7 64.8 48.1 49.5> wilcox.test(y, x, exact = TRUE, conf.int = TRUE, conf.level = 0.95)Wilcoxon rank sum exact testdata:  y and xW = 72, p-value = 0.09518alternative hypothesis: true location shift is not equal to 095 percent confidence interval: -1.0 18.4sample estimates:difference in location                   10.4

Здесь:с p-valueдоверительный интервалвключает истинное значение, точечная оценка равна, доверительная вероятность. По умолчанию p-value основано на точном распределении статистикидля малых выборок без хвостов. В случае выбора аргумента exact = FALSE и correct = FALSE (без корректировки на непрерывность) используется асимптотический метод, изложенный далее. В этом случае значение p-value равно.

> wilcox.test(y, x, exact = FALSE, correct = FALSE)Wilcoxon rank sum testdata:  y and xW = 72, p-value = 0.08738alternative hypothesis: true location shift is not equal to 0

6.Рассмотрим функцию, где функция(score функция) определена на интервалеи

Например,, где функция, обратная к функции cdf стандартной нормальной случайной величины. В этом случае преобразование посредством функции(Normal score function) ставит в соответствие рангам элементов выборки значения, которые ожидаемы тем, как если бы данные были получены из нормального распределения. Заметим, что в данном определении термин normal score употребляется в смысле rankit, а не в смысле standard score или z-score. Наряду с функцией normal score, возможны другие виды score функций, напримерзадаёт score функцию Уилкоксона.

Для выбранной модели для параметра положения определим следующую функцию переменной:

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

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

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

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

7.В следующей R сессии для данных из рассматриваемого в статье примера представлен расчёт асимптотической статистикии соответствующее ей значение p-value с использованием score функции Уилкоксона (пакет Rfit).

> x <- c(76.6, 41.0, 59.3, 34.9, 29.1, 45.0, 42.6, 31.1, 32.4, 52.5, 47.9)> y <- c(58.3, 47.2, 40.1, 45.8, 62.0, 58.7, 64.8, 48.1, 49.5)> # Объединим выборки x и y>   z = c(x, y)>   n1 = length(x)>   n2 = length(y)>   n = n1 + n2> # Выберем в качестве score функции функцию Уилкоксона>   scores = Rfit::wscores > # Найдём соответствующие score значения для рангов выборки z >   rs = rank(z)/(n + 1)>   asg = Rfit::getScores(scores, rs)> # Найдём значение статистики Sphi при нулевой гипотезе  >   Sphi = sum(asg[(n1 + 1):n])> # Найдём дисперсию Sphi>   asc = Rfit::getScores(scores, 1:n/(n + 1))>   varphi = ((n1 * n2)/(n * (n - 1))) * sum(asc^2)> # Найдём значения zphi и p-value >   zphi = Sphi/sqrt(varphi)>   alternative = "two.sided">   pvalue <-+     switch(+     alternative,+     two.sided = 2 * (1 - pnorm(abs(zphi))),+     less = pnorm(zphi),+     greater = 1 - pnorm(zphi)+   )> # Выведем результаты>   res <- list(Sphi = Sphi, statistic = zphi, p.value = pvalue)>   with(res, cat("statistic = ", statistic, ", p-value = ", p.value, "\n"))statistic =  1.709409 , p-value =  0.08737528

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

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

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

В следующей R сессии представлены результаты ранговой регрессии для имеющихся данных.

> z = c(x, y)> ci <- c(rep(0, n1), rep(1, n2))> fit <- Rfit::rfit(z ~ ci, scores = Rfit::wscores)> coef(summary(fit))            Estimate Std. Error  t.value      p.value(Intercept)     41.8   4.400895 9.498068 1.960951e-08ci              10.4   5.720594 1.817993 8.574801e-02

Таблица с результатами, наряду с самой оценкойсодержит её стандартную ошибку. Используя это, можно, например, построить приблизительныйдоверительный интервал для сдвига, основанный на квантиле уровня t-распределения сстепенями свободы:

> conf.level <- 0.95> estse <- coef(summary(fit))[2, 1:2]> alpha <- 1 - conf.level> alternative = "two.sided"> tcvs <- switch(+   alternative,+   two.sided = qt(1 - alpha / 2, n - 2) * c(-1, 1),+   less = c(-Inf, qt(1 - alpha, n - 2)),+   greater = c(qt(alpha, n - 2), Inf)+ )> conf.int <- estse[1] + tcvs * estse[2]> cat(100 * conf.level, " percent confidence interval:\n", conf.int)95  percent confidence interval: -1.618522 22.41852 

Построенные доверительный интервалнесколько шире своего дискретного аналога , найденного ранее.

Мы рассмотрели два подхода (процедура Манна-Уитни-Уилкоксона и метод score функций) для проверки гипотезы о наличии сдвига у параметров положения двух генеральных совокупностей. Для величины сдвига построены точечная оценка и два типа доверительных интервала.

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

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

Отметим, что в оговоренных условиях, имеем

откуда, в частности, следует

2.Интересующая нас гипотеза имеет вид

где. Далее рассматривается, основанная на рангах, процедура для проверки указанной гипотезы тест Флигнера-Киллина (Fligner-Killeen test), а также строятся точечная оценка и доверительный интервал для параметра.

3.Так как тест для проверки гипотезыдолжен не зависеть от параметров положения, то следует центрировать имеющиеся наблюдения. В параметрическом F-тесте (var.test {stats}), основанном на отношении выборочных дисперсий, для этой цели используются выборочные средние; вместо них мы применим, равные им в силу симметрии, медианы:

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

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

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

где ранги рассчитываются по объединенной выборке. Данная статистика может быть использована для проверки исследуемой нулевой гипотезыо величине параметра[см. пункты 6, 7 прошлой статьи].

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

Это позволяет использовать стандартизованную статистику

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

Решение вопроса о том, какую score функцию использовать при расчете статистикизависит от вида функции. В распространенном случае, когдаимеет нормальное распределение, оптимальной является квадратичная нормальная score функция (squared-normal score) вида

где cdf стандартного нормального распределения, FK аббревиатура, соответствующая работе Fligner и Killeen (1976) [1], где данная функция была предложена.

5.Для демонстрации изложенной выше методики рассмотрим данные конкретного эксперимента (Doksum and Sievers (1976) [2]) о влиянии озона на прирост в весе у крыс. Контрольная группа крысв течениедней находилась в обычной среде, в то время как экспериментальная группа крысв течение 7 дней находилась в среде с повышенной концентрацией озона. В конце недели изменение в весе каждой крысы были взяты в качестве отклика.

Следующая программа на зыке R содержит алгоритм расчета статистикии соответствующего значения p-value для проверки двухсторонней гипотезы о различии параметров масштаба двух выборок.

> x <-+   c(+     41.0,    38.4,    24.4,    25.9,    21.9,+     18.3,    13.1,    27.3,    28.5,   -16.9,+     26.0,    17.4,    21.8,    15.4,    27.4,+     19.2,    22.4,    17.7,    26.0,    29.4,+     21.4,    26.6,    22.7)> y <-+   c(+     10.1,     6.1,    20.4,     7.3,    14.3,+     15.5,    -9.9,     6.8,    28.2,    17.9,+     -9.0,  - 12.9,    14.0,     6.6,    12.1,+     15.7,    39.9,   -15.9,    54.6,   -14.7,+     44.1,    -9.0)> > # Зададим squared-normal score функцию>   fkscores = new(+     "scores",+     phi = function(u) {+       (qnorm((u + 1) / 2)) ^ 2 - 1+     },+     Dphi = function(u) {+       qnorm((u + 1) / 2) / dnorm(qnorm((u + 1) / 2))+     }+   )> # Проверим гипотезу о равенстве 0 параметра Delta>   zed = c(abs(x - median(x)), abs(y - median(y)))>   n1 = length(x)>   n2 = length(y)>   n = n1 + n2> > # В расчётe используется стандартизированная score функция>   rz = rank(zed)/(n + 1)>   afk = Rfit::getScores(fkscores, rz)> > # Расчет статистики Sfk, zfk, p-value>   Sfk = sum(afk[(n1+1):n])> >   varfk = ((n1 * n2)/(n*(n-1))) * sum(afk^2)>   zfk = Sfk/sqrt(varfk)>   see = "two.sided">   if (see == "two.sided") {+     pvalue = 2 * (1 - pnorm(abs(zfk)))+   }>   if (see == "greater") {+     pvalue = 1 - pnorm(zfk)+   }>   if (see == "less") {+     pvalue = pnorm(zfk)+   }> > # Вывод результатов>   res <- list(statistic = zfk, p.value = pvalue)>   with(res, cat("statistic = ", statistic, ", p-value = ", p.value, "\n"))statistic =  2.095976 , p-value =  0.03608434

Стандартизованная статистика теста Флигнера-Киллина имеет значениеc p-value, то есть можно сказать, что наличие озона влияет на вариабельность в приросте веса у крыс.

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

вектор-индикатор, с первымиэлементами равными нулю и последнимиэлементами равными единице. Тогда регрессионная задача [см. пункт 8 прошлой статьи] имеет вид

Применяя алгоритм построения ранговой регрессии (rfit {Rfit}) с использованием квадратичной нормальной score функцийможно найти оценкупараметра. Откуда оценка для параметраполучается в виде.

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

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

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

> # Подготовим данные   >   xs = abs(x - median(x))>   xs = xs[xs != 0]>   xstarl = log(xs)>   ys = abs(y - median(y))>   ys = ys[ys != 0]>   ystarl = log(ys)>   n1s = length(xs)>   n2s = length(ys)>   ns = n1s + n2s> > # Построим уравнение регрессии>   cvec = c(rep(0, n1s), rep(1, n2s))>   zed = c(xstarl, ystarl)>   fitz = rfit(zed ~ cvec, scores = fkscores)>   (sumf = summary(fitz))Call:rfit.default(formula = zed ~ cvec, scores = fkscores)Coefficients:            Estimate Std. Error t.value   p.value    (Intercept)  1.42288    0.38905  3.6573 0.0007044 ***cvec         0.86586    0.42783  2.0238 0.0493802 *  ---Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1   1Multiple R-squared (Robust): 0.07815811 Reduction in Dispersion Test: 3.56096 p-value: 0.06608 >   delta = coef(sumf)[2, 1]> > # Оценка отношения параметров масштаба>   (eta = exp(delta))[1] 2.377049> > # Построим доверительный интервал>   conf.level = 0.95>   se = coef(sumf)[2, 2]>   tc = abs(qt(((1 - conf.level)/2), ns - 2))>   lb = delta - tc * se>   ub = delta + tc * se>   conf.int = c(exp(lb), exp(ub))>   cat(100 * conf.level, " percent confidence interval:\n", conf.int)95  percent confidence interval: 1.002462 5.636488

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

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

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

• когда распределение генеральной совокупности относится к классу засорённых нормальных распределений (skewed contaminated normal distributions). В этом случае случайная величина получается по закону, гдеиимеюти распределения соответственно, имеет распределение Бернулли с вероятностью успехаи,инезависимы в совокупности.

Ссылки на литературу:

1. Fligner, M. A. and Killeen, T. J. (1976), Distribution-free two-sample test for scale, Journal of the American Statistical Association, 71, 210-213.

2. Doksum, K. A. and Sievers, G. L. (1976), Plotting with confidence: Graphical comparisons of two populations, Biometrika, 63, 421-434.

Сегодня, 4 мая, в день Звездных войн мы подготовили для Вас подробный гайд по основным функциям библиотеки dplyr. Почему именно в день Звездных войн? А потому что разбирать мы все будем на примере датасета starwars.

Ну что, начнем!

Если Вы хотите получать еще больше интересных материалов по программированию, Data Science и математике, то подписывайтесь на нашу группу ВК, канал в Телеграме и Инстаграм. Каждый день мы публикуем полезный контент и вопросы с реальных собеседований.

Кстати говоря, а Вы знаете, почему день Звездных войн отмечается именно 4 мая? Все очень просто - знаменитая фраза May the fource be with you крайне созвучна с May, the 4th, т.е. 4 мая :)

Знакомство с датасетом

Для начала, давайте подключим библиотеку dplyr. Делается это с помощью функции library.

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

1. name - имя или прозвище героя вселенной Звездных войн. Например, Оби-Ван Кеноби.

2. height - высота персонажа

3. mass - масса персонажа

4. hair_color - цвет волос

5. skin_color - цвет кожи

6. eye_color - цвет глаз

7. birth_year - год рождения (до битвы на Явине)

8. sex - биологический пол (есть существа без пола и гермафродиты)

9. gender - поведенческий пол персонажа (например, на какой пол запрограммированы дроиды)

10. homeworld - из какой вселенной существо

11. species - биологический вид

12. films - список фильмов, в которых появилось существо

13. vehicles - список транспорта, которым существо управляло

14. starships - список космических кораблей, которыми существо управляло

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

Знакомство с dplyr

dplyr - это один из основных пакетов семейства tidyverse. Его ближайший аналог в Python - библиотека Pandas. Библиотека dplyr служит для манипуляции с данными: фильтрация, сортировка, создание вычислимых столбцов и так далее.

По своему функционалу библиотека dplyr очень похожа на стандартный синтаксис SQL. Немного ранее мы вместе с Алексеем Селезневым из Netpeak делали карточки: сравнение глаголов dplyr и операторов SQL. Вы можете посмотреть их здесь.

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

1. Каждая переменная находится в отдельном столбце

2. Каждое измерение находится в отдельной строке

3. Каждое значение находится в отдельной ячейке

Кстати говоря, наш датасет starwars не совсем соответствует этим правилам. Нам это не помешает, но сможете ли Вы найти, в чем именно несоответствие? ;)

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

Итак, вернемся к глаголам dplyr. Пакет dplyr позволяет делать несколько основных операций:

• Отбор столбцов

• Фильтрация строк

• Сортировка строк

• Группировка

• Агрегирование

• Создание вычислимых столбцов

• Объединение таблиц

Давайте переходить к практике - хватит теории!

Отбор столбцов

Давайте начнем с самого простого - как выбрать только определенные столбцы из таблицы? Очень просто - с помощью функции select.

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

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

• contains: название столбца содержит

• ends_with: название столбца заканчивается на

• matches: название столбца соответствует регулярному выражению

• num_range: поиск занумерованных столбцов, например, V1, V2, V3...

• one_of: название столбца соответствует одному из вариантов

• starts_with: название столбца начинается с

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

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

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

Обратите внимание, все наши запросы возвращали таблицу tibble. А что, если мы хотим отобрать столбец и сразу начать работать с ним, как с вектором? Мало кто знает, но для этого в dplyr есть специальная функция pull. Она возвращает не таблицу, как остальные глаголы dplyr, а вектор.

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

Фильтрация строк

Фильтрация строк по значениям - это аналог привычного оператора WHERE в SQL. В синтаксисе dplyr же для этого используется глагол filter (как неожиданно, правда?).

Использовать filter очень просто - задаем некое логическое выражение и функция вернет только те строки, для которых это выражение True. Например:

Вы можете комбинировать несколько условий с помощью & и |:

Логические выражения Вы можете конструировать не только с помощью >/<, но и с помощью других логических операторов:

• >=

• <=

• is.na

• !is.na

• %in%

• !

Например:

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

Другая функция, sample_n, отбирает n случайных строк.

Функция slice же, например, позволяет отбирать строки по их индексу:

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

Итак, с фильтрацией строк мы тоже разобрались, переходим к следующему разделу.

Сортировка строк

За сортировку строк в SQL отвечает оператор ORDER BY. В dplyr для этого существует функция arrange.

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

Чтобы отсортировать строки по убыванию, достаточно добавить функцию desc.

А если отсортировать нужно по нескольким столбцам? Легко, просто указываем их названия :)

Кстати говоря, с arrange Вы можете также использовать вспомогательные глаголы, которые мы обсуждали в блоке с select. Для этого нужно использовать функцию across. Например:

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

Группировка и агрегатные функции

Группировка - базовая операция, которая необходима для расчета различных характеристик - средних значений, медиан, сумм, количества строк в группе и так далее. В SQL для этого используется оператор GROUP BY и агрегатные функции sum, min, max и так далее. В dplyr же все то же самое :) Ну, почти.

Давайте для начала сгруппируем наши строки по полю eye_color:

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

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

Обратите внимание, здесь мы дополнительно воспользовались функцией drop_na из пакета tidyr, чтобы удалить строки, в которых есть пропуски. Сделали мы это, чтобы при расчете наших максимальных/минимальных/других агрегатных значений не вылезали значения NA.

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

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

• n_distinct - считает количество уникальных элементов в группе

• last - возвращает последнее значение в группе

• nth - возвращает n-ое значение из группы

• quantile - возвращает заданную квантиль

• IQR - межквартильный размах, inter-quartile range

• mad - медианное абсолютное отклонение, median absolute deviation

• sd - стандартное отклонение

• var - вариация

Хорошо, с базовой группировкой мы разобрались. А что если пойти дальше

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

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

Без проблем - достаточно применить уже известную функцию across.

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

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

Вычислимые столбцы

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

Для начала базовая вещь - создадим вычислимый столбец отношения веса к росту.

А если нам нужно применить одну и ту же функцию к нескольким столбцам? Без проблем - снова нас выручит across. Например, давайте умножим на 10 все столбцы с численным типом данных. При этом мы не будем создавать новые столбцы - мы модифицируем старые.

Видно, что в столбцах с весом, ростом и возрастом все значения умножились на 10.

Давайте и с текстом немного поработаем - ко всем строковым значениям добавим суффикс _new. Для этого нам понадобится библиотека stringr все из того же семейства tidyverse.

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

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

Видно, что в конец таблицы добавился новый столбец rnk с рангами для каждой строки.

Таких функций, на самом деле, масса. Вот некоторые из них:

• lag

• lead

• cumsum

• dense_rank

• ntile

• row_number

• case_when

• coalesce

Это, пожалуй, самые популярные и все они на 100% перекликаются с SQL. Приведем несколько примеров.

Ну что, давайте переходить к последнему пункту - к объединению таблиц.

Объединение таблиц

Пожалуй, последняя необходимая для полноценной работы с датафреймами операция - объединение таблиц. Операция объединения в dplyr тесно связана с джоинами в SQL. Вот перечень основных функций:

• left_join

• right_join

• inner_join

• full_join

Аналогичные глаголы присутствуют и в стандартном SQL. Давайте создадим новую таблицу - возьмем датасет starwars и оставим только те строки, где вид существа - человек. А после попробуем различные виды джоинов.

Обратите внимание, мы воспользовались функцией rename и переименовали поле name. Сделали мы это намеренно, чтобы показать работу аргумента by (аналог ON в SQL) для связи столбцов при джоине.

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

Если сделать full_join аналогичным образом, то строк в итоговой таблице будет 87, т.к. в таблице starwars их 87.

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

1 вариант:

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

2 вариант:

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

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

• bind_rows - помещает одну таблицу под другой

• bind_cols - ставит одну таблицу справа от другой

• intersect - находит пересекающиеся строки

• setdiff - разность таблиц, т.е. строки из первой таблицы, которых нет во второй

• union - возвращает строки, которые есть в любой из таблиц (дубликаты исключаются)

• union_all - аналогично union, но оставляет дуликаты

Эпилог

Мы с Вами рассмотрели все основные операции библиотеки dplyr и почти все основные функции. Все остальное - практика, практика и еще раз практика. Если у Вас будут вопросы - рады будем помочь и ответить в комментариях :)

Если Вы хотите получать еще больше интересных материалов по программированию, Data Science и математике, то подписывайтесь на нашу группу ВК, канал в Телеграме и Инстаграм. Каждый день мы публикуем полезный контент и вопросы с реальных собеседований.

А, и совсем забыли. May the fource be with you!

P.S. Здесь Вы можете найти официальную шпаргалку по всем функциям dplyr. Все удобно и компактно собрано в одном месте :)

Зачем нужен Rust в R?

Первый вопрос, который должен возникнуть у читателя -- а зачем вообще использовать Rust вместе с R? Ответ довольно прост: Rust -- новый системный язык программирования, спроектированный специально для написания безопасного и легко распараллеливаемого кода. Rust довольно сложен в освоении (в сравнении с другими языками), но при этом предоставляет отличные инструменты для разработки. Rust имеет довольно неплохую ООП систему и очень много заимствует из функциональных языков программирования. Несмотря на дополнительную сложность из-за функциональных/ООП компонентов, Rust позиционируется как zero-cost abstraction язык, так же как и C++.

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

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

Что нужно, чтобы R код мог вызвать Rust-библиотеку?

На самом деле -- не так уж много. R-пакеты могут содержать директорию src/, в которой находится исходный код на одном из компилируемых языков. С помощью src/Makevars или src/Makevars.win файлов (вариация make) можно контролировать процесс сборки, например, вызвав на одном из шагов cargo (см. пример здесь):

cargo build --release --manifest-path=rustlib/Cargo.toml

При этом Rust -библиотека должна собираться как crate-type = ["staticlib"]. Кроме непосредственной компиляции Rust-кода, нужно предоставить C-обертки к экспортируемым функциям, а так же добавить несколько магических вызовов специальных R-функций, которые объясняют R, какие именно функции и какого типа экспортируются из данной библиотеки (например, вот так).

Основная проблема -- C-обертки и преобразование типов из R SEXP (фактически, специальный указатель) во что-то, совместимое с Rust, учитывая при этом специфику управления памятью в R (все эти ваши PROTECT, UNPORTECT, и т. д.). Как результат -- легко создать примитивный прототип без функционала, практически невозможно написать достаточно большой проект.

Интегрируем R и Rust: три простых шага

Шаг первый: баиндинги для заголовочных файлов R

Взаимодействие с R происходит через специализированный API, доступный обычно ввиде C/ C++ заголовочных файлов (см. $R_HOME\include\). Разумеется, вызывать эти методы можно практически из любого языка, но это неудобно -- загловочные файлы невозможно подключить напрямую к Rust. К счастью, у этой проблемы уже давно есть решение: rust-bindgen (rust-lang/rust-bindgen). bindgen позволяет автоматически генерировать Rust-обертки из заголовочных файлов, и делает это довольно эффективно.

Так появился крейт libR-sys, который предоставляет баиндинги ко всем необходимым внутренним R функциям. Генерация баиндингов -- вещь нетривиальная, bindgen зависит от clangи сложен в конфигурировании, поэтому мы предоставляем pre-computed (заранее сгененрированные) баиндинги для большинства платформ, поддерживающих R. Список включает в себя linux-x64 (созданный с помощью Ubuntu-20.04), win-x86/x64 (с помощью msys2, x86 может иметь проблемы в каких-то пограничных случаях), macOS включая 11 версию (по возможности), x64 и экспериментально arm64 (честно я не знаю, есть ли arm64 сборка R под macOS). Для каждой из упомянутых платформ/архитектур мы стараемся предоставить три версии баиндингов: oldrel, release, и devel, что соответствует "прошлой", "текущей" (сейчас это 4.1.0) и "находящейся в разработке" версиям R.

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

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

Шаг второй: автоматизируем преобразование типов и экспорт функций

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

Прежде чем продолжить, я хочу сделать небольшое отступление. Вся идея проекта extendr и, в особенности, имплементация большей части Rust-крейтов, принадлежит Энди Томасону (@andy-thomason). Без его вклада, на мой субъективный взгляд, extendr в том виде, в котором он существует сейчас, был бы невозможен.

Вернемся обратно к коду. Как избавиться от боилерплейта? Легко, надо всего лишь распарсить исходный код Rust. Например, используя syn и подобные крейты. Моей экспертизы недостаточно, чтобы детально описать процесс парсинга и кодогенерации, но для конечного пользователя экспорт Rust функции становится невероятно простым. Во-первых, нужно пометить функции с помощью аттрибута #[extendr]:

#[extendr]fn add_i32(x : i32, y : i32) -> i32 { x + y }#[extendr]fn add_vec(x : &[i32], y : &[i32]) -> Vec<i32> {     x.iter().zip(y.iter()).map(|v| v.0 + v.1).collect()}

Во-вторых, нужно явно объявить экспортируемые функции:

extendr_module! {  mod extendrtest;  fn add_i32;  fn add_vec;}

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

К сожалению, остается одно небольшое ограничение при интеграции Rust-кода в проект. Дело в том, что если в папке src/ отсутствуют файлы-исходники, то стандартная процедура компиляции R попросут игнорирует все остальное и библиотека не компилируется. Чтобы обойти это, в src/ добавляется единственный файл entrypoint.c, примерно следующего содержания:

void R_init_extendrtest_extendr(void *dll);void R_init_extendrtest(void *dll) {  R_init_extendrtest_extendr(dll);}

Здесь R_init_extendrtest_extendr генерируется автоматически с помощью Rust-крейта, а R_init_extendrtest -- непосредственно вызывается из R. Мы пока что не нашли способа избавиться от этого ограничения.

Некоторых изменений требуют и Makevars-файлы. Вот пример из одного из тестовых проектов:

LIBDIR = ./rust/target/releaseSTATLIB = $(LIBDIR)/libextendrtest.aPKG_LIBS = -L$(LIBDIR) -lextendrtestall: C_clean$(SHLIB): $(STATLIB)$(STATLIB):cargo build --lib --release --manifest-path=./rust/Cargo.tomlC_clean:rm -Rf $(SHLIB) $(STATLIB) $(OBJECTS)clean:rm -Rf $(SHLIB) $(STATLIB) $(OBJECTS) rust/target

Фактически, мы отдельно компилируем Rust-крейт, а потом создаем совместимую с R библиотеку используя Rust-библиотеку и результат компиляции entrypoint.c.

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

extendr выполняет не только кодогенерацию на стороне Rust, он еще генерирует обертки на стороне R. Если предположить, что приведенный выше Rust код является частью пакета {extendrtest}, то становятся досутпны следующие функции:

extendrtest::add_i32(4L, 11L)# [1] 15extendrtest::add_vec(1:10, 10:1)#  [1] 11 11 11 11 11 11 11 11 11 11

Да, насктолько просто.

Шаг третий: user-friendliness

В своей работе мы вдохновлялись такими проектами как {cpp11} - header-only пакет для интеграции C++11 кода. Так появился на свет {rextendr}, R - пакет без Rust-зависимости, который решает три основные задачи:

• Создание шаблона пакета, использующего extendr, наподобие {usethis};

• Компиляция и исполнение Rust - кода на лету, прямо в R-сессии. Именно это демонстрирует Анимация Для Привлечения Внимания;

• Предоставление специальных knitr-модулей (engines), а именно {extendr} и {extendrsrc}, которые позволяют включать фрагменты Rust-кода (и результаты его выполнения) в ваш Rmarkdown прямо рядом с R-кодом, обеспечивая их взаимодействие.

Сейчас {rextendr} отправился на проверку в CRAN, и мы ждем результатов. Я думаю, самое время продемонстроровать несколько примеров именно с использованием {rextendr}. Сразу оговорюсь, для упрощения и воспроизводимости, я буду использовать{reprex}.

Самый простой пример это, конечно же,

rextendr::rust_function("fn hello_r() -> &'static str { \"Hello R!\" }")#> i build directory: 'C:\Users\...\AppData\Local\Temp\Rtmp259cVM\file10186cb44264'#> v Writing 'C:/Users/.../AppData/Local/Temp/Rtmp259cVM/file10186cb44264/target/extendr_wrappers.R'.hello_r()#> [1] "Hello R!"

Пример со сложением я уже показывал, но что будет, если явно передать NA?

rextendr::rust_function("fn add_i32(x : i32, y : i32) -> i32 { x + y }")#> i build directory: 'C:\Users\...\AppData\Local\Temp\Rtmp2P2cnQ\file2f7c65e8269a'#> v Writing 'C:/Users/.../AppData/Local/Temp/Rtmp2P2cnQ/file2f7c65e8269a/target/extendr_wrappers.R'.add_i32(42L, NA)#> Error in add_i32(42L, NA): unable to convert R object to primitive

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

rextendr::rust_function("fn add_i32_opt(x : Option<i32>, y : Option<i32>) -> Option<i32> {    match (x, y) {        (Some(a), Some(b)) => Some(a + b),        _ => None    }}")#> i build directory: 'C:\Users\...\AppData\Local\Temp\Rtmpyg3uPw\file6587a897d2a'#> v Writing 'C:/Users/.../AppData/Local/Temp/Rtmpyg3uPw/file6587a897d2a/target/extendr_wrappers.R'.add_i32_opt(NA, 42L)#> [1] NAadd_i32_opt(42L, 100L)#> [1] 142

Хотите еще больше магии? Макрос R! выполняет внутри R-код, возвращая результат если операция была успешной. Как насчет

x <- 42L y <- 100Lrextendr::rust_eval("R!(x)? * 2 + R!(y)? * 3")#> i build directory: 'C:\Users\...\AppData\Local\Temp\RtmpKeC23J\file32ec53677fc9'#> v Writing 'C:/Users/.../AppData/Local/Temp/RtmpKeC23J/file32ec53677fc9/target/extendr_wrappers.R'.#> [1] 384

Можно попробовать смешать переменные из R и Rust.

library(tibble)x <- 10:1 # Эта переменная на стороне Rrextendr::rust_eval("call!(\"tibble\", x = R!(x), y = 1..=10)")#> i build directory: 'C:\Users\...\AppData\Local\Temp\RtmpcDWhlk\file45802f52dc5'#> v Writing 'C:/Users/.../AppData/Local/Temp/RtmpcDWhlk/file45802f52dc5/target/extendr_wrappers.R'.#> # A tibble: 10 x 2#>        x     y#>    <int> <int>#>  1    10     1#>  2     9     2#>  3     8     3#>  4     7     4#>  5     6     5#>  6     5     6#>  7     4     7#>  8     3     8#>  9     2     9#> 10     1    10

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

Для безопасной печати в Rout существует отдельный макрос : rprintln!.

x <- 42Lrextendr::rust_eval("rprintln!(\"Hello from Rust! x = {}\", R!(x)?.as_integer().unwrap());")#> i build directory: 'C:\Users\...\AppData\Local\Temp\RtmpWQh3w0\file48e024f161ce'#> v Writing 'C:/Users/.../AppData/Local/Temp/RtmpWQh3w0/file48e024f161ce/target/extendr_wrappers.R'.#> Hello from Rust! x = 42

Пишем свой extendr-пакет

В этом разделе я просто приведу пример генерации пакета с использование {rextendr} и других стандартных инструментов:

pkg <- file.path(tempfile(), "myextendr")dir.create(pkg, recursive = TRUE)usethis::create_package(pkg)usethis::proj_activate(pkg)rextendr::use_extendr()rextendr::document()rextendr::document()hello_world()

pkg <- file.path(tempfile(), "myextendr")dir.create(pkg, recursive = TRUE)usethis::create_package(pkg)#> v Setting active project to 'C:/Users/.../AppData/Local/Temp/RtmpAVW4HZ/file122c180d1953/myextendr'#> v Creating 'R/'#> v Writing 'DESCRIPTION'#> Package: myextendr#> Title: What the Package Does (One Line, Title Case)#> Version: 0.0.0.9000#> Authors@R (parsed):#>     * First Last <first.last@example.com> [aut, cre] (YOUR-ORCID-ID)#> Description: What the package does (one paragraph).#> License: `use_mit_license()`, `use_gpl3_license()` or friends to pick a#>     license#> Encoding: UTF-8#> LazyData: true#> Roxygen: list(markdown = TRUE)#> RoxygenNote: 7.1.1#> v Writing 'NAMESPACE'#> v Setting active project to '<no active project>'usethis::proj_activate(pkg)#> v Setting active project to 'C:/Users/.../AppData/Local/Temp/RtmpAVW4HZ/file122c180d1953/myextendr'#> v Changing working directory to 'C:/Users/.../AppData/Local/Temp/RtmpAVW4HZ/file122c180d1953/myextendr/'rextendr::use_extendr()#> v Creating 'src/rust/src'.#> v Writing 'src/entrypoint.c'#> v Writing 'src/Makevars'#> v Writing 'src/Makevars.win'#> v Writing 'src/.gitignore'#> v Writing 'src/rust/Cargo.toml'.#> v Writing 'src/rust/src/lib.rs'#> v Writing 'R/extendr-wrappers.R'#> v Finished configuring extendr for package myextendr.#> * Please update the system requirement in 'DESCRIPTION' file.#> * Please run `rextendr::document()` for changes to take effect.rextendr::document()#> i Generating extendr wrapper functions for package: myextendr.#> ! No library found at 'src/myextendr.dll', recompilation is required.#> Re-compiling myextendr#>   -  installing *source* package 'myextendr' ...#>      ** using staged installation#>      ** libs#>      rm -Rf myextendr.dll ./rust/target/x86_64-pc-windows-gnu/release/libmyextendr.a entrypoint.o#>      "C:/rtools40/mingw64/bin/"gcc  -I"C:/PROGRA~1/R/R-41~1.0/include" -DNDEBUG          -O2 -Wall  -std=gnu99 -mfpmath=sse -msse2 -mstackrealign  -UNDEBUG -Wall -pedantic -g -O0 -c entrypoint.c -o entrypoint.o#>      cargo build --target=x86_64-pc-windows-gnu --lib --release --manifest-path=./rust/Cargo.toml#>              Updating crates.io index#>             Compiling winapi-build v0.1.1#>       Compiling winapi v0.3.9#>       Compiling winapi v0.2.8#>       Compiling proc-macro2 v1.0.27#>       Compiling unicode-xid v0.2.2#>       Compiling syn v1.0.73#>       Compiling extendr-engine v0.2.0#>       Compiling lazy_static v1.4.0#>             Compiling kernel32-sys v0.2.2#>             Compiling quote v1.0.9#>             Compiling extendr-macros v0.2.0#>             Compiling libR-sys v0.2.1#>             Compiling extendr-api v0.2.0#>             Compiling myextendr v0.1.0 (C:\Users\...\AppData\Local\Temp\RtmpAVW4HZ\file122c180d1953\myextendr\src\rust)#>              Finished release [optimized] target(s) in 33.09s#>      C:/rtools40/mingw64/bin/gcc -shared -s -static-libgcc -o myextendr.dll tmp.def entrypoint.o -L./rust/target/x86_64-pc-windows-gnu/release -lmyextendr -lws2_32 -ladvapi32 -luserenv -LC:/PROGRA~1/R/R-41~1.0/bin/x64 -lR#>      installing to C:/Users/.../AppData/Local/Temp/RtmpAVW4HZ/devtools_install_122c37bd1965/00LOCK-myextendr/00new/myextendr/libs/x64#>   -  DONE (myextendr)#> v Writing 'R/extendr-wrappers.R'.#> i Updating myextendr documentation#> i Loading myextendr#> Writing NAMESPACE#> Writing NAMESPACE#> Writing hello_world.Rdrextendr::document()#> i Generating extendr wrapper functions for package: myextendr.#> i 'R/extendr-wrappers.R' is up-to-date. Skip generating wrapper functions.#> i Updating myextendr documentation#> i Loading myextendr#> Writing NAMESPACE#> Writing NAMESPACEhello_world()#> [1] "Hello world!"

hello_world() написана на Rust и автоматически экспортируется в R. Обратите внимание, что hello_world.Rd был создан при вызове rextendr::document() (аналог devtools::document()). Дело в том, что rextendr-парсер воспринимает /// комментарии как R комментарии. Rust функция выглядит вот так

/// Return string `"Hello world!"` to R./// @export#[extendr]fn hello_world() -> &'static str {    "Hello world!"}

Что автоматически генерирует R обертку

#' Return string `"Hello world!"` to R.#' @exporthello_world <- function() .Call(wrap__hello_world)

и, как результат, обновляет документацию и NAMESPACE с помощью {roxygen2}.

Если этого мало

Здесь я хотел бы коротко описать последнюю важную фичу extendr. Крейт позволяет экспортировать не просто функции, а целые типы. Легким движением руки можно пробросить кастомный тип из Rust в R , а инстансы этого типа -- передавать в обе стороны как ссылки. Это позволяет заполучить ООП в R в традиционном (object-first) стиле, модицифируя in-place объекты, созданные и доступные из Rust:

rextendr::rust_source(code = "struct Counter {    n: i32,}#[extendr]impl Counter {    fn new() -> Self {        Self { n: 0 }    }        fn increment(&mut self) {        self.n += 1;    }        fn get_n(&self) -> i32 {        self.n    }}")#> i build directory: 'C:\Users\...\AppData\Local\Temp\RtmpWOu1pt\file5318783e2176'#> v Writing 'C:/Users/.../AppData/Local/Temp/RtmpWOu1pt/file5318783e2176/target/extendr_wrappers.R'.cntr <- Counter$new()cntr$get_n()#> [1] 0cntr$increment()cntr$increment()cntr$get_n()#> [1] 2

Вместо заключения

Статья получилась гораздо длиннее и сумбурней, чем я ожидал. Тем не менее, я не успел описать все возможности extendr. Этот проект амбициозный и еще далек от завершения, но я считаю, что давно пришло время добавить поддержку Rust в R, а главное сделать взаимодействие этих языков удобным. Мы осторожно надеемся, что в конечном итоге сможем добавить официальную поддержку Rust, на равне с C / C++. К сожалению, сейчас ее отсутствие накладывает на нас некоторые ограничения.

Отдельным вызовом было заставить эту систему работать на Windows. Мы столкнулись со множеством проблем, но на данный момент нам удалось справиться практически со всеми трудностями. Для запуска на Windows extendr требует стандартный Rust - тулчейн, stable-x86_64-pc-windows-msvc, с дополнительными целями (targets) x86_64-pc-windows-gnu и i686-pc-windows-gnu, а также Rtools40v2 (последняя версия на момент написания, отличается от Rtools40).

Скудную документацию можно найти здесь и в репозиториях проекта extendr.

Спасибо что дочитали до конца!

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

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

Зачем вообще что-то интегрировать? Почему нельзя просто использовать R (Julia)?

Кратко рассмотрим преимущества обоих языков.

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

Julia тоже удобна, тоже имеет приятный синтаксис, немного более строга к пользователю, но имеет свою особенную стройность. Из уникальных особенностей стоит отметить концепцию множественной диспетчеризации, возможность лёгкого доступа к генерируемому низкоуровневому коду, вплоть до ассемблерного, удобную систему макросов. Если бы этим качества и ограничивались, то она была бы просто аналогом R, но Julia при этом очень производительна, даже по меркам компилируемых языков. Вместе с С, C++ и Fortran Julia входит в престижный клуб языков с петафлопсной производительностью. Код на Julia можно запускать как на обычных многоядерных CPU, так и на GPU или вычислительных кластерах. Язык активно используется при проведении научных симуляций, обработке больших массивов данных, задачах машинного обучения. По сравнению с R аналогичные программы на Julia почти всегда в десятки-сотни раз быстрее. При этом благодаря JIT-компиляции команды выполняются почти на лету, как правило без серьёзных задержек. Не так молниеносно, как в R, задержки иногда всё же возникают, Julia как бы берёт некоторую паузу на подумать, особенно когда по ходу дела вдруг надо скомпилировать какой-нибудь пакет. Но подумав и перейдя к вычислениям, Julia работает с очень высокой производительностью. А существенные задержки, во-первых, случаются не так часто и поэтому не слишком раздражают, и к тому же при желании их можно максимально сократить, а во-вторых, ведётся работа над тем, чтобы повысить отзывчивость исполнения команд. Из-за сочетания производительности и удобства использования популярность Julia стремительно набирает обороты. В рейтингах TIOBE и PYPL язык уверенно движется к двадцатке лидеров и имеет вполне реальные шансы войти в неё уже в этом году. С августа 2018 года, когда вышла версия 1.0 и язык в значительной степени стабилизировался, прошло уже достаточно времени, чтобы он перестал ассоциироваться с чем-то сырым и постоянно меняющимся, появилась какая-то литература по актуальной версии. С другой стороны, многие рецепты из сети просто не работают, так как были опубликованы во времена старых версий и язык с тех пор уже претерпел многие изменения. Пока что по современным версиям языка очень мало литературы (на русском по версиям 1.0 и выше печатных изданий пока что вообще нет), нет такой прозрачной системы документации, как у R, далеко не для всех проблем можно найти в сети готовые красивые решения. С Julia постоянно чувствуешь себя первооткрывателем. Но это тоже приятное ощущение.

Что может дать использование Julia пользователям R

Julia явно быстрее, поэтому на ней разумно производить все ресурсоёмкие вычисления. Есть рецепты, как делать в R вставки на Fortran или C++, но по сравнению с ними Julia имеет гораздо более близкий к R синтаксис, вплоть до прямого заимствования отдельных команд из R. При этом по производительности она находится где-то на уровне близком к C и как правило обгоняет Fortran. В отличие от Fortran, С и C++ код в Julia компилируется на лету. Это позволяет сохранить высокую интерактивность работы, привычную после опыта в R. Вплоть до того, что в Atom и VS Code наиболее популярных средах разработки для Julia, можно точно так же, как в RStudio (наиболее популярной среде разработки для R) отправлять на исполнение произвольные куски кода по Ctrl+Enter. То есть общий подход к работе и ощущение от работы остаются примерно такими же, как и в R, но становится доступна очень высокая производительность.

У Julia есть потенциал, чтобы стать универсальным языком учёных будущего, каким раньше был Fortran. Похоже, ставка авторов Julia на сочетание скорости, удобства и открытости вполне оправдалась. Современный R несмотря на свою приятность сильно ограничен максимально возможной производительностью. Ради оптимизации под производительность приходится отказываться от циклов и писать все ресурсоёмкие процедуры в векторизованном стиле, из-за чего код становится тяжёлым для чтения и отладки. Но даже в максимально оптимизированном виде производительность далеко отстаёт от компилируемых языков. Ведётся работа над повышением производительности, но результат как правило нужен уже сейчас, поэтому приходится искать инструменты, которые прямо сейчас позволяют достичь его. Julia позволяет полноценно использовать циклы, предлагает простые средства для организации параллельных вычислений. Интересной особенностью экосистемы Julia являются интерактивные блокноты Pluto, позволяющие создавать красиво свёрстанные интерактивные научные презентации c формулами в LaTeX и выполняющимися на лету компьютерными симуляциями, меняющими своё поведение в зависимости от изменяемых в ходе показа параметров. Наверное, как-то так должны выглядеть научные публикации будущего.

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

Что может дать использование R пользователям Julia

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

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

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

Что для этого нужно сделать

Очевидно, есть два подхода к решению проблемы: вызывать код R из Julia и вызывать код Julia из R. Для вызова R из Julia есть специальный пакет Rcall. Для вызова Julia из R существует сразу несколько альтернативных решений: пакеты XRJulia, JuliaCall, RJulia и JuliaConnectoR. На данный момент наиболее прогрессивным из них выглядит JuliaConnectoR. Предлагаемый им подход немного тяжеловесен, и ориентирован скорее на использование отдельных функций из Julia, чем на полностью смешанное использование двух языков. Довольно любопытны доводы авторов в польу использования именно R, как связующей среды, опубликованные в их статье на arxiv.org.

Изначально я как раз планировал вызвать Julia из R. Потратил много времени, но так и не смог добиться работоспособности этой конструкции ни с XRJulia, ни с JuliaCall. Какой-либо информации про JuliaConnectoR тогда ещё не было, поэтому не было уверенности, что сама идея совместного использования этих двух языков - здравая. Всё шло к тому, что тяжёлые расчётные части придётся переписать с R на Fortran, с которым соприкасался когда-то в институте, но с тех пор никогда больше не использовал. Просто ради интереса решил попробовать вызвать R из Julia, и каково же было моё удивление, когда всё заработало практически из коробки и оказалось действительно удобным способом совместного использования двух языков. Возможно, это мой субъективный опыт, но именно на этом способе мне хочется сейчас остановиться. Про вызов Julia из R, возможно, расскажу как-нибудь в другой раз.

Итак, что нам понадобится:

1) Установленная Julia

2) Установленный R

3) В Julia необходимо установить пакет Rcall: Pkg.add("RCall")

4) Необходимо прописать путь до R. В Windows это удобнее всего сделать, задав переменную среды. В моём случае: R_HOME = C:\Program Files\R\R-4.0.3

Ну вот, собственно, и всё. Наш волшебный гибрид готов! Теперь в любом месте внутри кода Julia можно вызывать код на R и это будет работать. Можно даже не закрывать RStudio с запущенным сеансом в R. Вызов из Julia будет обрабатываться независимо от него.

Подключаем установленную библиотеку Rcall:

using RCall

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

R"plot(rnorm(10))"

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

R"var <- 5*8"

R"var"

R"var / 100"

Фактически это и есть сеанс в R, только обёрнутый в команды Julia. На каком языке в данном случае ведётся работа это скорее уже философский вопрос.

Естественно, дальше нам захочется как-то передавать переменные из одного языка в другой. Для этого удобно использовать команды @rput и @rget. Забрасываем какие-то переменные в R, что-то там считаем, забираем результаты обратно:

a = 1

@rput a

R"b <- a*2 + 1"

@rget b

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

R"""

"""

Например:

a = 1

@rput a

R"""

b <- a*2 + 1

c <- b*3

"""

@rget c

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

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

Есть ли здесь какие-то подводные камни? Я нашёл пока только два:

В текущей версии не поддерживаются кириллические комментарии в R. Придётся их поудалять или заменить на латинские.

В интерактивных блокнотах Pluto, силе и гордости экосистемы Julia, поддержка R работает странно. При первом запуске всё работает, но при изменении входных данных или содержимого R скрипта вместо обновления порождаемых им объектов происходит разрыв коннекта с R. Я пока что не нашёл способа, как обойти это ограничение.

Ещё небольшой комментарий: для того, чтобы комфортно работать с подобными двуязычными проектами, необходимо, чтобы среда разработки поддерживала оба языка и их совместное использование. Как Atom (прошлый фаворит Джулии, с которым она рассталась летом), так и VS Code (нынешний фаворит) при установке синтаксических дополнений для обоих языков такой режим поддерживают. Но очень может быть, что в каких-то средах такой поддержки не будет. Это может стать проблемой, потому что осуществлять навигацию по простыне кода, оформленного как один огромный комментарий, не очень удобно.

На этом всё! Надеюсь, это небольшая заметка была полезна

(при условии нормальности распределения)

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

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

Проведенное предварительное исследование (текст исследования доступен на GitHub) показало, что в зависимости от конкретной комбинации значений средних, дисперсии, и особенностей постановки задачи лучшим может быть чуть ли не любой из тестов, рассмотренных в статье "Cavus, M., Yazici, B. Testing the equality of normal distributed and independent groups' means under unequal variances by doex package / The R Journal. 2020. 2 (12). P. 134-155".

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

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

library("tibble")

library("AER")

library("WRS2")

library("doex")

data("GrowthDJ")

XX<-na.omit(GrowthDJ)

library("psych")

describeBy(XX$gdpgrowth, XX$inter)

shapiro.test(XX[XX$inter=='yes',6])

shapiro.test(XX[XX$inter=='no',6])

Получаем, что наши данные распределены нормальны - значит, тесты можно применять

Методика проверки

1. Задается два средних значения и два значения дисперсии (исходя из имеющихся данных по группам)

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

3. Далее с использованием выбранного теста по объединенным данным первой и третьей выборки проверяется гипотеза о равенстве средних на уровне значимости в 0.01. Если полученное p-значение критерия больше 0.01, то гипотеза отклоняется верно, если меньше 0.01 то гипотеза ошибочно принята. По объединенным данным первой и второй выборки также проверяется данная гипотеза. Если полученное p-значение критерия меньше 0.01, то гипотеза принимается верно, если больше 0.01 то гипотеза ошибочно отклонена. После этого данная симуляция проводится 100 раз по каждому тесту, и подводится количество разных исходов.

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

• accuracy (процент правильно сделанных выводов);

• selectivity (доля правильных выводов для выборок, в которых средние не равны);

• precision (доля правильно сделанных выводов о равенстве средних);

• recall (доля правильных выводов для выборок, в которых средние равны);

• FOR (доля ошибочно сделанных выводов о неравенстве средних);

• F-мера (гармоническое среднее precision и recall, общая характеристика точности качества).

Результаты расчетов (код представлен в файле Итого.R на гитхабе) позволили получить следующую таблицу

Выведем лучшие тесты

Таким образом, получаем:

• если нам нужен самый точный тест вообще, то необходимо использовать AF или FA-тест (они лучшие и по точности, и по значению F-score

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

• если необходимо максимально избегать ложно-положительного вывода (т.е. ложного вывода о равенстве средних), можно выбрать любой из 8 тестов (AF,BA,CF,FA,JF,MBF,SS,WA)

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

• если нам нужна максимальная точность вывода о неравенстве средних, то также необходимо применять RGF-тест

Итого - давайте применим AF-тест (Approximate F-test)

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

Проблематика

Современная мета управления продуктом подразумевает управление на основании данных.

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

Вводная часть

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

Первое, что приходит в голову посмотреть, что там с экономикой продукта. И как вообще ведет себя продукт.

Еще через пару дней совещание, на котором вас просят ответить на несколько вопросов:

1. Сейчас мы платим за привлечение клиента 695494. Какая стоимость привлечения для нас является приемлемой? Имеет ли смысл повысить стоимость привлечения на клиента, чтобы получить больший объем?

2. Насколько здоровой выглядит экономика портфеля и какая динамика здесь и сейчас?

3. Мы недавно изменили подход к размеру выдачи и стали в первые кредиты выдавать меньшие чеки. Как это сказалось на продукте?

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

Но есть проблема. Метрики, конечно, все посчитаны: LTV, CAC и прочее. Есть финансовая отчетность, где видны все расходы и доходы, и видно, что продукт в небольшом операционном плюсе.

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

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

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

Анализируй это

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

Для того, чтобы анализировать, понадобятся следующие инструменты: R ,Rstudio,dbeaver(или аналог). Далее рассматривается реальный пример, как могут выглядеть эти данные и что с ними можно сделать.

Итак, транзакции. Напишем в базу нечто вроде select * from transactions t и посмотрим результат.

Rows: 2,226,532Columns: 10$ borrower_id          2, 2, 2, 6, 6, 12, 12, 12, 12, 16, 20, 20, 20, 20, 22, 23, 23, 33, 33, 39, 39, 36, 36$ con_id               1, 1, 1, 2, 2, 4, 4, 4, 4, 5, 7, 7, 7, 7, 8, 9, 9, 10, 10, 11, 11, 12, 12, 12, 12, 13$ disbursement_date    2017-11-23, 2017-11-23, 2017-11-23, 2017-11-24, 2017-11-24, 2017-11-27, 2017-11-27, $ prolongations_count  1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1$ loan_type            "pdl", "pdl", "pdl", "pdl", "pdl", "pdl", "pdl", "pdl", "pdl", "pdl", "pdl", "pdl", "$ date                 2017-12-15, 2017-11-23, 2017-12-06, 2017-11-24, 2017-12-25, 2017-12-29, 2017-11-27, $ type                 "Payments::Transaction::ContractAddTransaction", "Payments::Transaction::DisburseTran$ amount               250000, 1000000, 1200000, 2500000, 3500000, 1040000, 1500000, 10000, 1470000, 1500000$ id                   325, 2, 127, 5, 587, 557500, 557499, 557504, 557507, 17, 182865, 182874, 182869, 1828$ deleted_at           NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N

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

У нас есть идентификаторы заемщика, контракта, и транзакции. Дата выдачи кредита (disbursment_date). Отметка о пролонгации кредита (prolongations_count) - перенос даты выплаты за отдельный платеж. Размер и дата транзакции. Отметка о "мягком" удалении.

Главное - тип транзакции(ContractAdd- возврат денег заемщиком, DisburseTransaction - выплата денег заемщику).

Теперь поговорим об экономике продукта. Принципиально кредитный бизнес работает просто(как и большинство бизнесов): заплатили, чтобы привести пользователя -> выдали кредит -> собрали деньги -> выдали еще кредит.

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

И в нашей табличке есть все данные, чтобы это посчитать.

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

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

• для любителей R, данных в табличке прилично, больше 2 млн. Для таких датасетов предпочитаю data table. Но можно и на tidyverse(по запросу покажу как)

Теперь посмотрим на структуру данных и займемся обработкой.

Чуть обработаем их и приведем данные в нужный нам формат. Назовем транзакции удобно и сделаем отрицательными транзакции выдачи кредитов (поля z_type и am):

lk %>% data.table()->lk1 #Создаем отдельный объект вида data tablelk1[,':='(z_type=z_type<-fifelse(#  создаем отдельную переменную - которая отвечает за тип транзакций type=='Payments::Transaction::ContractAddTransaction','add','disb'),am=amount*fifelse(z_type=='add',1,-1))][1:20,c(-1,-11,-8,-6)] #убираем ненужные колонкиlk1[disbursement_date<'2020-01-01' & date<='2020-03-01',c(-1,-11,-8,-6)]->lk1glimpse(lk1) #посмотрим на получившийся фаил

$ borrower_id          2, 2, 2, 6, 6, 12, 12, 12, 12, 16, 20, 20, 20, 20, 22, 23, 23, 33, 33, 39, 39, 36, 36$ con_id               1, 1, 1, 2, 2, 4, 4, 4, 4, 5, 7, 7, 7, 7, 8, 9, 9, 10, 10, 11, 11, 12, 12, 12, 12, 13$ disbursement_date    2017-11-23, 2017-11-23, 2017-11-23, 2017-11-24, 2017-11-24, 2017-11-27, 2017-11-27, $ prolongations_count  1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1$ date                 2017-12-15, 2017-11-23, 2017-12-06, 2017-11-24, 2017-12-25, 2017-12-29, 2017-11-27, $ amount               250000, 1000000, 1200000, 2500000, 3500000, 1040000, 1500000, 10000, 1470000, 1500000$ id                   325, 2, 127, 5, 587, 557500, 557499, 557504, 557507, 17, 182865, 182874, 182869, 1828$ z_type               "add", "disb", "add", "disb", "add", "add", "disb", "add", "add", "disb", "disb", "ad$ am                   250000, -1000000, 1200000, -2500000, 3500000, 1040000, -1500000, 10000, 1470000, -150

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

Для этого посмотрим транзакции за все время. И так как цифры крупные, посмотрим сразу в миллионах:

lk1[,.(total_in_mln=sum(am*for_ex)/1e6),.(z_type)]

Видим ясный и понятный результат: транзакции от заемщиков (add) значительно превышают транзакции от выдачи кредита (disb). А ведь часть кредитов наверняка только выдана и по ним даже не подошла дата выплаты.

И здесь же разберемся с темпом, с которым выдавали кредиты:

lk1[!is.na(disbursement_date)&z_type=='disb',.(sum=sum(am*for_ex*-1,na.rm = T)/1e6),.(date=floor_date(disbursement_date,'month',))][,ggplot(.SD,aes(date,sum,label=round(sum,2)))+geom_col(fill=polar_night[2])+ff+tt+  geom_text(aes(y=sum+0.1),col=aurora[1])+  labs(x='месяц выдачи',y='выдача в миллионах долларов',  title='Выдача кредитов в месяц')]

Аналитически видно, что масштабирование продукта в последние несколько месяцев остановилось. Выдаем столько же сколько и выдавали.

Вернемся к Юнит экономике

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

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

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

Удобнее работать с относительным временем, так сможем сравнивать пользователей, пришедших в разное время:

lk1[][,':='( min_date=min(disbursement_date)),.(borrower_id)][,#найдем дату начала работы каждого заемщикаc("gen",'dif'):=.(floor_date(min_date,'1 month'),as.numeric(date-min_date))][, #месяц в котором появился заемщик и растояние в днях между датой его рождения и датой транзакцииn:=uniqueN(borrower_id),][, #количество уникальных клиентов  .(sum=sum(am*for_ex),n=unique(n)),  .(dif)][order(dif)][,":="(bal=bal<-sum/n,cum=cumsum(bal))][,    ggplot(.SD,aes(dif,cum))+ #строим график   geom_line()+   tt+ff+        labs(x='дни жизни заемщика',y='доход в USD',col='поколение',        title='LTV клиента по поколениям и по годам')+   scale_y_continuous(breaks = seq(-200,200,20),   labels =paste0('$',seq(-200,200,20),'k' ))+   scale_x_continuous(breaks = seq(0,1000,20))  ]

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

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

Поэтому стоит посмотреть этот же график, но разбив по поколениям и посмотрев последние 24 поколения. Сразу же наложим сюда стоимость привлечения и разобьем по годам(чисто для удобства сравнения):

lk1[,':='( min_date=min(disbursement_date)),.(borrower_id)][,c("gen",'dif'):=.(floor_date(min_date,'1 month'),as.numeric(date-min_date))][,n:=uniqueN(borrower_id),.(gen)][,  .(sum=sum(am*0.00003),n=unique(n)),  .(gen,dif)# сгруппируем данные  по поколениями][order(gen,dif)][,":="(bal=bal<-sum/n)][,':='(cum=cumsum(bal),m_dif=max(dif)),.(gen)][dif<=m_dif-30 & gen %between% c('2018-01-01','2021-03-31')][, ggplot(.SD,aes(dif,cum,col=factor(gen)))+ geom_line()+ facet_wrap(~factor(year(gen),levels = c(2019,2018)),nrow=2)+tt+ff+ labs(x='дни жизни заемщика',y='доход в USD',col='поколение', title = 'LTV клиента по поколениям и по годам')+ scale_y_continuous(breaks = seq(-200,200,20), labels =paste0('$',seq(-200,200,20),'k' ))+ scale_x_continuous(breaks = seq(0,1000,20))+ geom_hline(yintercept = 695494*for_ex,color='red',size=1)+ geom_hline(yintercept = 0,color='dark red',linetype='dashed')+ geom_text(inherit.aes = F,aes(x=as.Date(600), y=695494*for_ex+3,group=1),label='Стоимость привлечения клиента = $20.86', col='red',size=6)+tt+ff+ theme(legend.text = element_text(size=20),       legend.title = element_text(size=25))+  guides(colour = guide_legend(override.aes = list(size=10)))]

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

1. Точка, из которой начинаются линии - размер первого кредита.

2. Чем выше траектория, тем больше денег на одного пользователя мы заработали в этом поколении.

3. Чем раньше траектория поколения пересекает 0, тем быстрее мы начинаем зарабатывать на привлеченном пользователе.

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

Выводы

1. Продукт становится лучше. Траектории более поздних поколений пересекают ноль и стоимость привлечения раньше, чем траектории старых поколений (150 дней против 180).

2. Разница в темпах окупаемости объясняется скорее стартовой точкой и первой парой кредитов. Например, поколения конца 2018 года окупались лучше, имея меньшую сумму первого кредита и большую сумму второго. Это можно понять по провалу на траектории на 10-40 днях. Так же в 2019 году: поколения с меньшей средней суммой первого кредита окупали себя намного быстрее.

3. На дистанции в один год мы заработаем на одном пользователе от 40 до 60 долларов. Максимальный размер заработка неизвестен. Ни одно поколение пока не вышло на плато, но экстраполируя, оценка в 120-150 долларов на дистанции в 3 года выглядит разумной

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

Мы тут ради ответов, а не картинки смотреть

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

1. Сейчас мы платим 695494. Какая стоимость привлечения для нас является приемлемой? Имеет ли смысл повысить стоимость привлечения на клиента, чтобы получить больший объем?

Ситуация выглядит следующим образом - у нас есть успешный продукт, который мы не масштабируем - и это проблема. Ожидания по доходности от клиента на дистанции год - 40-60 долларов. Текущие затраты на привлечение $20.8 (695494* 0.00003)

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

Если мы можем позволить себе тратить сейчас - чтобы заработать потом - эксперименты оправданы.

2. Насколько здоровой выглядит экономика портфеля и какая динамика здесь и сейчас?

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

Там существенных прорывов пока не было.

3. Мы недавно изменили подход к размеру выдачи и стали в первые кредиты выдавать меньшие чеки. Как это сказалось на продукте?

Вообще - видно что поколения в которых первый кредит был меньшим по размерам окупали себя чуть лучше. Четкий ответ - только через АБ тест.

Итог

Работа продукт менеджера это принимать решения. 2.2 млн строк рассказали о том, что происходит. Такой анализ занимает от 30 минут до пары часов, в зависимости от знания предметной области и грязи в данных. Для такого анализа не нужно ничего, кроме сырых данных и открытого ПО.

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

Из этих же данных несложно собирается еще несколько важных оценок и выводов , но о них в дуговой раз.

Проблематика

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

Какие есть варианты?

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

В чем смысл объединения?

1. Заходить в 5 досок и просматривать отдельно каждого, просто не хватит времени и сил.

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

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

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

Как синхронизировать доски 5 сотрудников в одной доске "дашборде" и при этом не оплачивать лицензию в Трелло или передаточных сервисах типа Zapier?

Решение задачи:

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

Ход решения:

Сначала получим ключи и секретный код для работы с API https://trello.com/app-key (предварительно необходимо авторизоваться под своей учеткой в Trello и открыть себе доски сотрудников с правами администратора)

Сохраняем ключ и токен, который можно найти внизу этой же страницы

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

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

• работы с API Trello trelloR

• работы со временем и временными периодами lubridate

• работы с таблицами и агрегации данных dplyr

Чтобы установить пакеты из основного репозитория CRAN примените базовую функцию install.packages, а для загрузки пакетов с github функцию install_github:

remotes::install_github("jchrom/trelloR")install.packages("lubridate", dependencies = TRUE)install.packages("dplyr ", dependencies = TRUE)

Подключаемся к API и получаем токен:

# Указываем путь к папке куда будет сохранен полученный tokensetwd("C:\\*********\\R_script\\trello")# Получаем tokenmy_token = get_token("my_app", key = "", secret = "",                     expiration = c( "never"))

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

trelloadd <- function(delcard = NULL,                      addcard = NULL,                      nlista = NULL){  # Находим исходный лист сотрудника в Трелло откуда будем брать карточки  для записи в дашборд  ishod_tab <- get_list_cards(addcard)    # Определяем доску куда пишем новые карточки  bid = get_id_board(delcard)    # Получаем все списки с листами на доске дашборда  lid <- get_board_lists(bid)$id[nlista]    # Получаем все карточки в необходимом листе куда пишем  cid<-get_list_cards(lid)    # Удаляем карточки, которые уже есть в листе сотрудника на дашборде  if (length(cid$id)>0) {    for (i in 1:length(cid$id)) {      delete_resource(resource = "card", id = cid$id[i])    }  }else{    print("no-del")  }    # Находим дату создания карточки в исходном листе  dateList<- data.frame(dateadd = NA)  for (i in 1:length(ishod_tab$id)) {    cardID <- ishod_tab$id[i]    dateList[i,1] <- strtoi(strtrim(cardID, 8), 16L)  }  dateList$dateadd <-as.POSIXct(dateList$dateadd, origin = "1970-01-01")    # Циклом собираем лист записи и создаем карточку для переноса в дашборд  if (length(ishod_tab$name)>0) {    for (i in 1:length(ishod_tab$name)) {        payload = list(        idList = lid,        name = ishod_tab$name[i],        desc = paste0(ishod_tab$desc[i],"Date Add: " ,dateList$dateadd[i], "                     В работе уже ", floor(as.vector(difftime(today(),dateList$dateadd[i], units='days'))), " день"),        start = ishod_tab$start [i],        due = ishod_tab$due [i],        pos = "bottom"      )      r <- create_resource("card", body = payload)    }  }else{    print("Ok")  }  if (nrow(bind_rows(ishod_tab$labels[]))>0) {      # Добавляем лейблы (метки)    bid = get_id_board(delcard)    lid <- get_board_lists(bid)$id[nlista]    cid <-get_list_cards(lid)        # Находим карточки без меток    nlab <- which( lapply(ishod_tab$labels, length)!=0 %in% T)    for (i in nlab) {      labl <- ishod_tab$labels[[i]]      for (xi in 1:nrow(labl)) {        r <-  add_label(cid$id[i], color = ishod_tab$labels[[i]][xi,4],                        name = ishod_tab$labels[[i]][xi,3] )        }    }   }else{    print("no_lable")  }}

В функции необходимо будет только заполнять переменные:

• delcard - Это id дашборда в который будет записываться информация из досок сотрудников

• addcard - Это id листа из которого будут браться карточки для переноса в дашборд

• nlista - номер листа в дашборде в который будут заноситься карточки

Получаем delcard

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

Получаем addcard

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

get_board_lists("https://trello.com/b/*****/сотрудник1", query = list(customFieldItems = "true"))

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

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

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

# Сотрудник1 ####trelloadd(delcard = "od*****W",           addcard = "600**********04",          nlista = 2)# Сотрудник2 ####trelloadd(delcard = "od*****W",           addcard = "5fc4********24",          nlista = 3)# Сотрудник3 ####trelloadd(delcard = "od*****W",           addcard = "5e94*********8ce",          nlista = 4)# Сотрудник4 ####trelloadd(delcard = "od*****W",          addcard = "5faa*********c522",          nlista = 5)# Сотрудник5 ####trelloadd(delcard = "od*****W",          addcard = "60744*******3394",          nlista = 6)# Сотрудник6 ####trelloadd(delcard = "od*****W",          addcard = "5e73******b07",          nlista = 7)

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

В итоге:

1. Мы получили полноценную синхронизацию всех досок любого числа сотрудников в trello

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

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

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

Один из подходов к решению данной задачи (чтобы метод кластеризации автоматически отсеивал выбросы) получил название "optimally tuned robust improper maximum likelihood estimator" и был описан вот в этой статье 2017 года (http://dx.doi.org/10.1080/01621459.2015.1100996), а недавно и получил реализацию на R. Поговорим о нем.

Минутка теории

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

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

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

но вот с такими ограничениями

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

Эксперимент

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

library(ggplot2)# Создаем базу данныхset.seed(739)n <- 500 # numer of pointsy <- abs(c(rnorm(n,5,2)))x <- c(rnorm(n,5,1))/sqrt(y)class<-rep(1,n)dat1 <- as.data.frame(cbind(x,y,class)) # Первый кластер - что-то типа кривоугольного прямоугольного треугольникаy <- c(rnorm(n,8,0.4))x <- rlnorm(n, meanlog = log(7), sdlog = 0.2) class<-rep(2,n)dat2 <- as.data.frame(cbind(x,y,class)) # Второй кластер - больше похож на форме на горизонтальный прямоугольникy <- c(runif(n, min = 2, max = 7))x <- c(rnorm(n,8,0.4))class<-rep(3,n)dat3 <- as.data.frame(cbind(x,y,class)) # Третий кластер - вертикальный прямоугольникy <- c(runif(n/10, min = 0, max = 10))x <- c(runif(n/10, min = 0, max = 10))class<-rep(0,n/10)dat4 <- as.data.frame(cbind(x,y,class)) # Шумdat <- rbind(dat1,dat2,dat3,dat4)colnames(dat) <- c("x","y", "class")dat$class <- as.factor(dat$class)dat_test <- as.data.frame(dat)p_base <- ggplot(dat,aes(x=x,y=y,color=class)) + geom_point()ggExtra::ggMarginal(p_base, groupColour = TRUE, groupFill = TRUE)

Получаем вот такую картинку

Пакет otrimle достаточно современен в том плане, что в нем уже есть функция, позволяющая сравнить результаты кластеризации в случае задания нескольких кластеров. Зададим число кластеров от 2 до 5, и попробуем посмотреть на ход выполнения.

library(otrimle)clus <- otrimleg(dat[,c(1,2)], G=2:5, monitor=1) # параметр monitor позволяет видеть ход выполнения

В результате выведется вот что-то такое

В теории, этого уже достаточно для определения оптимального числа кластеров (мы же решаем задачу на максимум, и, в теории, нужно взять число кластеров, для которого iloglik максимален. С другой же стороны, видно, что значимого преимущества в разбиении на 4 или 5 кластеров по сравнению с 3 нет - поэтому выбор количества кластеров в значительной степени эмпирика). Но желательно посмотреть еще и на количество наблюдений внутри кластеров. Поэтому, не отходя далеко, жмем

clus$solution

и получаем следующую информацию

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

clus_2 <-otrimle(dat[,c(1,2)], G = 5, npr.max = 0.01, erc = 20, monitor = TRUE)# npr.max - максимальная доля выбросов в выборке# erc - соотношение максимального/минимального собственных значенияclus_2$codeclus_2$flag

Если clus_2$code возвращает 0, это значит, что задача оптимизации вообще не решилась, если clus_2$code = 1, это означает, что EM-алгоритм не сошелся (и надо поиграться с параметрами), если clus_2$code = 2, то все хорошо, алгоритм все посчитал.

Параметр clus_2$flag дает информацию о внутренней структуре работы алгоритма. Если

clus_2$flag = 1, то присутствует вырождение апостерирорных вероятностей

Фактически, это говорит о наличии лишних кластеров, так как вероятность одного из компонент смеси близка к 0.

clus_2$flag = 2, то для сходимости алгоритма было нарушено предположение о доли шума

clus_2$flag = 3, то для при расчетах ограничение на долю шума было использовано

clus_2$flag = 4, то для при расчетах было использовано ограничение на отношение собственных значений

В нашем случае все хорошо (clus_2$code = 2, clus_2$flag = 4), продолжаем разбираться дальше.

clus_2$mean # центра кластеровhead(clus_2$tau) # вероятности принадлежности к кластерамhead(clus_2$cluster) # принадлежность к кластерам

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

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