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

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

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

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

Представьте, что мы управляем крупным продуктовым магазином и стремимся максимизировать прибыль для каждого продукта, регулируя их цены. Чтобы достичь этой цели, нам необходимо оценить кривую спроса (эластичность) для каждого товара, а затем оптимизировать функцию максимизацию прибыли. Первым шагом в этом процессе является оценка ценовой эластичности спроса, то есть определение того, насколько быстро спрос реагирует на изменение цены на 1%. Для этого у нас есть данные в динамике по i ∈ N продуктам, собранные за t ∈ T периодов. Ценовая эластичность спроса рассчитывается по следующей формуле:

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

Где γc(i), t — набор переменных, которые учитывают средний спрос в каждой категории за определенный период времени, а δi — набор переменных, которые отражают изменения спроса, не связанные со временем, для каждого отдельного продукта. Такая форма «фиксированного эффекта» является стандартной и широко используется во многих регрессионных моделях для учета ненаблюдаемых факторов, которые могут оказывать влияние. Эта (объединенная) регрессионная модель позволяет нам определить среднюю эластичность β для всех N товарных единиц. Это означает, что магазин может установить средний уровень цен на все товары в своем ассортименте, чтобы максимизировать прибыль:

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

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

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

Суть байесовских иерархических моделей

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

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

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

• Индивидуальные параметры, которые применяются к каждому отдельному объекту.

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

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

Давайте рассмотрим это на примере ценовой эластичности. В идеале, мы хотим получить ценовую эластичность на уровне товарной единицы. Мы оцениваем:

Где:

• βi ∼ Normal(βc(i),σi)

• βc(i) ∼ Normal(βg,σc(i))

• βg ∼ Normal(μ,σ)

Единственное отличие от первого уравнения заключается в том, что мы заменяем глобальный член β на индивидуальные продуктовые βi. Уточним, что эластичность на уровне товарных единиц βi определяется на основе нормального распределения, которое центрируется вокруг среднего значения эластичности на уровне категории βc(i), которое, в свою очередь, рассчитывается на основе общей глобальной эластичности βg для всех групп. Для дисперсий σ мы также можем предположить иерархическую структуру, но здесь мы просто установим для них стандартные априорные значения, чтобы не усложнять пример. В данном случае мы предполагаем, что априорные значения: μ =-2, σ=1, σc (i)=1, σi=1. Эта априорная формулировка предполагает, что глобальная эластичность является эластичной, 95% эластичностей находятся в диапазоне от -4 до 0, со среднеквадратическим отклонением 1 на каждом иерархическом уровне. Чтобы проверить, правильно ли мы определили эти априорные значения, мы бы провели априорные прогностические проверки (которые не рассматриваются в этой статье), чтобы увидеть, могут ли наши априорные предположения восстановить данные, которые мы наблюдаем.

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

Реализация

В этом разделе мы подробно рассмотрим процесс реализации описанной ранее модели с использованием пакета Numpyro на Python. Numpyro представляет собой легковесный вероятностный язык программирования, основанный на JAX, что обеспечивает autograd‑ и JIT‑компиляцию на графических GPU / TPU / CPU. Мы начнем с генерации наших синтетических данных, затем определим модель и подгоним ее под эти данные. В заключение мы визуализируем результаты нашей работы.

Процесс генерации данных

Мы моделируем данные о продажах, где спрос находится в лог‑линейной зависимости от цены, а эластичность на уровне продукта генерируется на основе гауссова распределения βi ∼ Normal (-2,0,7). В каждый период времени с вероятностью 50% мы добавляем случайное изменение цены, временные тренды для конкретной категории и случайный шум. Эти факторы мультипликативно суммируются для формирования нашего логарифмического ожидаемого спроса. Чтобы получить фактический спрос, мы экспонируем логарифмический ожидаемый спрос и получаем количество реализованных единиц, проданных согласно распределению Пуассона. Для повышения точности оценок мы отбираем только те единицы, которые были проданы более 100 раз. В результате получаем N = 11 798 товаров за T = 156 периодов (еженедельные данные за 3 года). Исходя из этого набора данных, истинная глобальная эластичность составляет βg=-1,6, при этом эластичности на уровне категории варьируются в пределах βc(i)∈[-1.68,-1.48].

