title: Probit regression using PyMC 3
body: <p>I have posted an python notebook here: <a href="http://nbviewer.ipython.org/gist/awellis/9067358">http://nbviewer.ipython.org/gist/awellis/9067358</a></p>
  
  <p>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.</p>
  
  <p>My model looks like this:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>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. </p>
  
  <p>Can anyone point me in the right direction? Is there a better/easier way of doing this?</p>