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

Постановка проблемы

Сегодня в эпоху цифровой трансформации A/B тестирование стало неотъемлемой частью стратегии принятия решений, превратившись из инновационного инструмента в стандарт индустрии. Напомню, что мы работаем в дирекции по развитию программы лояльности Х5, и зачастую проверяем гипотезы об увеличении наших бизнес-метрик с помощью балльных механик (например, «купи шоколад Х, получи Y баллов программы лояльности»). Тестировать дополнительные балльные механики на всем ассортименте магазина не есть хорошо для бизнеса, поэтому иногда мы сталкиваемся с проблемой выбора пилотных механик.

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

Дополнительный кэшбэк X% на покупку фруктов и овощей
Дополнительный кэшбэк X% на швейцарский шоколад
Дополнительный кэшбэк X% на хлебобулочные изделия собственного производства

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

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

Описание подхода

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

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

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

Хорошая новость, такая метрика есть – и это JS-дивергенция (Jensen-Shannon divergence).

Определение дивергенции

Формально JS-дивергенция является модернизацией метода KL-дивергенции (Kullback–Leibler divergence). Дивергенция KL измеряет разницу между двумя распределениями вероятностей. Она также известна как относительная энтропия. Строгое определение такое: KL дивергенция измеряет количество дополнительных битов, необходимых для кодирования выборок из одного распределения при использовании кода на основе другого распределения. Другими словами, это объём информации, потерянной, когда одно распределение используется для аппроксимации другого. При этом KL дивергенция обозначается как D(P||Q), где P и Q – сравниваемые распределения. При этом она обладает рядом следующих свойств:

Нессиметрична, $D_{KL} (P||Q) \ne D_{KL}(Q||P).$
Не имеет верхней границы и вообще интерпретируемых границ.

Покажем интуитивно, как работает KL дивергенция. Представим, что мы имеем две столбчатые диаграммы P и Q, в которых представлены по оси Y вероятности какого-то события А, В и С (рис. 1). Допустим, со временем изменились вероятности наступления каждого из событий.

Рис.1 Столбчатая диаграмма вероятностей в разные временные периоды

Попробуем оценить изменение размера бинов. В качестве первой метрики для исследования возьмем отношение бинов. Тогда для события А получится А , для зеленого бина В для синего А. Если взять среднее арифметическое всех трех значений, получим . Такая метрика обла��ает большим недостатком: на среднее арифметическое сильно влияет самое большое значение отношений бинов, а также важно, от какого распределения бины стоят в числителе, а какие в знаменателе – от этого конечный результат конечной метрики будет другим, отсюда вытекает несимметричность метрики. Для устранения этого недостатка сильного влияния больших чисел на среднее каждую дробь оборачивают в функцию логарифма. Формально KL дивергенция считает отношение бинов в степени логарифма, который дополнительно умножается на вес каждого бина в распределениях. Тогда итоговая формула для KL дивергенции:

Для дискретной записи:

$D_ {KL}(P||Q)=\sum P(x) * log \frac{P(x)} {Q(x)}$

Для непрерывных величин:

$D_ {KL}(P||Q)=\int P(x) * log \frac{P(x)} {Q(x)}$

Если мы сравнивали бы не P c Q, а Q с P:

$D_ {KL}(Q||P)=\int Q(x) * log \frac{Q(x)} {P(x)}, т.е. D_{KL}(P||Q) \ne D_{KL}(Q||P)$

Отличие JS дивергенции от KL дивергенции в том, что эта метрика сравнивает каждое из распределений P и Q не друг с другом, а со средневзвешенным распределением вероятности каждого из P и Q распределений:

$D_{JS}(P||Q)=\frac {1}{2}D_{KL}(P||q)+\frac{1}{2}D_{KL}(Q||q), q=\frac{P+Q}{2}$

Тогда итоговая формула для JS дивергенции:

$D_{JS}(P||Q)=\frac{1}{2} \int P(x) * log \frac{P(x)} {\frac{P(x)+Q(x)}{2}} +\frac{1}{2} \int Q(x) * log \frac{Q(x)}{\frac{P(x)+Q(x)}{2}} =D_{JS}(Q||P)$

Из последней формулы видно, что метрика становится симметричной. Математически можно показать, что максимальное значение такой метрики при основании логарифма 2 будет равно 1, а минимальное равно 0. Действительно, если мы имеем идентичные распределения, то P=Q, тогда:

$\frac{1}{2} \int P(x) * log_2 \frac{P(x)}{\frac{P(x)+Q(x)}{2}} =\frac{1}{2} \int P(x) * log_2 \frac{P(x)}{\frac{P(x)+P(x)}{2}}=\frac{1}{2} \int P(x) * 0=0$

Теперь определим максимальное значение. Так как и , то:

$\frac{1}{2} \int P(x) * log_2 \frac{P(x)}{\frac{P(x)+Q(x)}{2}} =\frac{1}{2} \int Px * (1 + log_2\frac{P(x)}{P(x)+Q(x)} )$ $\frac{1}{2}\int P(x) * (1 + log_2\frac{P(x)}{P(x)+Q(x)})≤ \frac{1}{2}\int P(x) ≤ \frac{1}{2}$

Так как имеем две аналогичные половинки в формуле JS-дивергенции, то максимальное значение JS дивергенции не будет превышать 1, когда . Таким образом, мы получили следующее свойство метрики:

JS дивергенция симметрична при сравнении распределений.
Метрика обладает интерпретируемыми границами [0,1].

Если JS дивергенция JS = 1, то распределения максимально отличаются между собой, и пилотная группа максимально специфична относительно генеральной совокупности. Если JS = 0, то мы имеем дело с идентичными распределениями. Кроме этого, можно заметить, что JS дивергенция отлично справляется с нулевыми значениями.

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

Практическая симуляция JS дивергенции

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

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

$bins=\max(len(data1), len(data2)) ^ \frac{2}{3},$

Код для вычисления KL дивергенции:

import numpy as np
from scipy.stats import entropy

def KL(data1, data2):
    """
    Расчет KL-дивергенции между двумя наборами данных.
    
    Parameters:
    - data1: Первый набор данных.
    - data2: Второй набор данных.
    
    Returns:
    - kl_divergence_data1_data2: KL-дивергенция от data1 к data2.
    - kl_divergence_data2_data1: KL-дивергенция от data2 к data1.
    """
    
    # Определение количества бинов согласно рекомендации 
    bins = max(len(data1), len(data2))
    bins = int(bins ** (2 / 3))
    
    # Вычисление гистограмм
    hist1, bin_edges = np.histogram(data1, bins=bins, density=True)
    hist2, _ = np.histogram(data2, bins=bin_edges, density=True)  # Используем одинаковые интервалы
    
    # Нормализация гистограмм
    p = hist1 / hist1.sum()
    q = hist2 / hist2.sum()
    
    # Добавляем небольшое значение для избежания деления на ноль
    p = np.clip(p, 1e-10, None)
    q = np.clip(q, 1e-10, None)
    
    # Расчет KL-дивергенции вручную
    kl_divergence_data1_data2 = round(np.sum(p * np.log(p / q)), 2)
    kl_divergence_data2_data1 = round(np.sum(q * np.log(q / p)), 2)
    
    return kl_divergence_data1_data2, kl_divergence_data2_data1

Код для вычисления JS дивергенции:

import numpy as np
from scipy.stats import entropy

def JS(data1, data2):
    """
    Расчет дивергенции Йенсена — Шеннона (Jensen-Shannon divergence) между двумя наборами данных.
    
    Parameters:
    - data1: Первый набор данных.
    - data2: Второй набор данных.
    
    Returns:
    - js_divergence: Дивергенция Йенсена — Шеннона.
    """
    
    # Определение числа бинов для гистограмм
    bins = max(len(data1), len(data2))
    bins = int(bins ** (2 / 3))
    
    # Создание гистограмм
    hist1, bin_edges = np.histogram(data1, bins=bins, density=True)
    hist2, _ = np.histogram(data2, bins=bin_edges, density=True)  # Используем одинаковые интервалы для бинов
    
    # Нормализация
    p = hist1 / hist1.sum()
    q = hist2 / hist2.sum()
    
    # Добавляем небольшое значение для избежания деления на ноль
    p = np.clip(p, 1e-10, None)
    q = np.clip(q, 1e-10, None)
    
    # Расчет среднего распределения
    m = 0.5 * (p + q)
    
    # Используем формулу для расчета дивергенции Йенсена — Шеннона
    js_divergence = 0.5 * (entropy(p, m, base=2) + entropy(q, m, base=2))
    
    # Округление результата до двух знаков после запятой
    js_divergence = round(js_divergence, 2)

    return js_divergence

# Пример использования
data1 = np.random.normal(0, 1, 1000)
data2 = np.random.normal(1, 1, 1000)