Обратите внимание, что в этом DGP (динамическая гамма‑пуассоновская модель — методе, который мы используем для моделирования) не учитываются многие сложные аспекты реального мира. Мы не рассматриваем факторы, которые могут одновременно влиять как на цены, так и на спрос, такие как рекламные акции. Также мы не учитываем факторы, которые могут мешать работе модели. Этот пример предназначен только для того, чтобы продемонстрировать, как мы можем восстановить эластичность для конкретного продукта в рамках модели, разработанной Уэллсом. Он не ставит своей целью объяснить, как правильно определить, что цена является независимой переменной. Тем не менее, я рекомендую заинтересовавшимся читателям полистать книгу Causal Inference for the Brave and True для более глубокого понимания методов причинно‑следственного анализа.

import numpy as np
import pandas as pd
def generate_price_elasticity_data(N: int = 1000,
                                   C: int = 10,
                                   T: int = 50,
                                   price_change_prob: float = 0.2,
                                   seed = 42) -> pd.DataFrame:

    """

    Генерируем синтетические данные для анализа ценовой эластичности спроса.
    Данные генерируются с помощью
    """
    
    if seed is not None:
        np.random.seed(seed)

    
    # Спроса и тренды категории
    category_base_demand = np.random.uniform(1000, 10000, C)
    category_time_trends = np.random.uniform(0, 0.01, C)
    category_volatility = np.random.uniform(0.01, 0.05, C)  # Случайная волатильность для каждой категории
    category_demand_paths = np.zeros((C, T))
    category_demand_paths[:, 0] = 1.0
    shocks = np.random.normal(0, 1, (C, T-1)) * category_volatility[:, np.newaxis]
    trends = category_time_trends[:, np.newaxis] * np.ones((C, T-1))
    cumulative_effects = np.cumsum(trends + shocks, axis=1)
    category_demand_paths[:, 1:] = category_demand_paths[:, 0:1] + cumulative_effects

    # Эффекты продукта
    product_categories = np.random.randint(0, C, N)
    product_a = np.random.normal(-2, .7, size=N)
    product_a = np.clip(product_a, -5, -.1)

 
    # Начальные цены на каждый продукт
    initial_prices = np.random.uniform(100, 1000, N)
    prices = np.zeros((N, T))
    prices[:, 0] = initial_prices

   
    # Генерация случайных значений и определение факта изменения цен
    random_values = np.random.rand(N, T-1)
    change_mask = random_values < price_change_prob

    
    # Генерация коэффициентов изменения (от -20% до +20%)
    change_factors = 1 + np.random.uniform(-0.2, 0.2, size=(N, T-1))

    
    # Создание матрицы для хранения множителей
    multipliers = np.ones((N, T-1))

    
    # Применение коэффициентов изменения только там, где должны произойти изменения 
    multipliers[change_mask] = change_factors[change_mask]

    
    # Применение изменений кумулятивно для распространения цен
    for t in range(1, T):
        prices[:, t] = prices[:, t-1] * multipliers[:, t-1]

    
    # Генерация множители для конкретного продукта
    product_multipliers = np.random.lognormal(3, 0.5, size=N)
    
    # Получение временных эффектов для каждой категории продукта (размерность: N x T)
    time_effects = category_demand_paths[product_categories][:, np.newaxis, :].squeeze(1)

    
    # Чтобы временные эффекты не становились отрицательными
    time_effects = np.maximum(0.1, time_effects)

    
    # Генерация периодического шума для всех продуктов и периодов времени
    period_noise = 1 + np.random.uniform(-0.05, 0.05, size=(N, T))

    
    # Получение базового спроса категории для каждого продукта
    category_base = category_base_demand[product_categories]

    
    # Расчет базового спроса
    base_demand = (category_base[:, np.newaxis] *
                   product_multipliers[:, np.newaxis] *
                   time_effects *
                   period_noise)

    # Логарифм спроса
    alpha_ijt = np.log(base_demand)

    # Логарифм цены
    log_prices = np.log(prices)

    # Логарифм ожидаемого спроса
    log_lambda = alpha_ijt + product_a[:, np.newaxis] * log_prices

    # Обратное преобразование из логарифма для получения параметров интенсивности
    lambda_vals = np.exp(log_lambda)

    # Генерация проданных товарных единиц
    units_sold = np.random.poisson(lambda_vals)  # Формат: (N, T)

    
    # Создание массивов индексов для всех комбинаций продуктов и периодов времени
    product_indices, time_indices = np.meshgrid(np.arange(N), np.arange(T), indexing='ij')
    product_indices = product_indices.flatten()
    time_indices = time_indices.flatten()

    
    # Получение категорий для всех продуктов
    categories = product_categories[product_indices]

    
    # Получение всех цен и проданных единиц
    all_prices = np.round(prices[product_indices, time_indices], 2)
    all_units_sold = units_sold[product_indices, time_indices]

    
    # Вычисление эластичности
    product_elasticity = product_a[product_indices]
    df = pd.DataFrame({
        'product': product_indices,
        'category': categories,
        'time_period': time_indices,
        'price': all_prices,
        'units_sold': all_units_sold,
        'product_elasticity': product_elasticity
    })

    
    return df

