4 способа количественно оценить толстые хвосты с помощью Python

4 метода измерить размер густых хвостов с использованием Python

Интуиция и пример кода

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

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

Толстый (кошачий) хвост. Изображение из Canva.

В первой статье этой серии мы представили понятие толстых хвостов, которое описывает степень, в которой редкие события определяют общую статистику распределения. Мы увидели крайний пример толстых хвостов с использованием распределения Парето, где, например, 80% продаж генерируются 20% клиентов (и 50% продаж генерируются всего 1% клиентов).

Парето, Степенные законы и Толстые хвосты

Что не объясняют в курсах по статистике

towardsdatascience.com

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

Спектр Толстоты хвостов. Изображение автора.

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

4 способа количественной оценки толстых хвостов

Хотя нет «истинной» меры толстоты хвостов, существует несколько эвристик (т.е. эмпирических правил), которые мы можем использовать на практике для количественной оценки толстоты хвостов данных. Здесь мы рассмотрим 4 таких эвристики. Мы начнем с представления каждой техники концептуально, а затем перейдем к примеру кода на Python.

Эвристика 1: Показатель Толстоты Хвоста Степенного Закона

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

Пример распределений Степенного Закона с различными значениями α. Изображение автора.

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

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

Эвристика 2: Коэффициент эксцесса (то есть негауссовость)

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

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

Определение коэффициента эксцесса согласно ref [1] и [2]. Изображение автора.

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

Эта мера хорошо работает, когда все моменты являются конечными [3]. Однако одно из основных ограничений заключается в том, что коэффициент эксцесса не определен для некоторых распределений, например, Парето с α <= 4, что делает его бесполезным для большинства данных с толстыми хвостами.

Эвристика 3: Лог-нормальное распределение σ

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

Функция плотности вероятности лог-нормального распределения [4]. Изображение автора.

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

Мы можем получить эту меру похожим образом, как и Эвристика 1. А именно, мы подгоняем лог-нормальное распределение к нашим данным и извлекаем значение σ подгонки. Хотя это простая процедура, она (как и Эвристика 1) терпит неудачу, когда лог-нормальная подгонка плохо объясняет исходные данные.

Эвристика 4: К

Предшествующие эвристики (H) начинались с определенного распределения в виду (то есть H1 – Степенной закон и H3 – Лог-нормальное распределение). Однако на практике наши данные редко точно следуют конкретному распределению.

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

Это стимулирует использование мер неспецифичных для распределения. Одна такая мера была предложена Талебом в ref [3]. Предложенная метрика (κ) определяется для унимодальных данных со средним значением и принимает значения между 0 и 1, где 0 означает, что данные являются максимально тонкими хвостами, а 1 означает, что данные являются максимально толстыми хвостами. Он определяется согласно выражению ниже.

Определение метрики κ Талеба [3]. Изображение автора.

Метрика сравнивает две выборки (скажем, Sₙ₀ и Sₙ), где Sₙ является суммой n выборок, взятых из определенного распределения. Например, если мы оцениваем нормальное распределение и выбираем n=100, мы бы взяли 100 выборок из нормального распределения и сложили их вместе, чтобы получить S₁₀₀.

M(n) в выражении выше обозначает среднее абсолютное отклонение, определенное по следующему уравнению. Эта мера разброса вокруг среднего значения обычно более надежна, чем стандартное отклонение [3] [5].

Определение среднего абсолютного отклонения по уравнению κ [3]. Изображение автора.

Для упрощения можно выбрать n₀=1, что дает нам нижеследующее выражение.

κ с n₀=1. Изображение автора.

Ключевым термином здесь является M(n)/M(1), где M(n) характеризует разброс вокруг среднего значения для суммы n выборок (некоторого распределения).

Для распределений с узкими хвостами, M(30) будет относительно близким к M(1), поскольку данные обычно находятся близко к среднему значению. Таким образом, M(30)/M(1) ~ 1.

Однако для данных с толстыми хвостами M(30) будет значительно больше, чем M(1). Таким образом, M(30)/M(1) >> 1. Это иллюстрируется ниже, где левый график показывает, как шкалируется разброс для суммы гауссовых распределений, а правый график показывает, как он шкалируется для Парето.

Шкалирование M(n) и M(1) для гауссовых (слева) и Парето 80-20 (справа) распределений. Обратите внимание на метки оси y. Примечание: шкала разброса гаусса увеличивается из-за суммирования n распределений. Изображение автора.

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

