PYMC3 Моделирование полиномиального выбора

Question

Обновление

Использование pm.Multinomial с n=1 и shape=11 вместо pm.Categorical работает, однако мне также нужно указать init='adapt_diag' при выборке, или это приведетв плохой энергии. Есть идеи по этому поводу? Я до сих пор не понимаю, почему Категориальный не работает.

Исходное сообщение

Я пытаюсь сделать байесовский вывод на смоделированных данных, т.е. я знаю процесс генерации данныхи истинные параметры.

Вы можете думать о процессе генерирования данных следующим образом:

- У клиента, например, на странице оформления заказа на Amazon, есть 11 вариантов выбора

- Вариант 1-9: Заказать и получить доставку в одном из 9 предварительно определенных временных окон

- Выбор 10: Заказывать, но не выбирать предварительно определенные временные окна

- Выбор 11: Оставитьоформить заказ и не заказывать

- всегда делается один из приведенных выше вариантов выбора

- приведенные ниже данные показывают выбор, сделанный для каждого клиента, каждая строка является клиентом, каждый столбец - фиктивным(например, клиент 1 выбирает заказ во временном окне 9)

- Однако не все варианты доступны всегда

- Варианты 10 и 11 всегда доступны

- Вариант 1-9 может быть или не быть доступным

- Если выбор недоступен, клиент выбирает из оставшихся вариантов с вероятностью, которая пропорциональна скрытому вектору предпочтений

. Данные генерируются следующим образом:

import numpy as np
import pandas as pd
import pymc3 as pm
from theano import tensor as T
from scipy.stats import bernoulli, logistic

# Number of customers
N = 3000

# Latent preferences for choices 1-10: If a choice is not available, 
# the probability of the other choice increase proportionally to preference_lat 
# preference_lat  shall later be estimated from the data and does not work yet in 
# my pymc3 implementation
preference_lat = np.array([0.25, 0.125, 0.125, 0.10, 0.10, 0.10, 0.10, 0.05, 0.025, 0.025])
# Probability that choices 1-10 are available. Probability of Choice 10 is set to 1, so
# it is always available. The estimation of this in my pymc3 implementation works 
p_slots_avail = np.array([0.1, 0.1, 0.1, 0.5, 0.5, 0.5, 0.8, 0.8, 0.8, 1])
# Draw from choice space
slots_available = bernoulli.rvs(p_slots_avail, size=(N, 10))
# Number of available choices per customer
no_slots = slots_available.sum(axis=1)

# Probabilty of choice 12: When the choice space is larger, the probability #of 
# actually ordering increases (the probability of not-ordering decreases). The 
# estimation of a and b in my pymc3 implementation through logistic regression
# works well
a = -2
b = 0.8
p_conv = logistic.cdf(a + b*no_slots)

# Keep only preferences for available choices
slot_prob = slots_available * preference_lat
# Normalize preferences to probability of choice 1 to 10 given that customer orders
slot_prob_normalized = slot_prob / (slots_available*preference_lat).sum(axis=1)[..., np.newaxis]
# Add the probability of not ordering (logistic output) and renormalize such that 
# all probabilities sum to 1
slot_prob_norm2 = p_conv[..., np.newaxis] * slot_prob_normalized
data_prob = np.concatenate((slot_prob_norm2, 1 - p_conv[..., np.newaxis]), axis=1)

# Draw from choice space. Each customer (row) faces has a different 
# choice distribution as the available choices different between each customer
# (according to p_slots_avail)
def mnomial(slc):
    return np.random.multinomial(1, slc)
data_simulated = np.apply_along_axis(mnomial, 1, data_prob)

data_simulated_df = pd.DataFrame(data_simulated)

def get_choice(row):
    for c in data_simulated_df.columns:
        if row[c]==1:
            return c

data_simulated_rev = np.asarray(data_simulated_df.apply(get_choice, axis=1))

Я пытаюсь сделать вывод, используя pymc3, как показано ниже:

with pm.Model() as model:
    # This part works as intended, the posteriors are around the true values of the 
    # choice availability distributions
    pm_p_slots_avail = pm.Uniform('pm_p_slots_avail', shape=10, lower=0, upper=1)
    pm_slots_available = pm.Bernoulli('pm_slots_available', p=pm_p_slots_avail, shape=10, 
                                       observed=slots_available)

    # This logistic regression part works as intended, the posteriors of pm_a and pm_b are 
    # around the true values of a and b
    pm_no_slots = pm.Deterministic('pm_no_slots', pm_slots_available.sum(axis=1))
    pm_a = pm.Uniform('pm_a', lower=-10, upper=10)
    pm_b = pm.Uniform('pm_b', lower=-10, upper=10)
    pm_p_conv = pm.Deterministic('pm_p_conv', pm.math.sigmoid(pm_a + pm_b*pm_no_slots))
    pm_order = pm.Bernoulli('pm_order', 
                            p=pm_p_conv,
                            observed=1-data_simulated[:, 10])

    # This part does not work / sampling for pm_choice is extremely slow
    pm_preference_lat = pm.Dirichlet('pm_preference_lat', a=np.ones(10), shape=10)
    pm_slot_prob = pm_slots_available * pm_preference_lat
    pm_slot_prob_normalized = pm_slot_prob / (pm_slots_available*pm_preference_lat).sum(axis=1)[..., np.newaxis] 
    pm_slot_prob_norm2 = pm_p_conv[..., np.newaxis] * pm_slot_prob_normalized
    pm_data_prob = T.concatenate((pm_slot_prob_norm2, 1 - pm_p_conv[..., np.newaxis]), axis=1)
    pm_choice = pm.Categorical('pm_choice', 
                               p=pm_data_prob,
                               observed=data_simulated_rev)

Однако последняя часть, которую мне нужно сделать, pm_preference_lat, работает не так, как задумано. Когда я пробую

with model:
    trace = sample(25000, cores=1)

, я получаю "плохую начальную энергию". Если я передам p=pm_data_prob + 0.01 вместо p=pm_data_prob в pm.Categorical -Call выше, семплеры будут работать. Однако, выборка очень медленная.

Есть ли у меня неправильное понимание того, как работает pm.Categorical? Модель как-то неверно определена?