andrew_brdk Feb 26 at 12:20

Байесовская оценка А/Б-тестов

43 min

1.9K

Математика*Big Data*Веб-аналитика*Аналитика мобильных приложений*Статистика в IT

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

- А/Б тесты
- Байесовское моделирование
- Конверсии
- Средние
- Выручка на пользователя
- Заказы на посетителя
- Заключение
- Ссылки

Репозиторий: https://github.com/andrewbrdk/Bayesian-AB-Testing .

Библиотеки

import numpy as np
import pandas as pd
import scipy.stats as stats
import plotly.graph_objects as go

from collections import namedtuple

np.random.seed(7)

А/Б тесты

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

Увеличение стоимости (справа) может привести к росту выручки, но падению конверсии. Эффект непредсказуем и требует измерения.

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

При оценке эффекта сильные изменения заметны в динамике метрики (см. падение на 30%), но слабые изменения могут быть незаметны на фоне колебаний (см. рост на 3%). Кроме того, изменение метрик после релиза не всегда связано с новой функциональностью.

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

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

Причинная диаграмма [CausalDAG] А/Б-тестов следующая. Метрики определяются действиями пользователей в сервисе. Действия зависят от версии сайта или приложения (например, доступных тарифных планов), внешних факторов (например, сезонности) и также отличаются между сегментами пользователей (новые, постоянные клиенты и др.). В А/Б-тесте версии запускают одновременно и пользователей случайно делят между вариантами. Внешние факторы сохраняются, но при сравнении за одинаковый период их влияние на группы одинаково. При случайном делении пользователей состав сегментов можно считать одинаковым. В итоге разница метрик между группами объясняется функциональностью приложения.

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

По итогам эксперимента нужно оценить метрики, эффект и выбрать "лучшую" группу. Точные значения метрик неизвестны. Их удобнее рассматривать как случайные величины. Распределения вероятностей этих величин нужно подобрать для наибольшей совместимости с экспериментальными данными. Сравнение распределений позволяет оценить эффект. Для презентации удобна точечная оценка метрик и интервал наибольшей плотности вероятности. Например, среднее значение метрики в группе А $p_A = 7.1 \pm 0.2$ , в группе Б $p_B = 7.4 \pm 0.3$ . Эффект можно охарактеризовать относительной разностью $(p_B - p_A) / p_A = 4.2 \pm 0.2 \%$ . Для выбора "лучшей" группы оценить с какой вероятностью метрика в группе Б больше метрики в группе А, например, $P(p_B > p_A) = 95\%$ . Вероятность здесь и далее понимается в субъективном смысле - как мера уверенности в определенном исходе процесса с несколькими возможными исходами [SubjProb].

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

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

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

Для оценки распределений метрик на основе экспериментальных данных используется байесовское моделирование [SR, SGBS].

Байесовское моделирование

Первый пример: если утром облачно, будет ли днем дождь? Для ответа можно посчитать отношение кол-ва облачных дней с дождем к общему кол-ву облачных дней $P(\mbox{Дождь | Облачно}) = (\mbox{Облачно, дождь})/(\mbox{Облачно})$ . Общее количество облачных дней складывается из облачных дней с дождем и облачных дней без дождя $(\mbox{Облачно}) = \mbox{(Облачно, дождь) + (Облачно, без дождя)}$ . Пусть дождливых дней в году $P(\mbox{Дождь}) = 20\%$ , вероятность утренней облачности в дождливый день $P(\mbox{Облачно | Дождь}) = 70\%$ , в день без дождя $P(\mbox{Облачно | Без дождя}) = 10\%$ . Количество облачных дней с дождем можно выразить через соответствующие вероятности $\mbox{(Облачно, дождь)} = (\mbox{Всего дней}) P(\mbox{Дождь}) P(\mbox{Облачно | Дождь})$ , аналогично для количества дней без дождя. После подстановки $P(\mbox{Дождь | Облачно}) = (0.7 \cdot 0.2)/(0.7 \cdot 0.2 + 0.1 \cdot 0.8) = 63.6 \%$ .

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

$\begin{split}& P(\mbox{Дождь | Облачно}) = \frac{\mbox{Облачно, дождь}}{\mbox{Облачно, дождь + Облачно, без дождя}} \\\\& = \frac{P(\mbox{Облачно | дождь})P(\mbox{Дождь})}{P(\mbox{Облачно | Дождь})P(\mbox{Дождь}) + P(\mbox{Облачно | Без дождя})P(\mbox{Без дождя})}\\\\& = \frac{0.7 \cdot 0.2}{0.7 \cdot 0.2 + 0.1 \cdot 0.8} = 63.6 \%\end{split}$

В оценке вероятности дождя при облачности $P(\mbox{Дождь | Облачно})$ важно учитывать долю дождливых дней $P(\mbox{Дождь})$ и вероятность облачности в день без дождя $P(\mbox{Облачно | Без дождя})$ . Их игнорирование ведет к ошибке базового процента [BaseFal].

Еще один пример. Вечером в парке вы видите непонятный предмет. Издалека видны только очертания. Вы пытаетесь угадать, что это. Формально можно использовать соотношение Байеса. Нужно предположить возможные варианты - листья, упавшая шапка, птица, лужа. Для каждого варианта оценить вероятность наблюдаемых очертаний $P(\mbox{Очертания | Листья})$ , $P(\mbox{Очертания | Шапка})$ и т.д. Также учесть распространенность вариантов - листья встречаются чаще потярянной шапки $P(\mbox{Листья}) > P(\mbox{Шапка})$ . С помощью соотношения Байеса определить вероятность каждого из вариантов $P(\mbox{Листья | Очертания}) \propto P(\mbox{Очертания | Листья})P(\mbox{Листья})$ и выбрать наиболее подходящий. Когда вы подходите ближе, предмет недовольно оглядывается и быстро залезает на дерево - это оказалась белка, вы давно не видели их в парке.

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

На этом примере видны основные элементы байесовского моделирования. Есть данные/наблюдения/факты $\mathcal{D}$ и объяснения/гипотезы/модели/предположения $\mathcal{H}_i$ . Для выбора одной гипотезы нужно оценить, насколько каждая из гипотез согласуется с данными. Для этого вычисляют вероятности получить данные в рамках каждой модели - функции правдоподобия $P(\mathcal{D}|\mathcal{H}_i)$ . Прошлый опыт учитывают уверенностью в модели - априорных вероятностях $P(\mathcal{H}_i)$ . По соотношению Байеса вычисляют обновленную уверенность с учетом данных - апостериорные вероятности $P(\mathcal{H}_i|\mathcal{D})$ . На их основе выбирают наиболее подходящую модель. Модели необходимо валидировать - хотя для объяснения данных удается выбрать лучшую из предложенных гипотез, ни одна из гипотез может не соответствовать реальности.

$\begin{split}P(\mathcal{H}_i | \mathcal{D}) &= \frac{ P(\mathcal{D} | \mathcal{H}_i) P(\mathcal{H}_i) }{P(\mathcal{D})}= \frac{ P(\mathcal{D} | \mathcal{H}_i) P(\mathcal{H}_i) }{\sum \limits_i P(\mathcal{D} | \mathcal{H}_i) P(\mathcal{H}_i) }\\P(\mathcal{H}_i | \mathcal{D}) &\mbox{ - апостериорное распределение вероятности} \\P(\mathcal{D} | \mathcal{H}_i) &\mbox{ - функция правдоподобия}\\P(\mathcal{H}_i) &\mbox{ - априорное распределение вероятности}\\P(\mathcal{D}) &\mbox{ - нормировочный коэффициент}\end{split}$

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

Следующий пример. На страницу сайта зашло N = 1000 человек, n_s = 100 из них нажали кнопку "Продолжить". Как выглядит распределение возможных значений конверсии ? Вероятность конверсии каждого пользователя можно считать одинаковой, все возможные априорные значения равновероятными.

Необходимо для заданных n_s и оценить вероятность $P(\mathcal{H} | \mathcal{D}) = P(p | n_s, N)$ . По соотношению Байеса $P(p | n_s, N) \propto P(n_s, N | p) P(p)$ . Пользователь делает клик с вероятностью и не делает с вероятностью 1-p . Клики пользователей можно моделировать последовательностью случайных величин с 2 исходами - схемой Бернулли [BernProc, SciPyBern]. Вероятность n_s конверсий из при шансе на успех задается биномиальным распределением $P(\mathcal{D} | \mathcal{H}) = P(n_s, N | p) = \mbox{Binom}(n_s, N; p)$ [BinomDist, SciPyBinom]. Т.к. все возможные априорные значения конверсий равновероятны, априорное распределение равномерно $P(\mathcal{H}) = P(p) = \mbox{Unif}(0, 1) = 1$ . Апостериорное распределение будет бета-распределением [BetaDist, SciPyBeta].

$\begin{split}P(\mathcal{D} | \mathcal{H}) = P(n_s, N | p) & = \mbox{Binom}(n_s, N; p) = C^{n_s}_{N} p^{n_s} (1 - p)^{N-n_s}\\\\P(\mathcal{H}) = P(p) & = \mbox{Unif}(0, 1) = 1\\\\P(\mathcal{H} | \mathcal{D}) = P(p | n_s, N) & = \frac{P(n_s, N | p) P(p)}{P(n_s, N)}= \frac{P(n_s, N | p) P(p)}{\int_0^1 d p P(n_s, N | p) P(p)}\\& = \frac{p^{n_s} (1 - p)^{N-n_s}}{\int_0^1 d p (1 - p)^{N-n_s} p^{n_s} }= \mbox{Beta}(p; n_s + 1, N - n_s + 1)\end{split}$

График апостериорного распределения P(p | n_s, N) ниже. Мода совпадает со средним в выборке n_s/N , наиболее вероятные значения лежат вблизи этого значения.

Оценка конверсии

ns = 100
ntotal = 1000

p_samp = ns / ntotal
p_dist = stats.beta(a=ns+1, b=ntotal-ns+1)

xaxis_max = 0.2
x = np.linspace(0, xaxis_max, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=p_dist.pdf(x), line_color='black', name='Распределение'))
fig.add_trace(go.Scatter(x=[p_samp, p_samp], y=[0, max(p_dist.pdf(x))], 
                         line_color='black', mode='lines', line_dash='dash', name='Среднее в выборке'))
fig.update_layout(title='$\mbox{Апостериорное распределение} \, P(p | n_s, N)$',
                  xaxis_title='$p$',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[0, xaxis_max],
                  hovermode="x",
                  height=500)
fig.show()

Плотность вероятности конверсий задается бета-распределением. Мода совпадает со средним в выборке.

Еще один пример. На версию А страницы веб-сайта зашло 1000 человек, 100 нажали кнопку "Продолжить". На версию Б зашло 1000 человек, 110 нажали кнопку продолжить. С какой вероятностью конверсия страницы Б выше страницы А?

Нужна вероятность P(p_B > p_A) . Апостериорное распределение конверсий каждой группы вычисляется как в предыдущем примере $P(p; n_s, N) = \mbox{Beta}(p; n_s + 1, N - n_s + 1)$ . Вероятность можно оценить аналитически, посчитав P(p_B - p_A > 0) [ProbConv], но чаще это делают численно. Для этого генерируют выборки из апострериорных распределений и сравнивают их между собой. Для заданных параметров $P(p_B > p_A) = 76.7 \%$ .

Сравнение конверсий

na = 1000
sa = 100
nb = 1000
sb = 110

npost = 50000

p_dist_a = stats.beta(a=sa+1, b=na-sa+1)
p_dist_b = stats.beta(a=sb+1, b=nb-sb+1)

xaxis_max = 0.2
x = np.linspace(0, xaxis_max, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=p_dist_a.pdf(x), line_color='black', name='А'))
fig.add_trace(go.Scatter(x=x, y=p_dist_b.pdf(x), line_color='black', opacity=0.3, name='Б'))
fig.update_layout(title='Апостериорные распределения',
                  xaxis_title='$p$',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[0, xaxis_max],
                  hovermode="x",
                  height=500)