Если все это оказалось для вас слишком математическим, вот главное: Большой κ = толстый хвост, маленький κ = узкий хвост.

Пример кода: Количественная оценка толстопузости (распределения хвоста) данных реальных социальных медиа

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

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

YouTube-Blog/power-laws/3-quantifying-fat-tails на main · ShawhinT/YouTube-Blog

Коды для дополнения видеороликов YouTube и блогов на VoAGI. – YouTube-Blog/power-laws/3-quantifying-fat-tails на main…

github.com

Начнем с импорта полезных библиотек.

import matplotlib.pyplot as pltimport pandas as pdimport numpy as npimport powerlawfrom scipy.stats import kurtosis

Затем мы загрузим каждый набор данных и сохраним их в словаре.

filename_list = ['ПодписчикиVoAGI', 'ЗаработокYT', 'ПоказыLI']df_dict = {}for filename in filename_list:    df = pd.read_csv('data/'+filename+'.csv')    df = df.set_index(df.columns[0]) # установить индекс    df_dict[filename] = df

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

for filename in filename_list:    df = df_dict[filename]        # построить гистограммы (функция ниже определена в блокноте на GitHub)    plot_histograms(df.iloc[:,0][df.iloc[:,0]>0], filename, filename.split('-')[1])    plt.savefig("images/"+filename+"_histograms.png")        # печать топ-5 записей    print("Топ-5 записей по процентам")    print((df.iloc[:,0]/df.iloc[:,0].sum()).sort_values(ascending=False)[:5])    print("")
Гистограммы для ежемесячных подписчиков VoAGI. Изображение автора.
Гистограммы для заработка на YouTube видео. Примечание: если вы заметите разницу с предыдущей статьей, это потому, что я обнаружил проблемную запись в данных (поэтому просмотр полезен 😅). Изображение автора.
Гистограммы для ежедневных показов на LinkedIn. Изображение автора.

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

Топ-5 записей по процентам для каждого набора данных. Изображение автора.

Согласно этому представлению, подписчики VoAGI являются наиболее “толстым” хвостом, приходящемуся на 60% от подписчиков только за 2 месяца. Доходы от YouTube также имеют ярко выраженный “толстый” хвост, где около 60% доходов приходится на всего 4 видео. Показы на LinkedIn, кажется, имеют меньше “толстого” хвоста.

Хотя мы можем более качественно понять “толстый” хвост, просто посмотрев на данные, давайте сделаем это более количественным с помощью наших 4 эвристик.

Эвристика 1: Индекс “толстого” хвоста степенного закона

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

for filename in filename_list:    df = df_dict[filename]    # выполнить аппроксимацию степенным законом    results = powerlaw.Fit(df.iloc[:,0])    # выводить результаты    print("")    print(filename)    print("-"*len(filename))    print("Аппроксимация степенным законом")    print("a = " + str(results.power_law.alpha-1))    print("xmin = " + str(results.power_law.xmin))    print("")
Результаты аппроксимации степенным законом. Изображение авторства.

Результаты выше соответствуют нашей качественной оценке, что подписчики VoAGI имеют наиболее широкий и хвостатый распределение, за ними следуют заработки на YouTube и охваты на LinkedIn (помните, меньшее α – это более широкий хвост).

Эвристика 2: Коэффициент эксцесса

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

for filename in filename_list:    df = df_dict[filename]    # выводить результаты    print(filename)    print("-"*len(filename))    print("коэффициент эксцесса = " + str(kurtosis(df.iloc[:,0], fisher=True)))    print("")
Значения коэффициента эксцесса для каждого набора данных. Изображение авторства.

Коэффициент эксцесса рассказывает другую историю, чем Эвристика 1. Ранжировка хвостатости по этому показателю следующая: LinkedIn > VoAGI > YouTube.

Однако следует относиться к этим результатам с некоторым скепсисом. Как мы видели с результатами аппроксимации степенным законом выше, все 3 набора данных аппроксимируются степенным законом с α < 4, что означает, что коэффициент эксцесса бесконечен. Таким образом, хотя вычисление возвращает значение, разумно быть подозрительным к этим числам.

Эвристика 3: σ для логнормального распределения

Мы также можем использовать библиотеку powerlaw, чтобы получить оценки σ, аналогичные Эвристике 1. Вот как это выглядит.

