Данные
Во-первых, обратите внимание, что общая вероятность повторных испытаний Бернулли является в точности биномиальной вероятностью, поэтому нет необходимости расширять отдельные исследования в ваших данных.Я бы также предложил использовать DataFrame Pandas для управления вашими данными - это помогает поддерживать порядок:
import pandas as pd
df = pd.DataFrame({
'n': [0, 5, 10, 15, 20, 25],
'trials': [120, 111, 78, 144, 280, 55],
'successes': [1, 2, 1, 3, 7, 1]
})
Решение
Это поможет упростить модель, но решение на самом деле заключается вдобавьте аргумент shape к случайной переменной p, чтобы PyMC3 знал, как интерпретировать одномерные параметры.Дело в том, что вам нужно разное p распределение для каждого имеющегося у вас n случая, поэтому здесь нет ничего концептуально неправильного.
with pm.Model() as model:
# conversion rate hyperparameters
alpha = pm.Uniform("alpha_n", 5, 13)
beta = pm.Uniform("beta_n", 1000, 1400)
# switchpoint prior
n_sat = pm.Gamma("n_sat", alpha=20, beta=2, testval=10)
a_gamma = pm.Gamma("a_gamma", alpha=18, beta=15)
b_gamma = pm.Gamma("b_gamma", alpha=18, beta=27)
# NB: I removed pm.Deterministic b/c (a|b)_slope[0] is constant
# and this causes issues when using ArViZ
a_slope = 1 + (df.n.values/n_sat)*(a_gamma-1)
b_slope = 1 + (df.n.values/n_sat)*(b_gamma-1)
a = pm.math.switch(df.n.values >= n_sat, a_gamma, a_slope)
b = pm.math.switch(df.n.values >= n_sat, b_gamma, b_slope)
# conversion rates
p = pm.Beta("p", alpha=alpha*a, beta=beta*b, shape=len(df.n))
# observations
pm.Binomial("observed", n=df.trials, p=p, observed=df.successes)
trace = pm.sample(5000, tune=10000)
Этот пример прекрасно
и дает разумные интервалы по коэффициентам конверсии
, но тот факт, что потомки для alpha_n и beta_n идут прямо до вашего предыдущегограницы немного касаются:
Я думаю, что причина этого заключается в том, что для каждого условия вы делаете только 55-280 испытаний, которые, если условия были независимыми (худшееcase), сопряженность говорит нам, что ваши бета-гиперпараметры должны быть в этом диапазонеПоскольку вы выполняете регрессию, тогда в наилучшем сценарии обмена информацией между испытаниями ваши гиперпараметры будут помещаться в диапазон суммы испытаний (788), но это верхний предел.Поскольку вы находитесь за пределами этого диапазона, проблема заключается в том, что вы заставляете модель быть более точной в своих оценках, чем у вас есть фактические данные для подтверждения.Тем не менее, это может быть оправдано, если предварительное обоснование основано на сильных независимых доказательствах.
В противном случае, я бы предложил расширить диапазоны для тех априоров, которые влияют на окончательные alpha*a и beta*b чисел (их суммы должны быть близки к вашим пробным счетам в задней части).
Альтернативная модель
Я бы, вероятно, сделал что-то по следующим направлениям, которые яУ Think более прозрачная параметризация, хотя она не полностью идентична вашей модели:
with pm.Model() as model_br_sp:
# regression coefficients
alpha = pm.Normal("alpha", mu=0, sd=1)
beta = pm.Normal("beta", mu=0, sd=1)
# saturation parameters
saturation_point = pm.Gamma("saturation_point", alpha=20, beta=2)
max_success_rate = pm.Beta("max_success_rate", 1, 9)
# probability of conversion
success_rate = pm.Deterministic("success_rate",
pm.math.switch(df.n.values > saturation_point,
max_success_rate,
max_success_rate*pm.math.sigmoid(alpha + beta*df.n)))
# observations
pm.Binomial("successes", n=df.trials, p=success_rate, observed=df.successes)
trace_br_sp = pm.sample(draws=5000, tune=10000)
Здесь мы отображаем пространство предиктора в пространство вероятности через сигмовидную форму, которая максимизируется с максимальной вероятностью успеха.Априорная точка на точке насыщения идентична вашей, а максимальная вероятность успеха слабо информативна (бета [1,9] - хотя я скажу, что она также работает на флете до).Это также хорошо,
и дает аналогичные интервалы (хотя точка переключения, кажется, доминирует больше):
Мыможно сравнить две модели и увидеть, что нет существенной разницы в их объяснительной силе:
import arviz as az
model_compare = az.compare({'Binomial Regression w/ Switchpoint': trace_br_sp,
'Original Model': trace})
az.plot_compare(model_compare)