fig.show()

samp_a = p_dist_a.rvs(size=npost)
samp_b = p_dist_b.rvs(size=npost)
p_b_gt_a = np.sum(samp_b > samp_a) / npost

print(f"P(p_b > p_a): {p_b_gt_a}")

Апостериорные распределения конверсий в обеих группах задаются бета-распределениями. Конверсия группы Б выше А с вероятностью 76.7%.

Конверсии

В А/Б-тестах часто сравнивают конверсии в заказ, клик по кнопке и другие действия. Если пользователь не влияет на других, для моделирования можно использовать процесс Бернулли. При конверсии вероятность, что n_s человек из совершит целевое действие, задается биномиальным распределением $P(n_s, N | p) = \mbox{Binom}(n_s, N | p)$ . Априорное распределение конверсий удобно задать бета-распределением $P(p) = \mbox{Beta}(p; \alpha, \beta)$ . Зависимость бета-распределения от без учета нормировочных коэффициентов $\mbox{Beta}(p; \alpha, \beta) \propto p^{\alpha-1}(1-p)^{\beta-1}$ . Эта же форма сохранится для произведения правдоподобия на априорное распределение $P(p | n_s, N) \propto \mbox{Binom}(p, N) \mbox{Beta}(p; \alpha, \beta) \propto p^{n_s + \alpha - 1} (1-p)^{N - n_s + \beta - 1}$ . Важна только зависимость от , остальные множители войдут в нормировочный коэффициент. Таким образом апостериорное распределение также будет бета-распределением, но с другими параметрами $P(p | n_s, N) = \mbox{Beta}(p; \alpha + n_s, \beta + N - n_s)$ . Априорные распределения с подобным свойством называют сопряженными априорными распределениями [ConjPrior].

$P(\mathcal{H} | \mathcal{D}) \propto P(\mathcal{D} | \mathcal{H}) P(\mathcal{H})$ $P(\mathcal{D} | \mathcal{H}) = P(n_s, N | p) = \mbox{Binom}(n_s, N | p) = C_{N}^{n_s} p^{n_s} (1-p)^{N-n_s}$ $P(\mathcal{H}) = P(p) = \mbox{Beta}(p; \alpha, \beta) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)} p^{\alpha-1}(1-p)^{\beta-1}$ $\begin{split}P(\mathcal{H} | \mathcal{D}) & = P(p | n_s, N) \\& \propto \mbox{Binom}(n_s, N | p) \mbox{Beta}(p; \alpha, \beta)\\& \propto C_{N}^{n_s} p^{n_s} (1-p)^{N-n_s}\frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)} p^{\alpha-1}(1-p)^{\beta-1}\\& \propto p^{n_s + \alpha - 1} (1-p)^{N - n_s + \beta - 1}\\& = \mbox{Beta}(p; \alpha + n_s, \beta + N - n_s)\end{split}$

Бета-распределения для различных параметров приведены на графике ниже [BetaDist, SciPyBeta]. При $\alpha = 1, \beta=1$ распределение однородное - эти значения удобно выбрать как априорные. Для остальных случаев максимум распределения в точке $p = (\alpha-1) / (\alpha + \beta - 2)$ . При увеличении $\alpha$ и $\beta$ распределение сужается и приближается к нормальному. Вместо $\alpha = 1, \beta=1$ начальные значения параметров $\alpha, \beta$ можно подобрать по историческим данным, чтобы априорное распределение конверсий совпадало с историческим значением.

Бета-распределение

x = np.linspace(0, 1, 1000)
fig = go.Figure()
a, b = 1, 1
fig.add_trace(go.Scatter(x=x, y=stats.beta.pdf(x, a=a, b=b), 
                             mode='lines', line_color='black', line_dash='dash',
                             name=f'a={a}, b={b}'))
a, b = 1, 5
fig.add_trace(go.Scatter(x=x, y=stats.beta.pdf(x, a=a, b=b), 
                             mode='lines', line_color='black', line_dash='solid',
                             name=f'a={a}, b={b}'))
a, b = 3, 5
fig.add_trace(go.Scatter(x=x, y=stats.beta.pdf(x, a=a, b=b), 
                             mode='lines', line_color='black', line_dash='solid',
                             name=f'a={a}, b={b}'))
a, b = 25, 30
fig.add_trace(go.Scatter(x=x, y=stats.beta.pdf(x, a=a, b=b), 
                             mode='lines', line_color='black', line_dash='solid',
                             name=f'a={a}, b={b}'))
a, b = 150, 50
fig.add_trace(go.Scatter(x=x, y=stats.beta.pdf(x, a=a, b=b), 
                             mode='lines', line_color='black', line_dash='solid',
                             name=f'a={a}, b={b}')) 
fig.add_trace(go.Scatter(
    x=[0.98, 0.08, 0.32, 0.55, 0.87],
    y=[1.35, 5.00, 2.80, 6.20, 12.0],
    mode="text",
    name=None,
    showlegend=False,
    text=["a=1, b=1", "a=1, b=5", "a=3, b=5", "a=25, b=30", "a=150, b=50"],
    textposition="top left"
))
fig.update_layout(title='Бета-распределение Beta(a, b)',
                  xaxis_title='x',
                  yaxis_title='Плотность вероятности',
                  showlegend=False,
                  xaxis_range=[0, 1],
                  height=550)
fig.show()

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

Оценка конверсии

def posterior_dist_binom(ns, ntotal, a_prior=1, b_prior=1):
    a = a_prior + ns
    b = b_prior + (ntotal - ns) 
    return stats.beta(a=a, b=b)
    
p = 0.1
nsample = 1000

exact_dist = stats.bernoulli(p=p)
data = exact_dist.rvs(nsample)
post_dist = posterior_dist_binom(ns=np.sum(data), ntotal=len(data))

x = np.linspace(0, 1, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=post_dist.pdf(x), line_color='black', name='Апостериорное'))
fig.add_trace(go.Scatter(x=[data.mean(), data.mean()], y=[0, max(post_dist.pdf(x))], 
                         line_color='black', mode='lines', line_dash='dash', name='Среднее в выборке'))
fig.add_trace(go.Scatter(x=[exact_dist.mean(), exact_dist.mean()], y=[0, max(post_dist.pdf(x))*1.05], 
                         line_color='red', mode='lines', line_dash='dash', name='Точное p'))
fig.update_layout(title='Апостериорное распределение',
                  xaxis_title='p',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[p-0.1, p+0.1],
                  hovermode="x",
                  height=500)
fig.show()

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

Для примера сравнения двух групп задается конверсия p_A , конверсия p_B выбирается на 5% больше. Для каждой группы генерируется nsample точек, по ним строятся апостериорные распределения конверсий. Сэмплированием из распределений оценивается вероятность конверсии группы Б больше группы А P(p_B > p_A) . На графике $P(p_B > p_A) = 84.0 \%$ . Т.к. количество точек nsample относительно небольшое, значения могут поменяться при повторном запуске.

Сравнение групп

def prob_pb_gt_pa(post_dist_A, post_dist_B, post_samp=100_000):
    sa = post_dist_A.rvs(size=post_samp)
    sb = post_dist_B.rvs(size=post_samp)
    b_gt_a = np.sum(sb > sa)
    return b_gt_a / post_samp

p_A = 0.1
p_B = p_A * 1.05
nsample = 1000

exact_dist_A = stats.bernoulli(p=p_A)
exact_dist_B = stats.bernoulli(p=p_B)
data_A = exact_dist_A.rvs(nsample)
data_B = exact_dist_B.rvs(nsample)

post_dist_A = posterior_dist_binom(ns=np.sum(data_A), ntotal=len(data_A))
post_dist_B = posterior_dist_binom(ns=np.sum(data_B), ntotal=len(data_B))

x = np.linspace(0, 1, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=post_dist_A.pdf(x), line_color='black', name='A'))
fig.add_trace(go.Scatter(x=x, y=post_dist_B.pdf(x), line_color='black', opacity=0.2, name='Б'))
fig.add_trace(go.Scatter(x=[exact_dist_A.mean(), exact_dist_A.mean()], y=[0, max(post_dist_A.pdf(x))*1.05], 
                         mode='lines', line_dash='dash', line_color='black', name='Точное A'))
fig.add_trace(go.Scatter(x=[exact_dist_B.mean(), exact_dist_B.mean()], y=[0, max(post_dist_A.pdf(x))*1.05], 
                         mode='lines', line_dash='dash', line_color='black', opacity=0.2, name='Точное Б'))
fig.update_layout(title='Апостериорные распределения',
                  xaxis_title='p',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[p_A/2, p_A*2],
                  hovermode="x",
                  height=500)
fig.show()

print(f'P(pB > pA): {prob_pb_gt_pa(post_dist_A, post_dist_B)}')

Апостериорные распределения конверсии в группах. Конверсия группы Б выше А с вероятностью 84%.

Динамика оценок конверсий и по мере набора данных иллюстрируется следующим примером. Сравнивается две группы, задается конверсия p_A , конверсия второй группы p_B выбирается на 5% больше. В каждой группе генерируется npoints точек nstep раз (всего N = npoints * nstep точек). По накопленным данным на каждом шаге строятся апостериорные распределения и считается вероятность P(p_B > p_A) . Также в обеих группах строятся и отображаются на графике центральные области апостериорных распределений, в которых лежит 95% распределения. По мере набора данных интервалы сужаются, вероятность стремится к 1 в пользу лучшей группы.

Динамика P(pb > pa)

def posterior_binom_approx_95pdi(post_dist):
    lower = post_dist.ppf(0.025)
    upper = post_dist.ppf(0.975)
    return lower, upper

pa = 0.1
pb = pa * 1.05

npoints = 1000
nstep = 150
sa = stats.binom.rvs(p=pa, n=npoints, size=nstep)
sb = stats.binom.rvs(p=pb, n=npoints, size=nstep)

df = pd.DataFrame()
df['npoints'] = [npoints] * nstep
df['sa_step'] = sa
df['sb_step'] = sb
df['N'] = df['npoints'].cumsum()
df['sa'] = df['sa_step'].cumsum()
df['sb'] = df['sb_step'].cumsum()
df['pa'] = df.apply(lambda r: posterior_dist_binom(r['sa'], r['N']).mean(), axis=1)
df[['pa_lower', 'pa_upper']] = df.apply(lambda r: posterior_binom_approx_95pdi(posterior_dist_binom(r['sa'], r['N'])), axis=1, result_type="expand")
df['pb'] = df.apply(lambda r: posterior_dist_binom(r['sb'], r['N']).mean(), axis=1)
df[['pb_lower', 'pb_upper']] = df.apply(lambda r: posterior_binom_approx_95pdi(posterior_dist_binom(r['sb'], r['N'])), axis=1, result_type="expand")
df['pb_gt_pa'] = df.apply(lambda r: prob_pb_gt_pa(posterior_dist_binom(r['sa'], r['N']), posterior_dist_binom(r['sb'], r['N']), post_samp=10_000), axis=1)


fig = go.Figure()
fig.add_trace(go.Scatter(x=df['N'], y=df['pa'], name='A',
                         line_color='black'))
fig.add_trace(go.Scatter(x=list(df['N']) + list(reversed(df['N'])), 
                         y=list(df['pa_upper']) + list(reversed(df['pa_lower'])),
                         fill="toself", name='A, 95% PDI', marker_color='black', opacity=0.2))
fig.add_trace(go.Scatter(x=df['N'], y=df['pb'], name='B',
                         line_color='blue'))
fig.add_trace(go.Scatter(x=list(df['N']) + list(reversed(df['N'])), 
                         y=list(df['pb_upper']) + list(reversed(df['pb_lower'])),
                         fill="toself", name='B, 95% PDI', marker_color='blue', opacity=0.2))
fig.update_layout(title='$p_A, p_B$',
                  yaxis_tickformat = ',.1%',
                  xaxis_title='N')
fig.show()

fig = go.Figure()
fig.add_trace(go.Scatter(x=df['N'], y=df['pb_gt_pa'], name='P(pb > pa)',
                         line_color='black'))