for filename in filename_list:    df = df_dict[filename]    # выполнить аппроксимацию степенным законом    results = powerlaw.Fit(df.iloc[:,0])    # выводить результаты    print("")    print(filename)    print("-"*len(filename))    print("Аппроксимация логнормальным распределением")    print("mu = " + str(results.lognormal.mu))    print("sigma = " + str(results.lognormal.sigma))    print("")
Результаты аппроксимации логнормальным распределением. Изображение авторства.

Изучая значения σ выше, мы видим, что все аппроксимации подразумевают наличие широкого хвоста в данных, где оценки σ для подписчиков VoAGI и охватов на LinkedIn примерно одинаковы. Однако у заработков на YouTube значительно большее значение σ, что свидетельствует о более широком хвосте (гораздо).

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

Эвристика 4: κ по Талебу

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

def среднее_абсолютное_отклонение(S):    """    Вычисление среднего абсолютного отклонения от входной выборки S    """    M = np.mean(np.abs(S - np.mean(S)))    return Mdef генерация_n_выборки(X,n):    """    Функция для генерации n случайных выборок размера len(X) из массива X    """    # инициализировать выборку    S_n=0        for i in range(n):        # случайно выбираем len(X) наблюдений из X и добавляем их в выборку        S_n = S_n + X[np.random.randint(len(X), size=int(np.round(len(X))))]        return S_ndef каппа(X,n):    """    Метрика каппа Талеба от n0=1, описанная здесь: https://arxiv.org/abs/1802.05495        Примечание: K_1n = kappa(1,n) = 2 - ((log(n)-log(1))/log(M_n/M_1)), где M_n обозначает среднее абсолютное отклонение суммы n случайных выборок    """    S_1 = X    S_n = generate_n_sample(X,n)        M_1 = mean_abs_deviation(S_1)    M_n = mean_abs_deviation(S_n)        K_1n = 2 - (np.log(n)/np.log(M_n/M_1))    return K_1n

Первая функция mean_abs_deviation() вычисляет среднее абсолютное отклонение, как определено выше.

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

Наконец, я объединяю mean_abs_deviation(S) и generate_n_sample(X,n) для реализации рассчета κ, определенного ранее, и вычисляю его для каждого набора данных.

n = 100 # количество выборок для включения в расчет каппыfor filename in filename_list:    df = df_dict[filename]    # печать результатов    print(filename)    print("-"*len(filename))    print("kappa_1n = " + str(kappa(df.iloc[:,0].to_numpy(), n)))    print("")
Значения κ(1,100) для каждого набора данных. Изображение автора.

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

Следовательно, я запускаю тот же расчет 1000 раз и печатаю среднюю каппу(1,100) для каждого набора данных.

num_runs = 1_000kappa_dict = {}for filename in filename_list:    df = df_dict[filename]    kappa_list = []    for i in range(num_runs):        kappa_list.append(kappa(df.iloc[:,0].to_numpy(), n))        kappa_dict[filename] = np.array(kappa_list)    print(filename)    print("-"*len(filename))    print("mean kappa_1n = " + str(np.mean(kappa_dict[filename])))    print("")
Средние значения κ(1,100) из 1000 запусков для каждого набора данных. Изображение автора.

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

Примечание: можно сравнить эти значения с Таблицей III в [3], чтобы лучше понять каждое значение каппы. А именно, эти значения сравнимы с распределением Парето с α от 2 до 3.

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

Заключение

Хотя могло бы показаться соблазнительным размечивать данные как “толстохвостые” (или нет), толстохвостость существует на спектре. Здесь мы разбирали 4 эвристики для количественной оценки толстохвостости данных.

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

👉 Больше о степенных законах и толстых хвостах: Введение | Подгонка степенному закону

Ресурсы

Связь: Мой веб-сайт | Запишите звонок | Спросите меня что-нибудь

Социальные сети: YouTube 🎥 | LinkedIn | Twitter

Поддержка: Приготовьте мне кофе ☕️

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

Сообщество для предпринимателей в сфере данных. 👉 Присоединяйтесь к Discord!

VoAGI.com

[1] Scipy Kurtosis

[2] Scipy Moment

[3] arXiv:1802.05495 [stat.ME]

[4] https://en.wikipedia.org/wiki/Log-normal_distribution

[5] Pham-Gia, T., & Hung, T. (2001). The mean and median absolute deviations. Mathematical and Computer Modelling, 34(7–8), 921–936. https://doi.org/10.1016/S0895-7177(01)00109-1