figsize(12.5, 6)
subplot(211)
plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75,
            color="k", alpha=0.5)
plt.yticks([0, 1])
plt.ylabel("Damage Incident?")
plt.title("(Real) Defects of the Space Shuttle O-Rings vs temperature")

subplot(212)
n = challenger_data.shape[0]
plt.scatter(challenger_data[:, 0], stats.bernoulli.rvs(0.6, size=n),
            s=75, color="k", alpha=0.5)
plt.yticks([0, 1])
plt.ylabel("Damage Incident?")
plt.xlabel("Outside temperature (Farhenhit)")
plt.title("(Artificial) Defects of the Space Shuttle O-Rings vs \
temperature")

alpha_hat = alpha_samples[0, 0]  # select the first alpha sample
beta_hat = beta_samples[0, 0]  # select the first beta sample

p_hat = logistic(temperature, beta_hat, alpha_hat)
print "estimates of probability at observed temperature, defects: "
print np.array(zip(p_hat, temperature, D))[:3, :]
print "..."

print p_hat
_product = p_hat ** (D) * (1 - p_hat) ** (1 - D)
prob_given_model_1 = _product.prod()
print "probability of observations, model 1: %.10f" % prob_given_model_1

beta = 0
alpha = pm.Normal("alpha", 0, 0.001, value=0)


@pm.deterministic
def p(temp=temperature, alpha=alpha, beta=beta):
    return 1.0 / (1. + np.exp(beta * temperature + alpha))


observed = pm.Bernoulli("bernoulli_obs", p,
                        value=D, observed=True)

model = pm.Model({"observed": observed, "beta": beta, "alpha": alpha})

# mysterious code to be explained in Chapter 3
map_ = pm.MAP(model)
map_.fit()
mcmc = pm.MCMC(model)
mcmc.sample(260000, 220000, 3)
######

alpha_samples_model_2 = mcmc.trace("alpha")[:, None]
alpha_hat = alpha_samples_model_2[0]  # use the first 'model'
beta_hat = 0

p_hat = logistic(temperature, beta_hat, alpha_hat)
print "estimates of probability at observed temperature, defects: "
print np.array(zip(p_hat, temperature, D))[:3, :]
print
print "Notice the probability is constant for all temperatures?"

# compute the probability of observations given the model_2

_product = p_hat ** (D) * (1 - p_hat) ** (1 - D)
prob_given_model_2 = _product.prod()
print "probability of observations, model 2: %.10f" % prob_given_model_2

print "Bayes factor = %.3f" % (prob_given_model_1 / prob_given_model_2)

p = logistic(temperature[None, :], beta_samples, alpha_samples)
_product = p ** (D) * (1 - p) ** (1 - D)
prob_given_model_1 = _product.prod(axis=1).mean()
print "expected prob. of obs., given model 1: %.10f" % prob_given_model_1


p_model_2 = logistic(temperature[:, None],
                     np.zeros_like(temperature), alpha_samples_model_2)

_product = p_model_2 ** (D) * (1 - p_model_2) ** (1 - D)
prob_given_model_2 = _product.prod(axis=1).mean()
print "expected prob. of obs., given model 2: %.10f" % prob_given_model_2
print
print "Bayes factor: %.3f" % (prob_given_model_1 / prob_given_model_2)