# Сохраняем только единицы с >X продаж
def filter_dataframe(df, min_units = 100):
    temp = df[['product','units_sold']].groupby('product').sum().reset_index()
    unit_filter = temp[temp.units_sold>min_units]['product'].unique()
    filtered_df = df[df['product'].isin(unit_filter)].copy()
    
    # Сводка
    original_product_count = df['product'].nunique()
    remaining_product_count = filtered_df['product'].nunique()
    filtered_out = original_product_count - remaining_product_count

    
    print(f"Filtering summary:")
    print(f"- Original number of products: {original_product_count}")
    print(f"- Products with > {min_units} units: {remaining_product_count}")
    print(f"- Products filtered out: {filtered_out} ({filtered_out/original_product_count:.1%})")

    # Глобальная эластичность и эластичность по категориям
    global_elasticity = filtered_df['product_elasticity'].mean()
    filtered_df['global_elasticity'] = global_elasticity

    # Эластичность по категориям
    category_elasticities = filtered_df.groupby('category')['product_elasticity'].mean().reset_index()
    category_elasticities.columns = ['category', 'category_elasticity']
    filtered_df = filtered_df.merge(category_elasticities, on='category', how='left')

    # Сводка
    print(f"\nElasticity Information:")
    print(f"- Global elasticity: {global_elasticity:.3f}")
    print(f"- Category elasticities range: {category_elasticities['category_elasticity'].min():.3f} to {category_elasticities['category_elasticity'].max():.3f}")

    
    return filtered_df

df = generate_price_elasticity_data(N = 20000, T = 156, price_change_prob=.5, seed=42)
df = filter_dataframe(df)
df.loc[:,'cat_by_time'] = df['category'].astype(str) + '-' + df['time_period'].astype(str)
df.head()

Сводка по фильтрации:

• Первоначальное количество продуктов: 20 000

• Количество с более чем 100 единицами: 11 798

• Количество отфильтрованных продуктов: 8202 (41,0%)

Информация об эластичности:

• Глобальная эластичность: -1,598

• Диапазон эластичности по категориям: от -1,681 до -1,482

Модель

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

Для определения групп параметров в модели используется таблица основных параметров NumPyro:

• global_a: 1 глобальный параметр ценовой эластичности с априорным значением Normal(−2,1).

• category_a: C=10 эластичностей на уровне категории с априорными значениями, центрированными по глобальному параметру, и среднеквадратическим отклонением, равным 1.

• product_a: N = 11 798 эластичностей для конкретного продукта с априорными значениями, центрированными по параметрам соответствующей категории, и среднеквадратическим отклонением, равным 1.

