У нас часто есть данные, которые представляют собой нулевое завышенное распределение, в котором есть лишние нули. Это очень распространенный сценарий во многих отраслях. Например, у нас есть этот иллюстрированный пример ниже, который имеет ~ 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))

Вывод этого теста:

  1. Вариантные удары Контрольная группа со средней разницей средних значений колебалась от 15,0 до 17,0 на одного испытуемого, рост от 12% до 13% по сравнению с исходным уровнем
  2. Мы можем спрогнозировать инкремент, применив № 1 к целевому будущему исследуемому населению.
  3. Учитывая фактор Байеса > 1, мы можем завершить этот тест с убедительными доказательствами.