fig.update_layout(title='$P(p_B > p_A)$',
                  yaxis_range=[0, 1],
                  xaxis_title='N')
fig.show()

Оценки конверсий уточняются по мере набора данных.

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

Доля правильно угаданных вариантов демонстрируется следующим образом. В группе А задается конверсия p, в группе Б выбирается случайное значение в диапазоне $\pm 5\%$ от p. В группах генерируются данные с шагом n_samp_step. На каждом шаге считаются апостериорные распределения и P(p_B > p_A) . Эксперимент останавливается, если или P(p_A > p_B) достигает prob_stop=0.95 или сгенерировано максимальное количество точек n_samp_max. Процедура повторяется nexps раз, считается доля правильно угаданных групп во всех экспериментах. В данном случае в nexps = 100 правильно угадано 94. Точность 0.94 близка prob_stop = 0.95.

Nexp: 100, Correct Guesses: 94, Accuracy: 0.94

Доля правильно угаданных вариантов

cmp = pd.DataFrame(columns=['A', 'B', 'best_exact', 'exp_samp_size', 'A_exp', 'B_exp', 'best_exp', 'p_best'])

p = 0.1
nexps = 100
cmp['A'] = [p] * nexps
cmp['B'] = p * (1 + stats.uniform.rvs(loc=-0.05, scale=0.1, size=nexps))
cmp['best_exact'] = cmp.apply(lambda r: 'B' if r['B'] > r['A'] else 'A', axis=1)

n_samp_max = 10_000_000
n_samp_step = 10_000
prob_stop = 0.95

for i in range(nexps):
    pA = cmp.at[i, 'A']
    pB = cmp.at[i, 'B']
    exact_dist_A = stats.bernoulli(p=pA)
    exact_dist_B = stats.bernoulli(p=pB)
    n_samp_total = 0
    ns_A = 0
    ns_B = 0
    while n_samp_total < n_samp_max:
        dA = exact_dist_A.rvs(n_samp_step)
        dB = exact_dist_B.rvs(n_samp_step)
        n_samp_total += n_samp_step
        ns_A = ns_A + np.sum(dA)
        ns_B = ns_B + np.sum(dB)
        post_dist_A = posterior_dist_binom(ns=ns_A, ntotal=n_samp_total)
        post_dist_B = posterior_dist_binom(ns=ns_B, ntotal=n_samp_total)
        pb_gt_pa = prob_pb_gt_pa(post_dist_A, post_dist_B)
        best_gr = 'B' if pb_gt_pa >= prob_stop else 'A' if (1 - pb_gt_pa) >= prob_stop else None
        if best_gr:
            cmp.at[i, 'A_exp'] = post_dist_A.mean()
            cmp.at[i, 'B_exp'] = post_dist_B.mean()
            cmp.at[i, 'exp_samp_size'] = n_samp_total
            cmp.at[i, 'best_exp'] = best_gr
            cmp.at[i, 'p_best'] = pb_gt_pa
            break
    print(f'done {i}: nsamp {n_samp_total}, best_gr {best_gr}, P(b>a) {pb_gt_pa}')

cmp['correct'] = cmp['best_exact'] == cmp['best_exp']
display(cmp.head(10))
cor_guess = np.sum(cmp['correct'])
print(f"Nexp: {nexps}, Correct Guesses: {cor_guess}, Accuracy: {cor_guess / nexps}")

Средние

Байесовский подход требует предположения распределений сравниваемых величин. Выбор модели всегда произвольный и всегда остаются вопросы обоснования модели. Во многих случаях не нужно полное распределение, достаточно сравнивать средние: среднюю выручку на пользователя, среднюю длительностю просмотра и т.д. Для средних часто применима центральная предельная теорема [CLT], что позволяет использовать нормальное распределение [NormDist, SciPyNorm] в качестве функции правдоподобия даже при неизвестной форме исходного распределения.

Центральная предельная теорема формализует следующее наблюдение. Если взять произвольное распределение со средним значением $\mu$ и диспресий $\sigma^2$ , выбирать из него сэмплы длины и считать среднее в каждом сэмпле, то средние значения сэмплов будут распределены приблизительно нормально $Norm(\mu, \sigma^2/N)$ .

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

Есть несколько центральных предельных теорем [CLT]. Одна из формулировок следующая. Пусть есть последовательность независимых одинаково распределенных случайных величин $X_1, X_2, \dots$ с конечными математическим ожиданием $\mu$ и дисперсией $\sigma^2$ . Пусть $\overline{X}_N = (X_1 + \dots + X_N)/N$ их выборочное среднее. Тогда при стремящемся к бесконечности распределение центрированных и масштабированных выборочных средних сходится к нормальному распределению со средним значением 0 и дисперсией 1. Сходимость понимается как сходимость по распределению [RandVarsConv].

$\begin{split}P \left( \frac{\overline{X}_N - \mu}{\sigma / \sqrt{N}} = x \right) & \to Norm(x; 0, 1), \quad N \to \infty\\ Norm(x ; \mu, \sigma^2) & = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\tfrac{(x-\mu)^2}{2 \sigma^2} }\end{split}$

На графике приведено сравнение распределения выборочных средних с нормальным распределением на основе центральной предельной теоремы. Из гамма-распределения $P(x; \alpha, \beta) \propto x^{\alpha-1} \exp(-\beta x)$ [GammaDist, SciPyGamma] генерируется n_sample выборок по sample_len точек. В каждой выборке считаются средние, их распределение выводится на график. По точным значениям среднего и дисперсии исходного распределения определяются параметры нормального распределения центральной предельной теоремы clt_mu, clt_stdev. Это распределение также отображается на графике. Визуально распределение выборочных средних близко к нормальному распределению.

Выборочные средние гамма-распределения

a = 1
sample_len = 100
n_samples = 1000

exact_dist = stats.gamma(a=a)
samp = exact_dist.rvs(size=(n_samples, sample_len))
means = np.array([x.mean() for x in samp])
clt_mu = exact_dist.mean()
clt_stdev = exact_dist.std() / np.sqrt(sample_len)
means_stdev = means.std()

x = np.linspace(0, 10, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=exact_dist.pdf(x), 
                         mode='lines', line_color='black', line_dash='solid', name='Исходное распределение'))
fig.add_trace(go.Histogram(x=np.concatenate(samp), histnorm='probability density', name='Выборка', nbinsx=500,
                           marker_color='black', opacity=0.3))
fig.add_trace(go.Scatter(x=x, y=stats.norm.pdf(x, loc=clt_mu, scale=clt_stdev), 
                         mode='lines', line_color='black', line_dash='dash', name='$Norm(\mu, \sigma^2/N)$'))
fig.add_trace(go.Histogram(x=means, histnorm='probability density', name='Выборочные средние', nbinsx=50,
                           marker_color='green', opacity=0.5))
fig.update_layout(title='Выборочные средние',
                  xaxis_title='x',
                  yaxis_title='Плотность вероятности',
                  barmode='overlay',
                  hovermode="x",
                  height=550)
fig.update_layout(xaxis_range=[0, 5])
fig.show()

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

Центральная предельная теорема говорит о сходимости к нормальному распределению центрированных и масштабированных выборочных средних $\overline{X}_N$ при стремлении к бесконечности. Для конечного нормальное распределение не гарантируется. Отклонение распределения суммы конечного числа слагаемых от нормального дает теорема Берри-Эссеена [BerryEsseenTheorem]. Отличие зависит от количества слагаемых , дисперсии и коэффициента асимметрии исходного распределения. В приведенной формулировке центральная предельная теорема требует существования конечных среднего и дисперсии у исходного распределения. Примерами распределений, для которых эти свойства могут не выполнятся, являются распределение Парето [ParetoDist, SciPyPareto] и близкое к нему распределение Ломакса [LomaxDist, SciPyLomax]. Плотность вероятности последнего имеет вид

$P(x; c) = \frac{c}{(1 + x )^{c + 1}}, \quad x \ge 0, c > 0.$

При значениях параметра меньше или равном 2 дисперсия распределения Ломакса не является конечной. На графиках ниже приведена гистограмма n_samples выборочных средних с количеством слагаемых sample_len и нормальное распределение с параметрами из центральной предельной теоремы на основе точного распределения clt_mu, clt_stdev. Распределение выборочных средних скошено и сильнее отличается от нормального, чем в предыдущем случае. Перекос происходит из-за попадания в выборки больших значений из хвоста исходного распределения. На практике распределения ограничены, поэтому дисперсии и средние конечны. Центральная предельная теорема будет применима, но для приближения выборочных средних нормальным распределением потребуется большое количество точек.

Выборочные средние распределения Ломакса

c = 1.7
sample_len = 500
n_samples = 1000

exact_dist = stats.lomax(c=c)
samp = exact_dist.rvs(size=(n_samples, sample_len))
means = np.array([x.mean() for x in samp])
clt_mu = exact_dist.mean()
clt_stdev = exact_dist.std() / np.sqrt(sample_len)
means_stdev = means.std()

xaxis_max=10
x = np.linspace(0, xaxis_max, 2000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=exact_dist.pdf(x), 
                         mode='lines', line_color='black', line_dash='solid', name='Исходное распределение'))
fig.add_trace(go.Histogram(x=np.concatenate(samp)[np.concatenate(samp) < xaxis_max], histnorm='probability density', 
                           name='Выборка', nbinsx=500,
                           marker_color='black', opacity=0.3))
fig.add_trace(go.Scatter(x=x, y=stats.norm.pdf(x, loc=clt_mu, scale=means_stdev), 
                         mode='lines', line_color='black', line_dash='dash', name='$Norm(\mu, \sigma^2/N)$'))
fig.add_trace(go.Histogram(x=means, histnorm='probability density', name='Выборочные средние', nbinsx=150,
                          marker_color='green', opacity=0.5))
fig.update_layout(title='Выборочные средние',
                  xaxis_title='x',
                  yaxis_title='Плотность вероятности',
                  barmode='overlay',
                  hovermode="x",
                  xaxis_range=[0, xaxis_max],
                  height=550)
fig.show()

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

Для байесовской оценки параметров нормального распределения по выборке функция правдоподобия задается нормальным распределением $P(\mathcal{D} | \mathcal{H}) = Norm(x | \mu, \sigma_x^2)$ . У этой функции 2 параметра - $\mu$ и $\sigma_x$ . Для такой модели известно сопряженное априорное распределение [ConjPrior, Apx]. Ниже рассмотрен упрощенный вариант с одним параметром - подбирается только распределение $\mu$ , значение $\sigma_x$ фиксировано. Сопряженным априорным распределением $P(\mu)$ в этом случае также будет нормальное распределение $P(\mu) = Norm(\mu | \mu_0, \sigma_0^2)$ , но со своими параметрами $\mu_0$ , $\sigma_0$ . При расчете апостериорного распределения понадобится произведение функций правдоподобия по всем точкам x_i $P(\mathcal{H} | \mathcal{D}) \propto Norm(x_1 | \mu, \sigma_x^2) \cdots Norm(x_N | \mu, \sigma_x^2) Norm(\mu | \mu_0, \sigma_0^2)$ . В преобразованиях достаточно следить за зависимостью от $\mu$ , множители без $\mu$ войдут в нормировочный коэффициент итогового распределения. Апостериорное распределение сохранит нормальную форму $P(\mu | \mathcal{D}) = Norm(\mu | \mu_N, \sigma_N^2)$ с обновленными параметрами $\mu_N$ , $\sigma_N$ .

