Пробит-регрессия с использованием PyMC 3

Я разместил блокнот Python здесь: http://nbviewer.ipython.org/gist/awellis/9067358

Я пытаюсь создать модель пробит-регрессии с использованием PyMC 3, используя сгенерированные данные для восстановления известных параметров (см. Блокнот). Оценка точки пересечения почти в порядке, но оценка наклона далеко не соответствует действительности.

Моя модель выглядит так:

with pm.Model() as model:

    # priors
    alpha = pm.Normal('alpha', mu=0, tau=0.001)
    beta = pm.Normal('beta', mu=0, tau=0.001)

    # linear predictor
    theta_p = (alpha + beta * x)

    # logic transform (just for comparison - this seems to work ok)
#     def invlogit(x):
#         import theano.tensor as t
#         return t.exp(x) / (1 + t.exp(x))
#     theta = invlogit(theta_p)


    # Probit transform: this doesn't work
    def phi(x):
        import theano.tensor as t
        return 0.5 * (1 + t.erf(x / t.sqr(2)))
    theta = phi(theta_p)


    # likelihood
    y = pm.Bernoulli('y', p=theta, observed=y)

with model:
    # Inference
    start = pm.find_MAP() # Find starting value by optimization

    print("MAP found:")
    print("alpha:", start['alpha'])
    print("beta:", start['beta'])

    print("Compare with true values:")
    print("true_alpha", true_alpha)
    print("true_beta", true_beta)

with model:
        step = pm.NUTS()
        trace = pm.sample(2000,
                           step, 
                           start=start, 
                           progressbar=True) # draw posterior samples

Единственный способ, которым это работает, — это использовать Theano для определения phi(x) с использованием функции ошибки, аналогично примеру логистической регрессии из репозитория PyMC.

Может кто-то указать мне верное направление? Есть ли лучший/более простой способ сделать это?


python bayesian pymc
person Andrew Ellis    schedule 18.02.2014    source источник

Ответы (1)


arrow_upward
3
arrow_downward

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

Единственное отличие, которое у меня есть, это то, что я использовал функцию тензора sqrt(). Может просто опечатка с вашей стороны?

import theano.tensor as tsr

def probit_phi(x):
    """ Probit transform assuming 0 mean and 1 sd """
    mu = 0
    sd = 1
    return 0.5 * (1 + tsr.erf((x - mu) / (sd * tsr.sqrt(2))))
person jonsedar    schedule 06.10.2014