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

Question

Я отправил питона ноутбук здесь: 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 для определения фи (х), используя функцию ошибки, подобно логистический пример регрессионного из хранилища PyMC ,

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

Answer 1

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

Единственное отличие, которое я имею, это использовать функцию тензора 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))))