$\begin{split}P(\mathcal{D} | \mathcal{H}) & = Norm(x | \mu, \sigma_x^2) = \frac{1}{\sqrt{2 \pi \sigma_x^2}} e^{-\tfrac{(x - \mu)^2}{2 \sigma_x^2}}\\P(\mathcal{H}) & = Norm(\mu | \mu_0, \sigma_0^2) = \frac{1}{\sqrt{2 \pi \sigma_{0}^2}} e^{-\tfrac{(\mu-\mu_0)^2}{2 \sigma_{0}^2}} \\P(\mathcal{H} | \mathcal{D}) & \propto\prod_i^NNorm(x_i | \mu, \sigma_x^2)Norm(\mu | \mu_0, \sigma_0^2)\\& \propto_{\mu}\prod_i^Ne^{-\tfrac{(x_i - \mu)^2}{2 \sigma_x^2}}e^{-\tfrac{(\mu-\mu_0)^2}{2 \sigma_0^2}} \\& \propto_{\mu}e^{-\mu^2 \left[\tfrac{N}{2 \sigma_x^2} + \tfrac{1}{2 \sigma_0^2} \right] + 2\mu \left[\tfrac{\mu_0}{2 \sigma_0^2} + \tfrac{1}{2 \sigma_x^2} \sum_i^N x_i \right]}\\& \propto_{\mu}e^{-\tfrac{(\mu - \mu_N)^2}{2 \sigma_N^2}}= Norm(\mu | \mu_N, \sigma_N^2) \\ \sigma_N^2 &= \frac{\sigma_0^2 \sigma_x^2}{\sigma_x^2 + N \sigma_0^2},\quad\mu_N = \mu_0 \frac{\sigma_N^2}{\sigma_0^2} + \frac{\sigma_N^2}{\sigma_x^2} \sum_i^N x_i\end{split}$

Для проверки построения апостериорного распределения по данным задается нормальное распределение с параметрами mu, sigma и генерируется выборка размера nsample. Начальные параметры $\sigma_x, \sigma_0$ выбираются равными стандартному отклонению в выборке, $\mu_0$ - значению первой точки. Параметры $\mu_N, \sigma_N$ расчитываются по оставшимся точкам. Использовать всю выборку для задания начальных параметров некорректно, лучше использовать часть данных или исторические данные. На первом графике апостериорное распределение $\mu$ сравнивается с точным средним. Мода близка к среднему в выборке и точному среднему. На втором графике апостериорное прогнозное распределение сравнивается с исходным. Для построения распределения вначале нужно сгенерировать $\mu \sim Norm(\mu | \mu_N, \sigma_N^2)$ , после чего с этим $\mu$ сгенерировать $x \sim Norm(x | \mu, \sigma_x^2)$ . Гистограмма визуально близка исходному нормальному распределению. На последнем графике приведено сравнение распределений и $\mu$ . Это разные распределения - $P(x|\mathcal{D})$ приближает исходное нормальное распределение, $P(\mu|\mathcal{D})$ - оценка среднего. Распределение $\mu$ существенно уже. Его дисперсия задается $\sigma_N$ , тогда как дисперсия $P(x|\mathcal{D})$ определяется $\sigma_x$ и вариацией $\mu$ .

Оценка параметров нормального распределения

ConjugateNormalParams = namedtuple('ConjugateNormalParams', 'mu sigma sx')

def initial_params_normal(mu, sigma, sx):
    return ConjugateNormalParams(mu=mu, sigma=sigma, sx=sx)

def posterior_params_normal(data, initial_pars):
    N = len(data)
    sigma_n_2 = (initial_pars.sigma**2 * initial_pars.sx**2) / (initial_pars.sx**2 + N * initial_pars.sigma**2)
    mu_n = initial_pars.mu * sigma_n_2 / initial_pars.sigma**2 + np.sum(data) * sigma_n_2 / initial_pars.sx**2    
    return ConjugateNormalParams(mu=mu_n, sigma=np.sqrt(sigma_n_2), sx=initial_pars.sx)

def posterior_mu_dist(params):
    return stats.norm(loc=params.mu, scale=params.sigma)

def posterior_rvs(params, nsamp):
    mus = stats.norm.rvs(loc=params.mu, scale=params.sigma, size=nsamp)
    return stats.norm.rvs(loc=mus, scale=params.sx, size=nsamp)

mu = 3
sigma = 1
nsample = 1000
npostsamp = 100000

exact_dist = stats.norm(loc=mu, scale=sigma)
data = exact_dist.rvs(nsample)

sx = np.std(data)
mu0 = data[0]
sigma0 = np.std(data)
pars = initial_params_normal(mu=mu0, sigma=sigma0, sx=sx)
pars = posterior_params_normal(data[1:], pars)
post_mu = posterior_mu_dist(pars)
post_samp = posterior_rvs(pars, npostsamp)

x = np.linspace(0, 10, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=post_mu.pdf(x), line_color='black', name='$\mbox{Апостериорное }\mu$'))
fig.add_trace(go.Scatter(x=[data.mean(), data.mean()], y=[0, max(post_mu.pdf(x))], 
                         line_color='black', mode='lines', line_dash='dash', name='Среднее в выборке'))
fig.add_trace(go.Scatter(x=[exact_dist.mean(), exact_dist.mean()], y=[0, max(post_mu.pdf(x))*1.05], 
                         line_color='red', mode='lines', line_dash='dash', name='Точное среднее'))
fig.update_layout(title='$\mbox{Апостериорное распределение } \mu$',
                  xaxis_title='$\mu$',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[2, 4],
                  barmode='overlay',
                  hovermode="x",
                  height=500)                  
fig.show()

fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=exact_dist.pdf(x), line_color='red', line_dash='dash', name='Точное распределение'))
fig.add_trace(go.Histogram(x=post_samp, histnorm='probability density', name='Сэмпл апостериорного x', nbinsx=100,
                           marker_color='black', opacity=0.8))
fig.update_layout(title='Сэмпл апостериорного распределения x',
                  xaxis_title='x',
                  yaxis_title='Плотность вероятности',
                  #xaxis_range=[0, 10],
                  barmode='overlay',
                  hovermode="x",
                  height=500)                  
fig.show()

fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=post_mu.pdf(x), line_color='black', name='$\mbox{Апостериорное }\mu$'))
fig.add_trace(go.Scatter(x=x, y=exact_dist.pdf(x), line_color='red', line_dash='dash', name='Точное распределение'))
fig.add_trace(go.Histogram(x=post_samp, histnorm='probability density', name='Сэмпл апостериорного x', nbinsx=100,
                           marker_color='black', opacity=0.8))
fig.update_layout(title='$\mbox{Распределения } x \mbox{ и } \mu$',
                  xaxis_title='x',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[0, 6],
                  barmode='overlay',
                  hovermode="x",
                  height=500)                  
fig.show()

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

Апостериорное предиктивное распределение x близко исходному нормальному распределению.

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

Оценка среднего произвольного распределения показана для гамма-распределения. Выбирается nsample точек, сэмпл разбивается на части по nsplit штук, в каждой части считается среднее. Выборочные средние means предполагаются распределенными нормально, к ним применяется байесовское моделирование. Начальные параметры $\sigma_x$ и $\sigma_0$ задаются равными стандартному отклонению выборочных средних sx = np.std(means), sigma0 = sx, $\mu_0$ - первому выборочному среднему mu0 = means[0]. Количество точек nsplit = 100 выбрано произвольно, точнее можно оценить по теореме Берри-Эссеена. Вместо деления данных на части по nsplit можно посчитать среднее по всей выборке. В таком случае распределение $P(\mu|\mathcal{D})$ нужно будет оценивать по одной точке, что не позволит валидировать модель. Количество точек nsplit лучше считать гиперпараметром модели. На первом графике показано исходное распределение и распределение выборочных средних, которое ожидается приблизительно нормальным. На втором построено апостериорное распределение $\mu$ , мода близка среднему в выборке и точному среднему. На третьем графике приведено апостериорное прогнозное распределение выборочных средних. Оно визуально близко нормальному распределению с параметрами на основе центральной предельной теоремы. Как и в первом примере, распределение $\mu$ (второй график) уже распределение выборочных средних (третий график). Для сравнения средних нужно ориентироваться на распределения $\mu$ .

Оценка выборочных средних

def reshape_and_compute_means(sample, n_split):
    n_means = len(sample) // n_split
    samp_reshaped = np.reshape(sample[0 : n_means * n_split], (n_means, n_split))
    means = np.array([x.mean() for x in samp_reshaped])
    return means

def exact_clt_dist(exact_dist, n_split):
    clt_mu = exact_dist.mean()
    clt_stdev = exact_dist.std() / np.sqrt(n_split)
    return stats.norm(loc=clt_mu, scale=clt_stdev)

def sample_clt_dist(means):
    clt_mu = means.mean()
    clt_std = means.std()
    return stats.norm(loc=clt_mu, scale=clt_std)

nsample = 50000
nsplit = 100
npostsamp = 50000

a = 1
b = 2
exact_dist = stats.gamma(a=a, scale=1/b)
data = exact_dist.rvs(nsample)

means = reshape_and_compute_means(data, nsplit)
clt_dist_exact = exact_clt_dist(exact_dist, nsplit)

sx = np.std(means)
mu0 = means[0]
sigma0 = sx
pars = initial_params_normal(mu=mu0, sigma=sigma0, sx=sx)
pars = posterior_params_normal(means[1:], pars)
post_mu = posterior_mu_dist(pars)
post_samp = posterior_rvs(pars, npostsamp)

x = np.linspace(0, 10, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=exact_dist.pdf(x), 
                         mode='lines', line_color='black', line_dash='solid', name='Исходное распределение'))
fig.add_trace(go.Scatter(x=x, y=clt_dist_exact.pdf(x), 
                         mode='lines', line_color='black', line_dash='dash', name='$Norm(\mu, \sigma^2/N)$'))
fig.add_trace(go.Histogram(x=means, histnorm='probability density', name='Выборочные средние', nbinsx=50,
                           marker_color='green', opacity=0.5))
fig.update_layout(title='Выборочные средние',
                  xaxis_title='x',
                  yaxis_title='Плотность вероятности',
                  barmode='overlay',
                  hovermode="x",
                  height=550)
fig.update_layout(xaxis_range=[0, 3])
fig.show()

x = np.linspace(0, 4, 10000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=post_mu.pdf(x), line_color='black', name='$\mbox{Распределение }\mu$'))
fig.add_trace(go.Scatter(x=[data.mean(), data.mean()], y=[0, max(post_mu.pdf(x))], 
                         line_color='black', mode='lines', line_dash='dash', name='Среднее в выборке'))
fig.add_trace(go.Scatter(x=[exact_dist.mean(), exact_dist.mean()], y=[0, max(post_mu.pdf(x))*1.05], 
                         line_color='red', mode='lines', line_dash='dash', name='Точное среднее'))
fig.update_layout(title='$\mbox{Распределение }\mu$',
                  xaxis_title='$\mu$',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[exact_dist.mean()-0.1, exact_dist.mean()+0.1],
                  barmode='overlay',
                  hovermode="x",
                  height=500)                  
fig.show()

fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=exact_dist.pdf(x), line_dash='solid', line_color='black', name='Исходное распределение'))
fig.add_trace(go.Scatter(x=x, y=clt_dist_exact.pdf(x), 
                         mode='lines', line_color='black', line_dash='dash', name='$Norm(\mu, \sigma^2/N)$'))
fig.add_trace(go.Histogram(x=post_samp, histnorm='probability density', name='$\mbox{Апострериорное } \overline{X}_N$', nbinsx=300,
                           marker_color='black', opacity=0.2))
fig.update_layout(title='$\mbox{Апостериорное распределение } \overline{X}_N$',
                  xaxis_title='x',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[0, 3],
                  barmode='overlay',
                  hovermode="x",
                  height=500)                  
fig.show()

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

Оценка mu по выборочным средним. Мода распределения близка точному среднему гамма-распределения.

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

Для примера сравнения групп задаются 2 гамма-распределения. Параметры одинаковы, параметр группы Б на $5\%$ меньше А. Из каждого делается выборка размера nsample, считаются выборочные средние. По выборочным средним строятся апостериорные распределения $\mu$ . По этим распределениям вычисляется вероятность среднего в группе Б больше группы А $P(\mu_B > \mu_A)$ . На первом графике приведены исходные распределения и точные средние. На втором - распределения $\mu$ и точные средние. При выбранных сэмплах распределения $P(\mu|\mathcal{D})$ пересекаются слабо $P(\mu_B > \mu_A) = 1$ .