• product_effect: N=11 798 базовых эффектов спроса на конкретный продукт со среднеквадратическим отклонением 3.

• time_cat_effects: (T=156)⋅(C=10) эффектов, изменяющихся во времени, для каждой комбинации категории и времени со среднеквадратическим отклонением 3.

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

import jax
import jax.numpy as jnp
import numpyro
import numpyro.distributions as dist
from numpyro.infer.reparam import LocScaleReparam
def model(df: pd.DataFrame, outcome: None):
  
    # Определение индексов
    product_idx, unique_product = pd.factorize(df['product'])
    cat_idx, unique_category = pd.factorize(df['category'])
    time_cat_idx, unique_time_cat = pd.factorize(df['cat_by_time'])

    # Преобразование рядов цен и товарных единиц в массивы jax
    log_price = jnp.log(df.price.values)
    outcome = jnp.array(outcome) if outcome is not None else None

    # Генерация сопоставления
    product_to_category = jnp.array(pd.DataFrame({'product': product_idx, 'category': cat_idx}).drop_duplicates().category.values, dtype=np.int16)

    # Создание таблиц для хранения параметров
    category_plate = numpyro.plate("category", unique_category.shape[0])
    time_cat_plate = numpyro.plate("time_cat", unique_time_cat.shape[0])
    product_plate = numpyro.plate("product", unique_product.shape[0])
    data_plate = numpyro.plate("data", size=outcome.shape[0])

    # ОПРЕДЕЛЕНИЕ ПАРАМЕТРОВ МОДЕЛИ
    global_a = numpyro.sample("global_a", dist.Normal(-2, 1), infer={"reparam": LocScaleReparam()})

    with category_plate:
        category_a = numpyro.sample("category_a", dist.Normal(global_a, 1), infer={"reparam": LocScaleReparam()})

    with product_plate:
        product_a = numpyro.sample("product_a", dist.Normal(category_a[product_to_category], 2), infer={"reparam": LocScaleReparam()})
        product_effect = numpyro.sample("product_effect", dist.Normal(0, 3))

    with time_cat_plate:
        time_cat_effects = numpyro.sample("time_cat_effects", dist.Normal(0, 3))

    # Расчет ожидаемого спроса
    def calculate_demand():
        log_demand = product_a[product_idx]*log_price + time_cat_effects[time_cat_idx] + product_effect[product_idx]
        expected_demand = jnp.exp(jnp.clip(log_demand, -4, 20)) # Ограничение для стабильности 
        return expected_demand

    demand = calculate_demand()

    with data_plate:
        numpyro.sample(
            "obs",
            dist.Poisson(demand),
            obs=outcome

        )

   numpyro.render_model(
    model=model,
    model_kwargs={"df": df,"outcome": df['units_sold']},
    render_distributions=True,
    render_params=True,
)

Оценка

Хотя существует множество способов оценки данного уравнения, в этом конкретном случае мы используем стохастический вариационный вывод (SVI). SVI — это метод оптимизации на основе градиента, позволяющий минимизировать расхождение Кульбака‑Лейблера между предлагаемым апостериорным распределением и истинным апостериорным распределением путем минимизации нижней вариационной границы (ELBO). Этот метод отличается от метода Монте‑Карло с марковскими цепями (MCMC), который производит выборку непосредственно из истинного апостериорного распределения. В реальных приложениях SVI более эффективен и легко масштабируется до больших наборов данных. Для этого приложения мы устанавливаем случайное начальное значение, инициализируем guide (семейство апостериорного распределения, предполагается, что это диагональная нормаль), определяем скорость обучения и оптимизатор в Optax и запускаем оптимизацию для 1 000 000 (занимает ~ 1 час) итераций. Хотя модель, возможно, сходилась заранее, потери продолжают улучшаться на незначительную величину даже после выполнения оптимизации в течение 1 000 000 итераций. Наконец, мы отражаем потери (логарифмические) на графике.

from numpyro.infer import SVI, Trace_ELBO, autoguide, init_to_sample
import optax
import matplotlib.pyplot as plt