js_divergence = JS(data1, data2)
print(f'JS divergence: {js_divergence}')

Код для отрисовки графиков

def plot_divergences(klab_list, klba_list, js_list):
    """
    Функция строит три scatter-графика для списков дивергенций:
    - klab_list: KL-дивергенция A->B
    - klba_list: KL-дивергенция B->A
    - js_list: Jensen-Shannon дивергенция
    
    Добавляет горизонтальные линии 50-го и 95-го квантилей.
    """
    # Вычисление квантилей
    klab_q95 = np.quantile(klab_list, 0.95)
    klba_q95 = np.quantile(klba_list, 0.95)
    js_q95 = np.quantile(js_list, 0.95)

    klab_q50 = np.quantile(klab_list, 0.5)
    klba_q50 = np.quantile(klba_list, 0.5)
    js_q50 = np.quantile(js_list, 0.5)

    n_points = len(klab_list)
    x_vals = range(n_points)

    plt.figure(figsize=(10, 7))

    plt.subplot(3, 1, 1)
    sns.scatterplot(x=x_vals, y=klab_list, label='kl_ab')
    plt.title('KL-divergence (A to B)')
    plt.axhline(klab_q95, color='r', linestyle='--', label='95th Quantile')
    plt.axhline(klab_q50, color='b', linestyle='--', label='50th Quantile')
    plt.legend()

    plt.subplot(3, 1, 2)
    sns.scatterplot(x=x_vals, y=klba_list, label='kl_ba')
    plt.title('KL-divergence (B to A)')
    plt.axhline(klba_q95, color='r', linestyle='--', label='95th Quantile')
    plt.axhline(klba_q50, color='b', linestyle='--', label='50th Quantile')
    plt.legend()

    plt.subplot(3, 1, 3)
    sns.scatterplot(x=x_vals, y=js_list, label='JS-divergence')
    plt.title('Jensen-Shannon Divergence')
    plt.axhline(js_q95, color='r', linestyle='--', label='95th Quantile')
    plt.axhline(js_q50, color='b', linestyle='--', label='50th Quantile')
    plt.legend()

    plt.tight_layout()
    plt.show()

А/А тесты

В качестве А/А теста проведем измерения на двух нормальных распределениях с средним 0, дисперсией 1, количеством наблюдений 2000. Зафиксируем одно из распределений, а вто��ое нормальное распределение всегда будет генерироваться новое, но с теми же параметрами среднего и дисперсии (анимация 1). Сделаем 10 000 таких генераций и на каждом этапе будем считать три метрики: KL дивергенцию (напомним, что метрика несимметричная относительно сравниваемых распределений) и JS дивергенцию. Результаты представлены на рисунке 2.

*Анимация 1. Нормальное распределение. А/А тесты*

Код для отрисовки значений дивергенции

# Генерация данных
a = np.random.normal(0, 1, 2000)

klab_list = []
klba_list = []
js_list = []

# Вычисление KL и JS дивергенций для 10000 итераций
for k in tqdm(range(10000)):
    b = np.random.normal(0, 1, 2000)
    kl_ab, kl_ba = KL(a, b)
    js = JS(a, b)
    klab_list.append(kl_ab)
    klba_list.append(kl_ba)
    js_list.append(js)

# Отрисовка графиков
plot_divergences(klab_list, klba_list, js_list)

*Рисунок 2. Метрики дивергенции на А/А тестах*

По значениям видно, что 95 перцентиль всех значений JS дивергенции имеет значение 0,1. KL дивергенции – около 0,95 и 0,4 при сравнении распределений. Значения не изменяются, что позволяет сделать вывод о том, что метрика на АА тестах дает одинаковые значения. В идеальном случае мы бы хотели получить значение метрики, близкое к нулю, но анимация выше как раз показывает, что генерируемые распределения не всегда обладают полной площадью перекрытия, и предложенная метрика на это реагирует.

А/B тесты. Увеличение среднего.

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

*Анимация 2. А/B тесты, изменение среднего*

Код для отрисовки значений A/B теста по увеличению среднего

# Генерация данных
a = np.random.normal(0, 1, 2000)

klab_list = []
klba_list = []
js_list = []

# Вычисление KL и JS дивергенций для 10 000 итераций
for k in tqdm(range(10000)):
    b = np.random.normal(0 + k / 1000, 1, 2000)
    
    kl_ab, kl_ba = KL(a, b)
    js = JS(a, b)
    
    klab_list.append(kl_ab)
    klba_list.append(kl_ba)
    js_list.append(js)