Сравнение средних

def prob_pb_gt_pa(post_dist_A, post_dist_B, post_samp=100_000):
    sa = post_dist_A.rvs(size=post_samp)
    sb = post_dist_B.rvs(size=post_samp)
    b_gt_a = np.sum(sb > sa)
    return b_gt_a / post_samp

nsample = 30000
npostsamp = 50000
nsplit = 100

A, B = {}, {}
A['dist_pars'] = {'a': 1, 'b': 2}
B['dist_pars'] = {'a': 1, 'b': 2*0.95}
for g in [A, B]:
    g['exact_dist'] = stats.gamma(a=g['dist_pars']['a'], scale=1/g['dist_pars']['b'])
    g['data'] = g['exact_dist'].rvs(nsample)
    g['means'] = reshape_and_compute_means(g['data'], nsplit)
    g['post_pars'] = initial_params_normal(mu=g['means'][0], sigma=np.std(g['means']), sx=np.std(g['means']))
    g['post_pars'] = posterior_params_normal(g['means'][1:], g['post_pars'])
    g['post_mu'] = posterior_mu_dist(g['post_pars'])
    g['post_samp'] = posterior_rvs(g['post_pars'], npostsamp)

    
x = np.linspace(0, 3, 5000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=A['exact_dist'].pdf(x), line_color='black', opacity=0.2, name='A'))
fig.add_trace(go.Scatter(x=x, y=B['exact_dist'].pdf(x), line_color='black', name='Б'))
fig.add_trace(go.Scatter(x=[A['exact_dist'].mean(), A['exact_dist'].mean()], y=[0, max(A['exact_dist'].pdf(x))*1.05], 
                         mode='lines', line_dash='dash', line_color='black', opacity=0.2, name='Точное среднее A'))
fig.add_trace(go.Scatter(x=[B['exact_dist'].mean(), B['exact_dist'].mean()], y=[0, max(B['exact_dist'].pdf(x))*1.05], 
                         mode='lines', line_dash='dash', line_color='black', name='Точное среднее Б'))
fig.update_layout(title='Исходные распределения',
                  xaxis_title='x',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[0, 3],
                  hovermode="x",
                  height=500)
fig.show()

fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=A['post_mu'].pdf(x), line_color='black', opacity=0.2, name='A'))
fig.add_trace(go.Scatter(x=x, y=B['post_mu'].pdf(x), line_color='black', name='Б'))
fig.add_trace(go.Scatter(x=[A['exact_dist'].mean(), A['exact_dist'].mean()], y=[0, max(A['post_mu'].pdf(x))*1.05], 
                         mode='lines', line_dash='dash', line_color='black', opacity=0.2, name='Точное среднее A'))
fig.add_trace(go.Scatter(x=[B['exact_dist'].mean(), B['exact_dist'].mean()], y=[0, max(B['post_mu'].pdf(x))*1.05], 
                         mode='lines', line_dash='dash', line_color='black', name='Точное среднее Б'))
fig.update_layout(title='$\mbox{Распределения } \mu$',
                  xaxis_title='$\mu$',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[A['exact_dist'].mean()-0.1, A['exact_dist'].mean()+0.1],
                  hovermode="x",
                  height=500)
fig.show()

print(f"P(mu_B > mu_A): {prob_pb_gt_pa(A['post_mu'], B['post_mu'])}")

Исходные гамма-распределения и точные средние. Среднее Б больше А.

Оценки mu по выборочным средним. По собранным данным среднее Б больше А с вероятностью 1.

Для демонстрации доли правильно угаданных вариантов задается 2 группы с гамма-распределениями. В группе A параметры гамма-распределения фиксированы, в Б параметр меняется в пределах $\pm5\%$ группы A. Вместе с параметрами меняются средние. Из распределений генерируются данные с шагом n_samp_step. На каждом шаге считаются выборочные средние nsplit. По выборочным средним проводится оценка параметров $\mu$ . Распределения этих параметров сравниваются. Набор данных останавливается, если $P(\mu_B > \mu_A)$ или $P(\mu_A > \mu_B)$ достигает prob_stop = 0.95, или достигнуто предельное количество точек n_samp_max. Проводится nexps экспериментов, считается доля правильно угаданных групп с большим средним. В примере доля 0.97 близка prob_stop = 0.95.

Nexp: 100, Correct Guesses: 97, Accuracy: 0.97

Доля правильно угаданных вариантов

nexps = 100
prob_stop = 0.95
nsplit = 100
n_samp_max = 1_000_000
n_samp_step = 10_000

A = {'a': 1, 'b': 2}

cmp = pd.DataFrame(columns=['A_pars', 'B_pars', 'A_mean', 'B_mean', 'best_exact', 'exp_samp_size', 'A_exp', 'B_exp', 'best_exp', 'p_best'])
cmp['A_pars'] = [A] * nexps
cmp['B_pars'] = cmp['A_pars'].apply(lambda x: {'a': x['a'], 'b': x['b'] * (1 + stats.uniform.rvs(loc=-0.05, scale=0.1))})
cmp['A_mean'] = cmp['A_pars'].apply(lambda x: stats.gamma(a=x['a'], scale=1/x['b']).mean())
cmp['B_mean'] = cmp['B_pars'].apply(lambda x: stats.gamma(a=x['a'], scale=1/x['b']).mean())
cmp['best_exact'] = cmp.apply(lambda r: 'B' if r['B_mean'] > r['A_mean'] else 'A', axis=1)

for i in range(nexps):
    A_pars = cmp.at[i, 'A_pars']
    B_pars = cmp.at[i, 'B_pars']
    exact_dist_A = stats.gamma(a=A_pars['a'], scale=1/A_pars['b'])
    exact_dist_B = stats.gamma(a=B_pars['a'], scale=1/B_pars['b'])
    n_samp_total = 0
    dA = []
    dB = []
    while n_samp_total < n_samp_max:
        dA.extend(exact_dist_A.rvs(n_samp_step))
        dB.extend(exact_dist_B.rvs(n_samp_step))
        n_samp_total += n_samp_step
        means_A = reshape_and_compute_means(dA, nsplit)
        post_pars_A = initial_params_normal(mu=means_A[0], sigma=np.std(means_A), sx=np.std(means_A))
        post_pars_A = posterior_params_normal(means_A[1:], post_pars_A)
        post_mu_A = posterior_mu_dist(post_pars_A)
        means_B = reshape_and_compute_means(dB, nsplit)
        post_pars_B = initial_params_normal(mu=means_B[0], sigma=np.std(means_B), sx=np.std(means_B))
        post_pars_B = posterior_params_normal(means_B[1:], post_pars_B)
        post_mu_B = posterior_mu_dist(post_pars_B)
        pb_gt_pa = prob_pb_gt_pa(post_mu_A, post_mu_B)
        best_gr = 'B' if pb_gt_pa >= prob_stop else 'A' if (1 - pb_gt_pa) >= prob_stop else None
        if best_gr:
            cmp.at[i, 'A_exp'] = post_mu_A.mean()
            cmp.at[i, 'B_exp'] = post_mu_B.mean()
            cmp.at[i, 'exp_samp_size'] = n_samp_total
            cmp.at[i, 'best_exp'] = best_gr
            cmp.at[i, 'p_best'] = pb_gt_pa
            break
    print(f'done {i}: nsamp {n_samp_total}, best_gr {best_gr}, P(B>A) {pb_gt_pa}')


cmp['correct'] = cmp['best_exact'] == cmp['best_exp']
display(cmp.head(10))
cor_guess = np.sum(cmp['correct'])
print(f"Nexp: {nexps}, Correct Guesses: {cor_guess}, Accuracy: {cor_guess / nexps}")

Выручка на пользователя

Для оценки денежного эффекта сравнивают выручку на пользователя в группах $P_{пользователи}(x)$ . Удобно выделить выручку на платящего $P_{платящие}(x)$ . При конверсии в оплату распределение ненулевой выручки на пользователя $p P_{платящие}(x)$ , с вероятностью 1-p выручка нулевая.

$P_{пользователи}(x) = \begin{cases}1-p, \, x = 0\\p P_{платящие}(x), \, x > 0\end{cases}$

Оценка конверсии в оплату делалась ранее. Выручку на платящего можно моделировать логнормальным распределением [LognormDist, SciPyLognorm] или распределением Парето [ParetoDist, SciPyPareto] по аналогии с распределением богатства. Для транзакционных сервисов, в частности маркетплейсов, лучше подходит логнормальное распределение. Случайная величина логнормальная $X \sim Lognormal(\mu, s^2)$ , если логарифм распределен нормально $\ln(X) \sim Norm(\mu, s^2)$ . Плотность вероятности приведена ниже

$\begin{split}P(x) & = \frac{1}{x s \sqrt{2 \pi}} e^{-\tfrac{(\ln(x) - \mu)^2}{2 s^2}}\end{split}$

Логнормальное распределение

x = np.linspace(0, 20, 2000)
fig = go.Figure()
mu, s = 1, 1
fig.add_trace(go.Scatter(x=x, y=stats.lognorm.pdf(x, s=s, scale=np.exp(mu)), 
                             mode='lines', line_color='black', line_dash='solid',
                             name=f'$\mu={mu}, \, s={s}$'))
mu, s = 2, 1
fig.add_trace(go.Scatter(x=x, y=stats.lognorm.pdf(x, s=s, scale=np.exp(mu)), 
                             mode='lines', line_color='black', line_dash='solid',
                             name=f'$\mu={mu}, \, s={s}$'))
mu, s = 1, 2
fig.add_trace(go.Scatter(x=x, y=stats.lognorm.pdf(x, s=s, scale=np.exp(mu)), 
                             mode='lines', line_color='black', line_dash='solid',
                             name=f'$\mu={mu}, \, s={s}$'))
fig.add_trace(go.Scatter(
    x=[2.90, 11.2, 1.80],
    y=[0.25, 0.06, 0.48],
    mode="text",
    name=None,
    showlegend=False,
    text=["$\mu=1, \, s=1$", "$\mu=2, \, s=1$", "$\mu=1, \, s=2$"],
    textposition="middle center"
))
fig.update_layout(title='Логнормальное распределение',
                  xaxis_title='x',
                  yaxis_title='Плотность вероятности',
                  hovermode="x",
                  showlegend=False,
                  height=550)
fig.show()

Сопряженное априорное распределение к логнормальной функции правдоподобия $P(\mathcal{D} | \mathcal{H}) = Lognorm(x | \mu, s^2)$ строится аналогично нормальному распределению [ConjPrior]. Для упрощенной модели с одним параметром $\mu$ и фиксированным сопряженное априорное распределение нормальное $P(\mu) = Norm(\mu | \mu_0, \sigma_0^2)$ с параметрами $\mu_0$ , $\sigma_0$ . Апостериорное распределение нормальное $P(\mu | \mathcal{D}) = Norm(\mu | \mu_N, \sigma_N^2)$ с обновленными параметрами $\mu_N$ , $\sigma_N$ . В $\mu_N$ суммируются логарифмы точек выборки.

$\begin{split}P(\mathcal{D} | \mathcal{H}) & = Lognorm(x | \mu, s^2) = \frac{1}{x \sqrt{2 \pi s^2}} e^{-\tfrac{(\ln x - \mu)^2}{2 s^2}}\\P(\mathcal{H}) & = Norm(\mu | \mu_0, \sigma_0^2) = \frac{1}{\sqrt{2 \pi \sigma_{0}^2}} e^{-\tfrac{(\mu-\mu_0)^2}{2 \sigma_{0}^2}} \\P(\mathcal{H} | \mathcal{D}) & \propto\prod_i^NLognorm(x_i | \mu, s^2)Norm(\mu | \mu_0, \sigma_0^2)\\& \propto_{\mu}\prod_i^Ne^{-\tfrac{(\ln x_i - \mu)^2}{2 s^2}}e^{-\tfrac{(\mu-\mu_0)^2}{2 \sigma_0^2}} \\& \propto_{\mu}e^{-\mu^2 \left[\tfrac{N}{2 s^2} + \tfrac{1}{2 \sigma_0^2} \right] + 2\mu \left[\tfrac{\mu_0}{2 \sigma_0^2} + \tfrac{1}{2 s^2} \sum_i^N \ln x_i \right]}\\& \propto_{\mu}e^{-\tfrac{(\mu - \mu_N)^2}{2 \sigma_N^2}}= Norm(\mu | \mu_N, \sigma_N^2) \\ \sigma_N^2 &= \frac{\sigma_0^2 s^2}{s^2 + N \sigma_0^2},\quad\mu_N = \mu_0 \frac{\sigma_N^2}{\sigma_0^2} + \frac{\sigma_N^2}{s^2} \sum_i^N \ln x_i\end{split}$

