Probit regression using PyMC 3

Question

I have posted an python notebook here: http://nbviewer.ipython.org/gist/awellis/9067358

I am trying create a probit regression model using PyMC 3, using generated data to recover the known parameters (see notebook). The estimate for the intercept is just about ok, but the slope estimate is way off the mark.

My model looks like this:

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

The only way it seems to work is to use Theano to define phi(x), using the error function, similarly to the logistic regression example from the PyMC repository.

Can anyone point me in the right direction? Is there a better/easier way of doing this?

Answer 1

This may be long after the horse has bolted, but I've just tried implementing this myself for a simple hierarchical Binomial model and found results comparable to a logit function.

The only difference I have is I've used the tensor sqrt() function. Possibly just a typo on your part?

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))))