Probit

from scipy.stats import norm
from pymc import Normal, Lambda, Bernoulli
from numpy import mean, std

# Data
grade = (0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1)
sat = (525, 533, 545, 582, 581, 576, 572, 609, 559, 543, 576, 525, 574, 582, 574, 471, 595, 557, 557, 584, 599, 517, 649, 584, 463, 591, 488, 563, 553, 549)

# Probit transform
probit = norm.cdf

# Linear model parameters
alpha = Normal('alpha', mu=0.0, tau=0.001, value=0)
beta = Normal('beta', mu=0.0, tau=0.001, value=0)

# Probability of passing
p = Lambda('p', lambda a=alpha, b=beta, x=sat: probit(a+b*x))

# Data likelihood
y = Bernoulli('y', p=p, value=grade, observed=True)

if __name__ == '__main__':
    
    from pymc import MCMC, Matplot
    M = MCMC([alpha, beta, p, y])
    M.isample(20000, 10000, verbose=2)
    Matplot.plot(M)