Для примера построения апостериорного распределения по данным из логнормального распределения с параметрами mu, s генерируется выборка размера nsample. Считается логарифм точек выборки. Параметры и $\sigma_0$ выбираются равными стандартному отклонению логарифма точек выборки, $\mu_0$ задается равным значению первой точки. По оставшимся точкам расчитываются $\mu_N, \sigma_N$ . Апостериорное распределение $\mu$ приведено на первом графике. Среднее логнормального распределения задается выражением $E[x] = \exp(\mu + s^2/2)$ , поэтому $\mu + s^2/2$ должно быть оценкой логарифма точного среднего. Т.к. $\mu$ распределено нормально $P(\mu) = Norm(\mu_N, \sigma_N^2)$ , величина $\mu + s^2/2$ также распределена нормально $Norm(\mu_N + s^2/2, \sigma_N^2)$ . На первом графике мода распределения $\mu + s^2/2$ близка логарифму среднего в выборке и точного среднего. На втором графике апостериорное прогнозное распределение сравнивается с исходным. Гистограмма близка исходному.

$\begin{split}P(x) & = Lognorm(x | \mu, s^2)\\E[x] & = e^{\mu + s^2/2}\\P(\mu) & = Norm(\mu | \mu_N, \sigma_N^2)\\P_{\mu + s^2/2}(y) & = Norm(y | \mu_N + s^2/2, \sigma_N^2)\end{split}$

Оценка параметров логнормального распределения

ConjugateLognormalParams = namedtuple('ConjugateLognormalParams', 'mu sigma sx')

def initial_params_lognormal(mu, sigma, sx):
    return ConjugateLognormalParams(mu=mu, sigma=sigma, sx=sx)

def posterior_params_lognormal(data, initial_pars):
    N = len(data)
    lnx = np.log(data)
    sigma_n_2 = (initial_pars.sigma**2 * initial_pars.sx**2) / (initial_pars.sx**2 + N * initial_pars.sigma**2)
    mu_n = initial_pars.mu * sigma_n_2 / initial_pars.sigma**2 + np.sum(lnx) * sigma_n_2 / initial_pars.sx**2    
    return ConjugateLognormalParams(mu=mu_n, sigma=np.sqrt(sigma_n_2), sx=initial_pars.sx)

def posterior_mu_dist_lognormal(params):
    return stats.norm(loc=params.mu, scale=params.sigma)

def posterior_lognormal_rvs(params, nsamp):
    mus = stats.norm.rvs(loc=params.mu, scale=params.sigma, size=nsamp)
    return stats.lognorm.rvs(s=params.sx, scale=np.exp(mus), size=nsamp)

def posterior_mean_dist_lognormal(params):
    return stats.lognorm(scale=np.exp(params.mu + params.sx**2/2), s=params.sigma)

def posterior_ln_mean_dist_lognormal(params):
    return stats.norm(loc=params.mu + params.sx**2/2, scale=params.sigma)
    
s = 1
mu = 1.5
nsample = 1000

exact_dist = stats.lognorm(s=s, scale=np.exp(mu))
data = exact_dist.rvs(nsample)

lnx = np.log(data)
sx = np.std(lnx)
mu0 = lnx[0]
sigma0 = sx

pars = initial_params_lognormal(mu=mu0, sigma=sigma0, sx=sx)
pars = posterior_params_lognormal(data[1:], pars)
post_mu = posterior_mu_dist_lognormal(pars)
post_lnmeans = posterior_ln_mean_dist_lognormal(pars)
npostsamp = 10000
post_samp = posterior_lognormal_rvs(pars, npostsamp)

x = np.linspace(0, 4, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=post_mu.pdf(x), line_color='black', name='$\mbox{Распределение }\mu$'))
fig.add_trace(go.Scatter(x=x, y=post_lnmeans.pdf(x), line_color='black', opacity=0.2, name='$\mbox{Распределение } \mu + s^2/2$'))
fig.add_trace(go.Scatter(x=[np.log(data.mean()), np.log(data.mean())], y=[0, max(post_mu.pdf(x))], 
                         line_color='black', mode='lines', line_dash='dash', name='Логарифм среднего в выборке'))
fig.add_trace(go.Scatter(x=[np.log(exact_dist.mean()), np.log(exact_dist.mean())], y=[0, max(post_mu.pdf(x))*1.05], 
                         line_color='red', mode='lines', line_dash='dash', name='Логарифм точного среднего'))
fig.update_layout(title='$\mbox{Апостериорные распределения } \mu \mbox{ и } \mu + s^2/2$',
                  xaxis_title='$\mu$',
                  yaxis_title='Плотность вероятности',
                  #xaxis_range=[2, 4],
                  barmode='overlay',
                  hovermode="x",
                  height=500)                  
fig.show()

xaxis_max=20
x = np.linspace(0, xaxis_max, 10000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=exact_dist.pdf(x), line_dash='dash', line_color='red', name='Точное распределение'))
fig.add_trace(go.Histogram(x=post_samp[post_samp < xaxis_max], histnorm='probability density', name='Сэмпл апостериорного', nbinsx=100,
                          marker_color='black', opacity=0.8))
fig.update_layout(title='$\mbox{Апостериорное распределение } x$',
                  xaxis_title='$x$',
                  yaxis_title='Плотность вероятности',
                  #xaxis_range=[0, 10],
                  barmode='overlay',
                  hovermode="x",
                  height=500)                  
fig.show()

Апостериорное распределение mu + s^2/2 близко логарифму точного среднего.

Апостериорное распределение x близко исходному распределению.

Для большей ожидаемой выручки на платящего нужно сравнивать группы по $E[x]=\exp(\mu + s^2/2)$ . Достаточно сравнивать $\mu + s^2/2$ . Распределение этой величины нормальное $Norm(\mu + s^2/2 | \mu_N, \sigma_N)$ . В примере задаются 2 логнормальных распределения: одно с параметрами s, mu, во втором mu на 5% больше. Генерируется nsample точек. Строятся апостериорные распределения. Вероятность P(E[x]_B > E[x]_A) ожидаемой выручки на пользователя группы Б больше А в данном случае близка 1. На первом графике приведены исходные распределения и их точные средние. На втором - распределения $Norm(\mu + s^2/2 | \mu_N, \sigma_N)$ в группах.

Сравнение групп

def prob_pb_gt_pa(post_dist_A, post_dist_B, post_samp=100_000):
    sa = post_dist_A.rvs(size=post_samp)
    sb = post_dist_B.rvs(size=post_samp)
    b_gt_a = np.sum(sb > sa)
    return b_gt_a / post_samp

nsample = 3000
npostsamp = 50000

A, B = {}, {}
s = 1
mu = 2
A['dist_pars'] = {'s': s, 'scale': np.exp(mu)}
B['dist_pars'] = {'s': s, 'scale': np.exp(mu * 1.05)}
for g in [A, B]:
    g['exact_dist'] = stats.lognorm(s=g['dist_pars']['s'], scale=g['dist_pars']['scale'])
    g['data'] = g['exact_dist'].rvs(nsample)
    g['post_pars'] = initial_params_lognormal(mu=np.log(g['data'])[0], sigma=np.std(np.log(g['data'])), sx=np.std(np.log(g['data'])))
    g['post_pars'] = posterior_params_lognormal(g['data'][1:], g['post_pars'])
    g['post_ln_means_dist'] = posterior_ln_mean_dist_lognormal(g['post_pars'])
    
x = np.linspace(0, 30, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=A['exact_dist'].pdf(x), line_color='black', opacity=0.2, name='Исходное A'))
fig.add_trace(go.Scatter(x=x, y=B['exact_dist'].pdf(x), line_color='black', name='Исходное Б'))
fig.add_trace(go.Scatter(x=[A['exact_dist'].mean(), A['exact_dist'].mean()], y=[0, max(A['exact_dist'].pdf(x))*1.05], 
                         mode='lines', line_dash='dash', line_color='red', opacity=0.2, name='Точное среднее A'))
fig.add_trace(go.Scatter(x=[B['exact_dist'].mean(), B['exact_dist'].mean()], y=[0, max(B['exact_dist'].pdf(x))*1.05], 
                         mode='lines', line_dash='dash', line_color='red', name='Точное среднее Б'))
fig.update_layout(title='Исходные распределения',
                  xaxis_title='x',
                  yaxis_title='Плотность вероятности',
                  hovermode="x",
                  height=500)
fig.show()

x = np.linspace(0, 3, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=A['post_ln_means_dist'].pdf(x), line_color='black', opacity=0.2, name='A'))
fig.add_trace(go.Scatter(x=x, y=B['post_ln_means_dist'].pdf(x), line_color='black', name='Б'))
fig.add_trace(go.Scatter(x=[np.log(A['exact_dist'].mean()), np.log(A['exact_dist'].mean())], y=[0, max(A['post_ln_means_dist'].pdf(x))*1.05], 
                         mode='lines', line_dash='dash', line_color='red', opacity=0.3, name='Логарифм точного среднего А'))
fig.add_trace(go.Scatter(x=[np.log(B['exact_dist'].mean()), np.log(B['exact_dist'].mean())], y=[0, max(B['post_ln_means_dist'].pdf(x))*1.05], 
                         mode='lines', line_dash='dash', line_color='red', name='Логарифм точного среднего Б'))
fig.update_layout(title='$\mbox{Распределения } \mu + s^2/2$',
                  xaxis_title='$\mu$',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[2, 3],
                  hovermode="x",
                  height=500)
fig.show()

print(f"P(E[x]_B > E[x]_A): {prob_pb_gt_pa(A['post_ln_means_dist'], B['post_ln_means_dist'])}")

Исходные распределения и точные средние.

Распределения оценок логарифма среднего mu + s^2/2. Среднее Б больше А с вероятностью 1.

Доля правильно угаданных вариантов проверяется для выручки на пользователя $P_{пользователи}(x)$ . В группе А задается конверсия и выручка на платящего $\mu, s$ . В группе Б параметры и $\mu$ случайно выбираются в пределах $\pm 5\%$ от A. Параметры меняются независимо, на практике они чаще меняются согласованно в противоположные стороны. Группы сравниваются по большей средней выручке на пользователя $E_{пользователи}[x] = p \exp(\mu + s^2/2)$ . Для оценки используется бета-распределение $P(p) = \mbox{Beta}(p; \alpha + n_s, \beta + N - n_s)$ , где - общее количество пользователей, n_s - количество платящих. Выручка на платящего моделируется логнормальным распределением. Т.к. распределение $\mu + s^2/2$ нормальное, распределение $\exp(\mu + s^2/2)$ логнормальное $P_{\exp(\mu + s^2/2)}(y) = Lognorm(y | \mu_N + s^2/2, \sigma_N^2)$ . Распределение $p\exp(\mu + s^2/2)$ будет произведением бета- и логнормального распределений $P_{p\exp(\mu + s^2/2)} \sim \mbox{Beta}(p; \alpha + n_s, \beta + N - n_s) Lognorm(y ; \mu_N + s^2/2, \sigma_N^2)$ . Данные в каждом эксперименте добавляются с шагом n_samp_step. Эксперимент останавливается при достижении вероятности среднего в одной группе больше другой prob_stop или при наборе n_samp_max точек. При малом n_samp_step доля правильно угаданных вариантов может быть ниже prob_stop, что можно объяснить неточностью модели и попаданием выбросов. При достаточно большом n_samp_step доля правильно угаданных вариантов соответствует ожидаемой prob_stop.