# Отрисовка графиков
plot_divergences(klab_list, klba_list, js_list)

*Рисунок 3. Метрики дивергенции на А/В тестах, изменение среднего*

Видно, что 95-й перцентиль для JS дивергенции при увеличивающейся разнице не меняется, остается на уровне 0,1. KL дивергенция остается на уровне 0,5 и 0,6 соответственно. Таким образом, мы детектируем первый подводный камень метрики – она не чувствительна к изменениям центральных мер тенденции между распределениями.

A/B тесты. Влияние выбросов.

Посмотрим, как реагируют метрики на выбросы. По мере увеличения итерации расчета будем добавлять в выборку выбросы в левую и правую часть, от чего хвосты одного распределения будут становиться тяжелее, чем исходное (анимация 4). Сделаем 10 тыс. итераций, в конце расчета будет 100 выбросов в группе из 2000 наблюдений. Результаты расчета метрик представлены на рисунке 4.

*Анимация 3. А/B тесты, добавление выбросов*

Код для отрисовки значений A/B теста по влиянию выбросов

# Генерация данных
a = np.random.normal(0, 1, 2000)

klab_list = []
klba_list = []
js_list = []
for k in tqdm(range(10000)):

    b = np.random.normal(0, 1, 2000-int(k/100))
    
    b = np.append(b, np.random.normal(10, 1, int(k/200) ))  # Добавляем 5 выбросов в b
    b = np.append(b, np.random.normal(-10, 1, int(k/200) ))  # Добавляем 5 выбросов в b    
    
    kl_ab, kl_ba = KL(a, b)
    js = JS(a, b)
    klab_list.append(kl_ab)
    klba_list.append(kl_ba)
    js_list.append(js)

plot_divergences(klab_list, klba_list, js_list)

*Рисунок 4. Метрики дивергенции на А/В тестах, добавление выбросов*

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

А/B тесты. Сравнение двух разных форм распределений.

Теперь сравним два абсолютно разных по форме распределения. Нормальное и эскпоненциальное. При этом параметры средних значений и дисперсии у экспоненциального распределения не будут отличаться. То есть проведем А/А тест на двух разных по форме распределений (анимация 4). Результаты расчета метрик представлены на рисунке 5.

*Анимация 4. А/B тесты, сравнение двух разных по форме распределений*

Код для отрисовки значений A/B теста при сравнении распределений

# Генерация нормального распределения
a = np.random.normal(0, 1, 2000)

# Списки для хранения значений
klab_list = []
klba_list = []
js_list = []

# Основной цикл для вычисления KL-дивергенций и JS-дивергенции
for k in tqdm(range(10000)):
    # Генерация экспоненциального распределения
    b = np.random.exponential(1, 2000)
    
    # Вычисление KL-дивергенций и JS-диверсии
    kl_ab, kl_ba = KL(a, b)  # Предполагается, что функция KL определена
    js = JS(a, b)            # Предполагается, что функция JS определена
    
    # Добавление значений в списки
    klab_list.append(kl_ab)
    klba_list.append(kl_ba)
    js_list.append(js)

# Отрисовка графиков
plot_divergences(klab_list, klba_list, js_list)

*Рисунок 5. Метрики дивергенции на А/В тестах, сравнение нормального и экспоненциального распределений*

При сравнении двух совершенно разных по форме распределений можно сразу увидеть, как отреагировала метрика JS дивергенции, 95-й перцентиль около 0,8. При этом 95-й перцентиль значений KL дивергенции находится на уровне 12 и 10 соответственно. Зная интерпретируемые границы JS дивергенции можно заключить, что значение 0,8 достаточно близко к 1, значит по этой метрике можно сделать вывод, что распределения абсолютно разные. В отличие от KL дивергенции, которая несимметрична, и не совсем понятны её границы.

Заключение

Мы получили практические знания о метриках дивергенции:

KL/JS дивергенция не чувствительна к изменениям центральных мер тенденции между распределениями.
KL/JS дивергенция чувствительна к выбросам.
KL/JS дивергенция детектирует разницу между разными классами распределений.
JS дивергенция имеет строгие интерпретируемые границы и лежит в интервале [0,1].

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

Ну и напоследок раскроем имена авторов данного материала: Никита Сумнительный, Егор Карнаух

Дивергенция как мера специфичности аудитории вашей пилотной механики