rng_key = jax.random.PRNGKey(42)
guide = autoguide.AutoNormal(model, init_loc_fn=init_to_sample)
# Определение скорости обучения
learning_rate_schedule = optax.exponential_decay(
    init_value=0.01,
    transition_steps=1000,
    decay_rate=0.99,
    staircase = False,
    end_value = 1e-5,
)

# Определение оптимизатора
optimizer = optax.adamw(learning_rate=learning_rate_schedule)
svi = SVI(model, guide, optimizer, loss=Trace_ELBO(num_particles=8, vectorize_particles = True))

# Запуск SVI
svi_result = svi.run(rng_key, 1_000_000, df, df['units_sold'])
plt.semilogy(svi_result.losses);

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

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

predictive = numpyro.infer.Predictive(
    autoguide.AutoNormal(model, init_loc_fn=init_to_sample),
    params=svi_result.params,
    num_samples=1000
)

samples = predictive(rng_key, df, df['units_sold'])

# Извлечение мат ожидания и среднеквадратического отклонения
results = {}
excluded_keys = ['product_effect', 'time_cat_effects']
for k, v in samples.items():
    if k not in excluded_keys:
        results[f"{k}"] = np.mean(v, axis=0)
        results[f"{k}_std"] = np.std(v, axis=0)

# Оценки эластичности продуктов
prod_elasticity_df = pd.DataFrame({
    'product': df['product'].unique(),
    'product_elasticity_svi': results['product_a'],
    'product_elasticity_svi_std': results['product_a_std'],
})
result_df = df.merge(prod_elasticity_df, on='product', how='left')

# Оценки эластичности по категориям
prod_elasticity_df = pd.DataFrame({
    'category': df['category'].unique(),
    'category_elasticity_svi': results['category_a'],
    'category_elasticity_svi_std': results['category_a_std'],
})
result_df = result_df.merge(prod_elasticity_df, on='category', how='left')

# Оценки глобальной эластичности
result_df['global_a_svi'] = results['global_a']
result_df['global_a_svi_std'] = results['global_a_std']
result_df.head()

Результаты

В представленном коде показаны истинные и рассчитанные значения эластичности для каждого продукта. Каждая точка ранжирована в соответствии с их истинным значением эластичности, которое обозначено черным цветом. Также отображается расчетная эластичность, полученная с помощью модели. Можно заметить, что оцененная эластичность следует траектории истинной эластичности со средней абсолютной погрешностью, равной примерно 0.0724. Точки, окрашенные в красный цвет, представляют продукты, 95% доверительного интервала которых не включают истинную эластичность, в то время как синие точки обозначают продукты, 95% доверительного интервала которых совпадают с истинной эластичностью. Учитывая, что глобальное мат ожидание составляет -1.598, это означает, что средняя погрешность на уровне продукта составляет около 4.5%. Можно отметить, что оценки SVI точно следуют за схемой истинных эластичностей, хотя и с некоторым уровнем шума, особенно когда эластичности становятся более отрицательными. На верхней правой панели мы отображаем график зависимости между погрешностью расчетных значений эластичности и истинными значениями. С увеличением отрицательности истинной эластичности, наша модель начинает давать все более неточные результаты.

Чтобы получить оценки эластичности на уровне категории и на глобальном уровне, мы можем воспользоваться двумя методами. Первый — это увеличение эластичности на уровне каждого продукта в рамках категории. Второй метод — получение оценок на уровне категории непосредственно из апостериорных параметров. Если мы обратимся к оценкам эластичности на уровне категории, представленным в левом нижнем углу, то увидим, что как оценки, восстановленные из модели, так и оценки, полученные на основе загруженных выборок эластичности на уровне продукта, слегка смещены в сторону нуля, со значением средней абсолютной ошибки около 0.033. Однако доверительный интервал, задаваемый параметром уровня категории, охватывает истинный параметр, в отличие от оценок, полученных с использованием начальной загрузки. Это свидетельствует о том, что при определении эластичности на уровне группы следует напрямую использовать параметры на уровне группы, вместо того чтобы прибегать к начальной загрузке более детализированных оценок. При анализе на глобальном уровне оба метода дают точную оценку параметра в пределах 95% доверительного интервала. Однако глобальный параметр оказывается выше, чем первоначальная загрузка на уровне продукта, из‑за увеличения среднеквадратических погрешностей.