$\begin{split}P_{пользователи}(x) & = \begin{cases}1-p, \, x = 0\\p P_{платящие}(x), \, x > 0\end{cases}= \begin{cases}1-p, \, x = 0\\p Lognorm(x | s, \mu_N, \sigma_N), \, x > 0\end{cases}\\E_{пользователи}[x] & = p e^{\mu + s^2/2}\\P(p) & = \mbox{Beta}(p; \alpha + n_s, \beta + N - n_s),\\P_{\exp(\mu + s^2/2)}(y) & = Lognorm(y | \mu_N + s^2/2, \sigma_N^2)\\P_{p\exp(\mu + s^2/2)} & \sim \mbox{Beta}(p; \alpha + n_s, \beta + N - n_s) Lognorm(y ; \mu_N + s^2/2, \sigma_N^2)\end{split}$

Nexp: 100, Correct Guesses: 97, Accuracy: 0.97

Доля правильно угаданных вариантов

ConjugateRevPerUserParams = namedtuple('ConjugateRevPerUserParams', 'a b mu sigma sx')

def posterior_params_rev_per_user(data):
    d_paying = data[data != 0]
    d_paying_total = len(d_paying)
    d_total = len(data)
    a, b = posterior_params_binom(ns=d_paying_total, ntotal=d_total)
    post_pars = initial_params_lognormal(mu=np.log(d_paying)[0], sigma=np.std(np.log(d_paying)), sx=np.std(np.log(d_paying)))
    post_pars = posterior_params_lognormal(d_paying[1:], post_pars)
    return ConjugateRevPerUserParams(a=a, b=b, mu=post_pars.mu, sigma=post_pars.sigma, sx=post_pars.sx)

def posterior_params_binom(ns, ntotal, a_prior=1, b_prior=1):
    a = a_prior + ns
    b = b_prior + ntotal - ns
    return a, b

def rev_per_user_p_dist(params):
    return stats.beta(a=params.a, b=params.b)

def posterior_mean_rev_per_user_rvs(params, nsamples=100_000):
    p_dist = rev_per_user_p_dist(params)
    ps = p_dist.rvs(size=nsamples)
    means_dist = posterior_mean_dist_lognormal(params)
    means = means_dist.rvs(nsamples)
    return ps * means

def exact_rev_per_user_rvs(p, mu, s, nsamples):
    conv = stats.bernoulli.rvs(p=p, size=nsamples)
    rev = stats.lognorm.rvs(s=s, scale=np.exp(mu), size=nsamples)
    return conv * rev

def prob_pb_gt_pa_samples(post_samp_A, post_samp_B):
    if len(post_samp_A) != len(post_samp_B):
        return None
    b_gt_a = np.sum(post_samp_B > post_samp_A)
    return b_gt_a / len(post_samp_A)

nexps = 100
prob_stop = 0.95
n_samp_max = 3_000_000
n_samp_step = 50000
n_post_samp = 50000

A = {'p': 0.1, 'mu': 2, 's': 1}

cmp = pd.DataFrame(columns=['A_pars', 'B_pars', 'A_mean', 'B_mean', 'best_exact', 'exp_samp_size', 'A_exp', 'B_exp', 'best_exp', 'p_best'])
cmp['A_pars'] = [A] * nexps
cmp['B_pars'] = cmp['A_pars'].apply(lambda x: {'p': x['p'] * (1 + stats.uniform.rvs(loc=-0.05, scale=0.1)), 's': x['s'], 'mu': x['mu'] * (1 + stats.uniform.rvs(loc=-0.05, scale=0.1))})
cmp['A_mean'] = cmp['A_pars'].apply(lambda x: x['p'] * stats.lognorm(s=x['s'], scale=np.exp(x['mu'])).mean())
cmp['B_mean'] = cmp['B_pars'].apply(lambda x: x['p'] * stats.lognorm(s=x['s'], scale=np.exp(x['mu'])).mean())
cmp['best_exact'] = cmp.apply(lambda r: 'B' if r['B_mean'] > r['A_mean'] else 'A', axis=1)

for i in range(nexps):
    A_pars = cmp.at[i, 'A_pars']
    B_pars = cmp.at[i, 'B_pars']
    n_samp_total = 0
    dA = np.array([])
    dB = np.array([])
    while n_samp_total < n_samp_max:
        dA = np.append(dA, exact_rev_per_user_rvs(p=A_pars['p'], mu=A_pars['mu'], s=A_pars['s'], nsamples=n_samp_step))
        dB = np.append(dB, exact_rev_per_user_rvs(p=B_pars['p'], mu=B_pars['mu'], s=B_pars['s'], nsamples=n_samp_step))
        n_samp_total += n_samp_step
        post_pars_A = posterior_params_rev_per_user(dA)
        post_pars_B = posterior_params_rev_per_user(dB)
        post_samp_A = posterior_mean_rev_per_user_rvs(post_pars_A)
        post_samp_B = posterior_mean_rev_per_user_rvs(post_pars_B)
        pb_gt_pa = prob_pb_gt_pa_samples(post_samp_A, post_samp_B)
        best_gr = 'B' if pb_gt_pa >= prob_stop else 'A' if (1 - pb_gt_pa) >= prob_stop else None
        if best_gr:
            cmp.at[i, 'A_exp'] = post_samp_A.mean()
            cmp.at[i, 'B_exp'] = post_samp_B.mean()
            cmp.at[i, 'exp_samp_size'] = n_samp_total
            cmp.at[i, 'best_exp'] = best_gr
            cmp.at[i, 'p_best'] = pb_gt_pa
            break
    print(f'done {i}: n_samp {n_samp_total}, best_group {best_gr}, P(b>a) {pb_gt_pa}')

cmp['correct'] = cmp['best_exact'] == cmp['best_exp']
display(cmp.head(10))
cor_guess = np.sum(cmp['correct'])
print(f"Nexp: {nexps}, Correct Guesses: {cor_guess}, Accuracy: {cor_guess / nexps}")

Заказы на посетителя

Посетитель может сделать несколько заказов или не сделать ни одного. Для распределения количества заказов посетителя $P_{заказы}(n)$ , $n \in 0, 1, 2, \dots$ можно ожидать дискретный аналог лог-нормального или степенное распределение вроде распределения Ципфа $P(n ; s) \propto n^{-s}$ [ZipfDist, SciPyZipfian, SciPyZipf]. Важно точно промоделировать вероятности небольшого количества заказов. Распределение Ципфа в качестве функции правдоподобия может этого не позволить. Более гибкой моделью будет набор вероятностей p_i сделать заказов. Пусть максимальное количество заказов ограничено, n_i - количество пользователей с $i=0, 1, 2, \dots, N$ заказами. Функция правдоподобия задается мультиномиальным распределением $P(\mathcal{D} | \mathcal{H}) = Mult(n_0, \dots, n_N | p_0, \dots, p_N)$ [MultiDist, SciPyMult]. Сопряженным априорным распределением будет распределение Дирихле $P(\mathcal{H}) = Dir \left( p_{0}, \dots, p_{N}; \alpha_{0}, \dots, \alpha_{N} \right)$ [DirDist, SciPyDir]. В апостериорном распределении обновленные параметры $\alpha_i + n_i$ . Маржинальными распределениями каждого p_i будут бета-распределения, что согласуется с интерпретацией p_i как конверсией в заказов.

$\begin{split}P(\mathcal{D} | \mathcal{H}) & = Mult(n_0, \dots, n_N | p_0, \dots, p_N) = \frac{(n_0 + \dots + n_N)!}{n_{0}! \dots n_{N}!} p_{0}^{n_{0}} \dots p_{N}^{n_{N}} \\P(\mathcal{H}) & = Dir \left( p_{0}, \dots, p_{N}; \alpha_{0}, \dots, \alpha_{N} \right) = \dfrac{1}{B( \alpha_{0}, \dots, \alpha_{N} )} \prod_{i=0}^{N} p_{i}^{\alpha_{i}-1}, \\ &\sum_{i=0}^{N} p_i = 1,\qquad p_i \in [0, 1], \qquad B(\alpha_{0}, \dots, \alpha_{N}) = \frac{\prod \limits_{i=0}^{N} \Gamma( \alpha_{i} )}{\Gamma \left( \sum \limits_{i=0}^{N} \alpha_{i} \right)}\\P(\mathcal{H} | \mathcal{D}) & \propto Mult(n_0, \dots, n_N | p_0, \dots, p_N) Dir \left( p_{0}, \dots, p_{N}; \alpha_{0}, \dots, \alpha_{N} \right)\\& \propto p_{0}^{n_{0}} \dots p_{N}^{n_{N}} \prod _{i=0}^{N} p_{i}^{\alpha_{i}-1}\\& \propto \prod_{i=0}^{N} p_{i}^{n_{i} + \alpha_{i} - 1}\\& =Dir \left( p_{0}, \dots, p_{N}; \alpha_{0} + n_0, \dots, \alpha_{N} + n_N \right)\\P(p_i | \mathcal{D} ) & = \int dp_0 \dots dp_{i-1}dp_{i+i} \dots dp_N P(\mathcal{H} | \mathcal{D}) \\ &= Beta( p_i; \alpha_i + n_i, \sum_{k=0}^{N} (\alpha_k + n_k) - \alpha_i - n_i )\end{split}$

Для примера оценки параметров задается распределение Ципфа с параметрами s, Nmax. Из него делается выборка, по выборке строится апострериорное распределение. Также считаются распределения конверсий в заказов p_i . На графике приведены исходное распределение, выборка, апостериорное предиктивное распределение , оценки и 95%-интервалы p_i . Для большей части точные значения лежат внутри интервалов p_i .

Оценка распределения количества заказов

def initial_params_dir(N):
    return np.ones(N)

def posterior_params_dir(data, initial_pars):
    u, c = np.unique(data, return_counts=True)
    post_pars = np.copy(initial_pars)
    for k, v in zip(u, c):
        post_pars[k] = post_pars[k] + v
    return post_pars

def posterior_dist_dir(params):
    return stats.dirichlet(alpha=params)

def posterior_nords_dir_rvs(params, nsamp):
    nords = np.empty(nsamp)
    d = posterior_dist_dir(params)
    probs = d.rvs(size=nsamp)
    for i, p in enumerate(probs):
        nords[i] = np.argmax(stats.multinomial.rvs(n=1, p=p))
    return nords

def marginal_pi_dist_dir(i, params):
    return stats.beta(a=params[i], b=np.sum(params) - params[i])

def posterior_pi_mean_95pdi(i, params):
    p = marginal_pi_dist_dir(i, params)
    m = p.mean()
    lower = p.ppf(0.025)
    upper = p.ppf(0.975)
    return m, lower, upper

Nmax = 30
s = 1.5
nsample = 1000

Npars = Nmax + 1
exact_dist = stats.zipfian(a=s, n=Npars, loc=-1)
data = exact_dist.rvs(nsample)
pars = initial_params_dir(Npars)
pars = posterior_params_dir(data, pars)
post_samp = posterior_nords_dir_rvs(pars, 100000)
pi = [posterior_pi_mean_95pdi(i, pars) for i in range(Npars)]

x = np.arange(0, Npars+1)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=exact_dist.pmf(x), name='Точное распределение Ципфа', 
                         line_color='black'))
fig.add_trace(go.Histogram(x=data, histnorm='probability', name='Выборка', nbinsx=round(Nmax*2),
                         marker_color='black'))
fig.add_trace(go.Histogram(x=post_samp, histnorm='probability', name='$\mbox{Апостериорные } n_i$', 
                         marker_color='black', opacity=0.2, nbinsx=round(Nmax*2)))
fig.add_trace(go.Scatter(x=x, 
                         y=[p[0] for p in pi],
                         error_y=dict(type='data', symmetric=False, array=[p[2] - p[0] for p in pi], arrayminus=[p[0] - p[1] for p in pi]), 
                         name='$\mbox{Оценки } p_i$',
                         mode='markers',
                         line_color='red',
                         opacity=0.8))
