У нас часто есть данные, которые представляют собой нулевое завышенное распределение, в котором есть лишние нули. Это очень распространенный сценарий во многих отраслях. Например, у нас есть этот иллюстрированный пример ниже, который имеет ~ 30% нулей в данных.
С точки зрения бизнеса, мы хотели бы завершить эту оценку тремя целями:
A. Можем ли мы сделать вывод, превосходит ли вариантная группа базовую контрольную группу?
B. Каков диапазон повышения и прироста по сравнению с группой вариантов по сравнению с исходным уровнем?
C. Должны ли мы продолжить или завершить тест?
Сначала построим распределение данных. Мы можем проверить форму и масштаб данных. Совершенно очевидно, что мы не можем использовать стандартный метод статистики для оценки теста, поскольку распределение представляет собой длинный хвост с избыточным нулем. Мы также хотели бы использовать байесовский метод с его гибкостью, которая поможет нам ответить на поставленные выше вопросы.
Затем мы используем Fitter, чтобы помочь нам выбрать 5 лучших кандидатов, которых мы можем включить в наш эксперимент. В этом случае мы выбираем Inverse Gaussian. Вы можете узнать больше информации о Fitter здесь: https://github.com/cokelaer/fitter
Преобразование из Scipy invgauss в TFD InverseGaussian Priors
Прежде чем мы создадим априорные данные TFD, мы используем scipy stats MLE для получения априорных значений. Scipy invgauss возвращает форму, местоположение и масштаб. Нам нужно преобразовать три априора в спецификацию с двумя априорами, используемую в TFD InverseGaussian.
_sh, _lo, _sc = stats.invgauss.fit(sample_data)
mu = _sh * _sc + _lo
con = _sc
Целевое совместное распространение
С приведенным выше преобразованием у нас есть локальные и концентрационные априорные значения для установки слабых априорных значений для InverseGaussian. Нам также нужно добавить 1-й априор, используемый в качестве вероятности смеси, чтобы представить его как нулевое завышенное обратное значение Гаусса.
tgt_zi_model = tfd.JointDistributionSequential([
tfd.Sample(tfd.Uniform(0, 1.), 1), # mixture prob
tfd.Uniform(mu-100, mu+100), # loc weak prior
tfd.Uniform(con-1500, con+1500), # concentration weak prior
tfd.Sample(lambda co, lo, _: tfd.InverseGaussian(loc=lo, concentration=co)), sample_data.shape[0]])
Чтобы вычислить нулевую завышенную логарифмическую вероятность, мы используем двухчастный подход модели. Первая часть — это Бернулли для оптимизации нуля/не нуля. Вторая часть — непрерывная логарифмическая вероятность обратного Гаусса. Затем мы суммируем обе логарифмические вероятности как совместную вероятность.
sample_data_zmask = tf.cast(sample_data == 0, 'float32')
sample_data_nzmask = sample_data[sample_data > 0]
# ZI two part model
def target_log_prob_fn(*params):
zero_lp = tfd.Bernoulli(probs=params[0]).log_prob(sample_data_zmask)
t_zimodel = init_zimodel(sample_data_nzmask.shape[0])
nzero_lp = t_zimodel.log_prob(*params + (sample_data_nzmask, ))
return tf.reduce_sum(zero_lp) + nzero_lp
Пример инициализированного состояния MCMC
Чтобы запустить MCMC, мы выбираем априорные значения начального состояния, которые оптимизируют и находят его траекторию с помощью интегратора чехарды.
_init_state_weak = tgt_zi_model.sample(1000)[:-1]
_init_state_weak = tf.unstack(tf.reduce_mean(_init_state_weak, axis=1))
MCMC одинарная цепь
chain_state, sampler_stat = mcmc.sample_chain(
num_results=num_steps,
current_state=_init_state_weak,
num_burnin_steps=500,
kernel=mcmc.DualAveragingStepSizeAdaptation(
inner_kernel=mcmc.TransformedTransitionKernel(
inner_kernel=mcmc.HamiltonianMonteCarlo(
target_log_prob_fn=target_log_prob_fn, # MLE optimize
step_size=0.1, # leapfrog step size
num_leapfrog_steps=4, # leapfrog inner loop
state_gradients_are_stopped=True),
bijector=[tfb.Sigmoid(), tfb.Softplus(), tfb.Softplus()]),
num_adaptation_steps=40))
Трассировка MCMC для контрольной группы и группы вариантов сзади
Чтобы вычислить распределение средней разности, мы можем получить апостериорное значение из трассировки chain_state MCMC. Параметр loc представляет среднее апостериорное значение. Взяв разницу между вариантной группой апостериорного мю и контрольной группой апостериорного мю, мы можем рассчитать распределение средних различий, а также распределение подъема, показанное ниже.
Имея три апостериорных значения, мы можем вывести байесовский коэффициент между тестовой и контрольной группой. Мы вычисляем P(вариантные данные | H1 = вариантные апостериорные значения) \div P(вариантные данные | H0 = контрольные апостериорные значения). С байесовским коэффициентом > 1 мы можем считать этот эксперимент завершенным с убедительными доказательствами.
p_bi_c = tf.reduce_mean(tfd.Bernoulli(tf.reduce_mean(bi_c_trace, axis=1)).prob(outcome_bin_t)) p_bi_t = tf.reduce_mean(tfd.Bernoulli(tf.reduce_mean(bi_t_trace, axis=1)).prob(outcome_bin_t)) p_cont_c = tf.reduce_mean(tfd.InverseGaussian(loc=tf.reduce_mean(cont_c_trace[0], axis=0), concentration=tf.reduce_mean(cont_c_trace[1], axis=0)).prob(outcome_t)) p_cont_t = tf.reduce_mean(tfd.InverseGaussian(loc=tf.reduce_mean(cont_t_trace[0], axis=0), concentration=tf.reduce_mean(cont_t_trace[1], axis=0)).prob(outcome_t))
Вывод этого теста:
- Вариантные удары Контрольная группа со средней разницей средних значений колебалась от 15,0 до 17,0 на одного испытуемого, рост от 12% до 13% по сравнению с исходным уровнем
- Мы можем спрогнозировать инкремент, применив № 1 к целевому будущему исследуемому населению.
- Учитывая фактор Байеса > 1, мы можем завершить этот тест с убедительными доказательствами.