Соображения

1. Байесовский иерархический метод недооценивает апостериорную дисперсию: Одним из недостатков использования SVI для оценки является то, что он недооценивает апостериорную дисперсию. Хотя мы подробно рассмотрим эту тему в следкющей статье, целевая функция для SVI учитывает только разницу в математическом ожидании нашего предполагаемого распределения и истинного распределения. Это означает, что он не учитывает полную структуру корреляции между параметрами в апостериоре. Аппроксимация среднего поля, обычно используемая в SVI, предполагает условную (при построении предыдущей иерархии) независимость между параметрами, которая игнорирует любые ковариации между продуктами в пределах одной и той же иерархии. Это означает, что при наличии каких‑либо побочных эффектов (таких как каннибализация или перекрестная ценовая эластичность) они не будут учитываться в пределах доверительных интервалов. Из‑за этого допущения о среднем поле оценки неопределенности обычно слишком надежны. Это приводит к слишком узким доверительным интервалам, которые не могут точно отразить истинные значения параметров с необходимой интенсивностью. На верхнем левом рисунке мы видим, что только 9.7% эластичностей на уровне продукта соответствуют их истинным значениям. В следующем статье мы поговорим о возможных решениях этой проблемы.

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

def elasticity_plots(result_df, results=None):
    # Создаем фигуру с сеткой 2х2
    fig = plt.figure(figsize=(12, 10))
    gs = fig.add_gridspec(2, 2)

 
    # Эластичность продукта
    ax1 = fig.add_subplot(gs[0, 0])

    
    # Подготовка данных
    df_product = result_df[['product','product_elasticity','product_elasticity_svi','product_elasticity_svi_std']].drop_duplicates()
    df_product['product_elasticity_svi_lb'] = df_product['product_elasticity_svi'] - 1.96*df_product['product_elasticity_svi_std']
    df_product['product_elasticity_svi_ub'] = df_product['product_elasticity_svi'] + 1.96*df_product['product_elasticity_svi_std']
    df_product = df_product.sort_values('product_elasticity')
    mae_product = np.mean(np.abs(df_product.product_elasticity-df_product.product_elasticity_svi))
    colors = []
    for i, row in df_product.iterrows():
        if (row['product_elasticity'] >= row['product_elasticity_svi_lb'] and 
            row['product_elasticity'] <= row['product_elasticity_svi_ub']):
            colors.append('blue')  # В пределах доверительного интервала
        else:
            colors.append('red')   # За пределами доверительного интервала

   
    # Процентное количество точек в пределах допустимых значений
    within_bounds_pct = colors.count('blue') / len(colors) * 100

    
    # Данные графика
    ax1.scatter(range(len(df_product)), df_product['product_elasticity'], 
                color='black', label='True Elasticity', s=20, zorder=3)

  
    ax1.scatter(range(len(df_product)), df_product['product_elasticity_svi'], 
                color=colors, label=f'SVI Estimate (MAE: {mae_product:.4f}, Coverage: {within_bounds_pct:.1f}%)', 
                s=3, zorder=2)
    ax1.set_xlabel('Product Index (sorted by true elasticity)')
    ax1.set_ylabel('Elasticity Value')
    ax1.set_title('SVI Estimates of True Product Elasticities')
    ax1.legend()
    ax1.grid(alpha=0.3)

   
    # Взаимосвязь между средней абсолютной ошибкой и истинной эластичностью
    ax2 = fig.add_subplot(gs[0, 1])

    
    # Расчет средней абсолютной ошибки для каждого продукта
    temp = result_df[['product','product_elasticity', 'product_elasticity_svi']].drop_duplicates().copy()
    temp['product_error'] = temp['product_elasticity'] - temp['product_elasticity_svi']
    temp['product_mae'] = np.abs(temp['product_error'])
    correlation = temp[['product_mae', 'product_elasticity']].corr()

   
    # Данные графика
    ax2.scatter(temp['product_elasticity'], temp['product_error'], alpha=0.5, s=5, color = colors)
    ax2.set_xlabel('True Elasticity')
    ax2.set_ylabel('Error (True - Estimated)')
    ax2.set_title('Relationship Between True Elasticity and Estimation Accuracy')
    ax2.grid(alpha=0.3)
    ax2.text(0.5, 0.95, f"Correlation: {correlation.iloc[0,1]:.3f}", 
                transform=ax2.transAxes, ha='center', va='top',
                bbox=dict(boxstyle='round', facecolor='white', alpha=0.7))

    
    # Эластичность категория 
    ax3 = fig.add_subplot(gs[1, 0])

    
    # Уникальные категории и эластичности
    category_data = result_df[['category', 'category_elasticity', 'category_elasticity_svi', 'category_elasticity_svi_std']].drop_duplicates()
    category_data = category_data.sort_values('category_elasticity')

    
    # Загрузка мат ожиданий и эластичностей продуктов в каждой категории
    bootstrap_means = []
    bootstrap_ci_lower = []
    bootstrap_ci_upper = []

   
    for cat in category_data['category']:
        # Получение эластичности продукта для этой категории
        prod_elasticities = result_df[result_df['category'] == cat]['product_elasticity_svi'].unique()

        
        # Загрузка мат ожиданий
        boot_means = [np.mean(np.random.choice(prod_elasticities, size=len(prod_elasticities), replace=True)) 
                     for _ in range(1000)]

       
        bootstrap_means.append(np.mean(boot_means))
        bootstrap_ci_lower.append(np.percentile(boot_means, 2.5))
        bootstrap_ci_upper.append(np.percentile(boot_means, 97.5))

    
    category_data['bootstrap_mean'] = bootstrap_means
    category_data['bootstrap_ci_lower'] = bootstrap_ci_lower
    category_data['bootstrap_ci_upper'] = bootstrap_ci_upper

    
    # Вычисление средней абсолютной ошибки
    mae_category_svi = np.mean(np.abs(category_data['category_elasticity_svi'] - category_data['category_elasticity']))
    mae_bootstrap = np.mean(np.abs(category_data['bootstrap_mean'] - category_data['category_elasticity']))

        
    # Вывод данных на график
    left_offset = -0.2
    right_offset = 0.2
    x_range = range(len(category_data))
    ax3.scatter(x_range, category_data['category_elasticity'], 
                color='black', label='True Elasticity', s=50, zorder=3)

    
    # Загрузка эластичности продукта
    ax3.scatter([x + left_offset for x in x_range], category_data['bootstrap_mean'], 
                color='green', label=f'Bootstrapped Product Estimate (MAE: {mae_bootstrap:.4f})', s=30, zorder=2)
    for i in x_range:
        ax3.plot([i + left_offset, i + left_offset], 
                [category_data['bootstrap_ci_lower'].iloc[i], category_data['bootstrap_ci_upper'].iloc[i]], 
                color='green', alpha=0.3, zorder=1)

    
    # Оценки SVI на уровне категорий
    ax3.scatter([x + right_offset for x in x_range], category_data['category_elasticity_svi'], 
                color='blue', label=f'Category SVI Estimate (MAE: {mae_category_svi:.4f})', s=30, zorder=2)
    for i in x_range:
        ci_lower = category_data['category_elasticity_svi'].iloc[i] - 1.96 * category_data['category_elasticity_svi_std'].iloc[i]
        ci_upper = category_data['category_elasticity_svi'].iloc[i] + 1.96 * category_data['category_elasticity_svi_std'].iloc[i]
        ax3.plot([i + right_offset, i + right_offset], [ci_lower, ci_upper], color='blue', alpha=0.3, zorder=1)

    ax3.set_xlabel('Category Index (sorted by true elasticity)')
    ax3.set_ylabel('Elasticity')
    ax3.set_title('Comparison with True Category Elasticity')
    ax3.legend()
    ax3.grid(alpha=0.3)

    # Глобальная эластичность
    ax4 = fig.add_subplot(gs[1, 1])
    temp = result_df[['product','product_elasticity_svi','global_elasticity']].drop_duplicates()
    bootstrap_means = [np.mean(np.random.choice(np.array(temp['product_elasticity_svi']), 100)) for i in range(10000)]
    global_means = np.random.normal(result_df['global_a_svi'].iloc[0], result_df['global_a_svi_std'].iloc[0], 10000)
    true_global = np.unique(temp.global_elasticity)[0]
    p_mae = np.abs(np.mean(bootstrap_means) - true_global)
    g_mae = np.abs(np.mean(global_means) - true_global)

   
    # Данные графика
    ax4.hist(bootstrap_means, bins=100, alpha=0.3, density=True, 
             label=f'Product Elasticity Bootstrap (MAE: {p_mae:.4f})')
    ax4.hist(global_means, bins=100, alpha=0.3, density=True, 
             label=f'Global Elasticity Distribution (MAE: {g_mae:.4f})')
    ax4.axvline(x=true_global, color='black', linestyle='--', linewidth=2, 
                label=f'Global Elasticity: {true_global:.3f}', zorder=0)
    ax4.set_xlabel('Elasticity')
    ax4.set_ylabel('Frequency')
    ax4.set_title('Global Elasticity Comparison')
    ax4.legend()
    ax4.grid(True, linestyle='--', alpha=0.7)

    # Отображение
    plt.tight_layout()
    plt.show()

elasticity_plots(result_df)

Заключение

Альтернативные варианты использования: Помимо оценки ценовой эластичности спроса, байесовские иерархические модели также находят множество других применений в Data Science. В розничной торговле они могут прогнозировать спрос на товары в уже существующих магазинах и решать проблему «холодного запуска» новых торговых точек. Для этого модели заимствуют информацию из уже открытых магазинов или сетей, которые группируются в рамках единой иерархии. В рекомендательных системах модели HB позволяют оценить предпочтения пользователей на основе сочетания характеристик как самих пользователей, так и товаров. Эта структура позволяет давать релевантные рекомендации новым пользователям, основываясь на их поведении в группах, и постепенно переходить к индивидуальным рекомендациям по мере накопления истории покупок. Если же нет возможности сформировать группы по когортам, то для группировки схожих единиц на основе их характеристик может быть использован алгоритм K‑средних.

В заключение, эти модели также могут быть полезны для объединения результатов экспериментальных и наблюдательных исследований. Ученые могут использовать исторические оценки увеличения объема наблюдений и дополнить их данными недавно проведенных A/B‑тестов, чтобы сократить необходимый объем выборки для экспериментов, используя уже имеющиеся знания. Такой подход позволяет создать основу для непрерывного обучения, где каждый новый эксперимент строится на предыдущих результатах, а не начинается с нуля. Для команд, сталкивающихся с ограниченными ресурсами, это означает более быстрое время до понимания (особенно в сочетании с суррогатными моделями) и более эффективные конвейеры экспериментов.

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

Тем не менее, этот обзор должен стать практическим руководством для внедрения иерархического байесовского метода в ваш инструментарий. Эти модели предоставляют основу для обработки сложных многоуровневых структур данных, которые часто встречаются в реальных бизнес‑задачах. Когда вы начнете применять эти подходы, я с нетерпением жду возможности услышать о вашем опыте, проблемах, успехах и новых возможностях использования этого класса моделей. Пожалуйста, не стесняйтесь присылать вопросы, идеи или примеры по электронной почте или через LinkedIn. Если у вас есть отзывы об этой статье или вы хотите обсудить другую тему в разделе причинно‑следственного вывода и машинного обучения, не стесняйтесь обращаться. Спасибо за внимание!

Углубляетесь в ML и ищете что‑то мощнее стандартных моделей?

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

Перед уроком — короткий тест: проверьте, насколько вы уже готовы к следующему шагу в ML и «потяните» ли программу продвинутого курса по ML.