fig.update_layout(title='Заказы на посетителя',
                  xaxis_title='$i$',
                  yaxis_title='Вероятность',
                  xaxis_range=[-1, Nmax+1],
                  hovermode="x",
                  barmode="group",
                  height=550)
fig.show()
fig.write_image("./figs/ch6_postdist.png", scale=2)

Распределение заказов на посетителя. Точные конверсии в i заказов лежат внутри оцененных интервалов.

Распределение количества заказов позволяет оценить среднее количество заказов $E[n] = \sum i p_i$ , конверсию в заказ 1-p_0 , конверсию в 2 и более заказов 1-p_0-p_1 . Ниже приведен пример сравнения среднего количества заказов E[n] . Задается 2 распределения Ципфа, параметр s группы Б на 5% меньше группы А. На первом графике приведены точные распределения, точные средние и оценки p_i . Для большей части точные значения лежат внутри интервалов p_i . На втором графике приведено апостериорное распределение среднего количества заказов. Для выбранных параметров вероятность, что среднее группы Б выше А $P(E[n]_B > E[n]_A) = 98.2\%$ .

Сравнение групп

def posterior_nords_mean_rvs(params, nsample):
    ns = np.arange(len(params))
    probs = stats.dirichlet.rvs(alpha=params, size=nsample)
    means = np.sum(ns * probs, axis=1)
    return means

def prob_pb_gt_pa_samples(post_samp_A, post_samp_B):
    if len(post_samp_A) != len(post_samp_B):
        return None
    b_gt_a = np.sum(post_samp_B > post_samp_A)
    return b_gt_a / len(post_samp_A)

nsample = 3000
Nmax = 30
Npars = Nmax + 1

post_samp_len = 100000
A, B = {}, {}
s = 1.5
A['dist_pars'] = {'s': s}
B['dist_pars'] = {'s': s * 0.95}
for g in [A, B]:
    g['exact_dist'] = stats.zipfian(a=g['dist_pars']['s'], n=Npars, loc=-1)
    g['data'] = g['exact_dist'].rvs(nsample)
    g['post_pars'] = initial_params_dir(Npars)
    g['post_pars'] = posterior_params_dir(g['data'], g['post_pars'])
    g['post_nords'] = posterior_nords_dir_rvs(g['post_pars'], post_samp_len)
    g['post_means'] = posterior_nords_mean_rvs(g['post_pars'], post_samp_len)
    g['pi'] = [posterior_pi_mean_95pdi(i, g['post_pars']) for i in range(Npars)]

x = np.arange(0, Npars)
fig = go.Figure()
fig.add_trace(go.Bar(x=x, y=A['exact_dist'].pmf(x), name='Точное распределение A',
                        marker_color='black', opacity=0.2))
fig.add_trace(go.Bar(x=x, y=B['exact_dist'].pmf(x), name='Точное распределение Б',
                        marker_color='black', opacity=0.8))
fig.add_trace(go.Scatter(x=[A['exact_dist'].mean(), A['exact_dist'].mean()], 
                         y=[0, np.max(A['exact_dist'].pmf(x))*1.1],
                         name='Точное среднее A', 
                         mode='lines', line_dash='dash',
                         line_color='black', opacity=0.3))
fig.add_trace(go.Scatter(x=[B['exact_dist'].mean(), B['exact_dist'].mean()], 
                         y=[0, np.max(B['exact_dist'].pmf(x))*1.1],
                         name='Точное среднее Б', 
                         mode='lines', line_dash='dash',
                         line_color='black'))
fig.add_trace(go.Scatter(x=x - 0.1, 
                         y=[p[0] for p in A['pi']],
                         error_y=dict(type='data', symmetric=False, array=[p[2] - p[0] for p in A['pi']], arrayminus=[p[0] - p[1] for p in A['pi']]), 
                         name='$p_i, A$',
                         line_color='black', opacity=0.3,
                         mode='markers'
                    ))
fig.add_trace(go.Scatter(x=x + 0.1, 
                         y=[p[0] for p in B['pi']],
                         error_y=dict(type='data', symmetric=False, array=[p[2] - p[0] for p in B['pi']], arrayminus=[p[0] - p[1] for p in B['pi']]), 
                         name='$p_i, Б$',
                         line_color='black',
                         mode='markers'))
fig.update_layout(title='Заказы на посетителя',
                  xaxis_title='$n$',
                  yaxis_title='Вероятность',
                  xaxis_range=[-1, Npars+1-20],
                  hovermode="x",
                  barmode="group",
                  height=550)
fig.show()

x = np.arange(0, Npars)
fig = go.Figure()
fig.add_trace(go.Scatter(x=[A['exact_dist'].mean(), A['exact_dist'].mean()], 
                         y=[0, np.max(A['exact_dist'].pmf(x))*1.1],
                         name='Точное среднее A', 
                         mode='lines', line_dash='dash',
                         line_color='black', opacity=0.3))
fig.add_trace(go.Scatter(x=[B['exact_dist'].mean(), B['exact_dist'].mean()], 
                         y=[0, np.max(B['exact_dist'].pmf(x))*1.1],
                         name='Точное среднее Б', 
                         mode='lines', line_dash='dash',
                         line_color='black'))
fig.add_trace(go.Histogram(x=A['post_means'], histnorm='probability', name='Среднее n, A', 
                           marker_color='black', opacity=0.3, nbinsx=round(Nmax*2)))
fig.add_trace(go.Histogram(x=B['post_means'], histnorm='probability', name='Среднее n, Б', 
                           marker_color='black', nbinsx=round(Nmax*2)))
fig.update_layout(title='Среднее количество заказов',
                  xaxis_title='$n$',
                  yaxis_title='Вероятность',
                  xaxis_range=[-1, Npars+1-20],
                  hovermode="x",
                  barmode="group",
                  height=550)
fig.show()

print(f"P(E[n]_B > E[n]_A): {prob_pb_gt_pa_samples(A['post_means'], B['post_means'])}")

Точные распределения, точные средние количества заказов и оценки конверсий.

Оценки среднего количества заказов. Среднее Б выше А с вероятностью 98.2%.

Количество правильно угаданных "лучших" групп проверяется в nexps экспериментах. В группе A количество заказов на пользователя задается распределением Ципфа с параметром s, в группе Б параметр в пределах $\pm 5\%$ от A. Группы сравниваются по среднему количеству заказов. В экспериментах в выборки пошагово добавляется n_samp_step точек, считаются параметры апостериорных распределений и вероятность P(E[n]_B > E[n]_A) . Эксперимент останавливается при достижении вероятности среднего одной из групп больше другой prob_stop или наборе максимального количества точек n_samp_max. Доля правильно угаданных групп 0.94 близка ожидаемой prob_stop = 0.95.

Nexp: 100, Correct Guesses: 94, Accuracy: 0.94

Доля правильно угаданных вариантов

cmp = pd.DataFrame(columns=['A', 'B', 'best_exact', 'exp_samp_size', 'A_exp', 'B_exp', 'best_exp', 'p_best'])

s = 1.5
Nmax = 30
Npars = Nmax + 1
nexps = 100
cmp['A'] = [s] * nexps
cmp['B'] = s * (1 + stats.uniform.rvs(loc=-0.05, scale=0.1, size=nexps))

n_samp_max = 200000
n_samp_step = 5000

prob_stop = 0.95
for i in range(nexps):
    s_a = cmp.at[i, 'A']
    s_b = cmp.at[i, 'B']
    exact_dist_a = stats.zipfian(a=s_a, n=Npars, loc=-1)
    exact_dist_b = stats.zipfian(a=s_b, n=Npars, loc=-1)
    cmp.at[i, 'best_exact'] = 'A' if exact_dist_a.mean() > exact_dist_b.mean() else 'B'
    n_samp_total = 0
    pars_a = initial_params_dir(Npars)
    pars_b = initial_params_dir(Npars)
    while n_samp_total < n_samp_max:
        data_a = exact_dist_a.rvs(n_samp_step)
        data_b = exact_dist_b.rvs(n_samp_step)
        n_samp_total += n_samp_step
        pars_a = posterior_params_dir(data_a, pars_a)
        pars_b = posterior_params_dir(data_b, pars_b)
        post_samp_len = 10000
        post_means_a = posterior_nords_mean_rvs(pars_a, post_samp_len)
        post_means_b = posterior_nords_mean_rvs(pars_b, post_samp_len)
        pb_gt_pa = prob_pb_gt_pa_samples(post_means_a, post_means_b)
        best_gr = 'B' if pb_gt_pa >= prob_stop else 'A' if (1 - pb_gt_pa) >= prob_stop else None
        if best_gr:
            cmp.at[i, 'A_exp'] = post_means_a.mean()
            cmp.at[i, 'B_exp'] = post_means_b.mean()
            cmp.at[i, 'exp_samp_size'] = n_samp_total
            cmp.at[i, 'best_exp'] = best_gr
            cmp.at[i, 'p_best'] = pb_gt_pa
            break
    print(f'done {i}: nsamp {n_samp_total}, best_gr {best_gr}, P(B>A) {pb_gt_pa}')

cmp['correct'] = cmp['best_exact'] == cmp['best_exp']
display(cmp.head(10))
cor_guess = np.sum(cmp['correct'])
print(f"Nexp: {nexps}, Correct Guesses: {cor_guess}, Accuracy: {cor_guess / nexps}")

Заключение

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

Ссылки

[Apx] - Bayesian AB Testing, Appendices, GitHub.
[BaseFal] - Base Rate Fallacy, Wikipedia.
[BernProc] - Bernoulli Process, Wikipedia.
[BerryEsseenTheorem] - Berry-Esseen Theorem, Wikipedia.
[BetaDist] - Beta Distribution, Wikipedia.
[BinomDist] - Binomial Distribution, Wikipedia.
[CLT] - Central Limit Theorem, Wikipedia.
[CausalDAG] - Causal Graph, Wikipedia.
[ConjPrior] - Conjugate Prior, Wikipedia.
[DirDist] - Dirichlet Distribution, Wikipedia.
[GammaDist] - Gamma Distribution, Wikipedia.
[LognormDist] - Log-normal Distribution, Wikipedia.
[LomaxDist] - Lomax Distribution, Wikipedia.
[MicroExp] - R. Kohavi et al, Online Experimentation at Microsoft.
[MultiDist] - Multinomial Distribution, Wikipedia.
[NormDist] - Normal Distribution, Wikipedia.
[ParetoDist] - Pareto Distribution, Wikipedia.
[ProbConv] - Convolution of Probability Distributions, Wikipedia.
[RandVarsConv] - Convergence of Random Variables, Wikipedia.
[SGBS] - B. Lambert, A Student’s Guide to Bayesian Statistics (Textbook, Student Resources).
[SR] - R. McElreath, Statistical Rethinking: A Bayesian Course with Examples in R and STAN (Textbook, Video Lectures, Course Materials).
[SciPyBern] - scipy.stats.bernoulli, SciPy Reference.
[SciPyBeta] - scipy.stats.beta, SciPy Reference.
[SciPyBinom] - scipy.stats.binom, SciPy Reference.
[SciPyDir] - scipy.stats.dirichlet, SciPy Reference.
[SciPyGamma] - scipy.stats.gamma, SciPy Reference.
[SciPyLognorm] - scipy.stats.lognorm, SciPy Reference.
[SciPyLomax] - scipy.stats.lomax, SciPy Reference.
[SciPyMult] - scipy.stats.multinomial, SciPy Reference.
[SciPyNorm] - scipy.stats.norm, SciPy Reference.
[SciPyPareto] - scipy.stats.pareto, SciPy Reference.
[SciPyPareto] - scipy.stats.pareto, SciPy Reference.
[SciPyZipf] - scipy.stats.zipf, SciPy Reference.
[SciPyZipfian] - scipy.stats.zipfian, SciPy Reference.
[SubjProb] - Probability Interpretations, Wikipedia.
[ZipfDist] - Zipf's Law, Wikipedia.

Tags:

Hubs: