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 = mc.Normal( "alpha", 0, 0.001, value = 0 ) @mc.deterministic def p( temp = temperature, alpha = alpha, beta = beta): return 1.0/( 1. + np.exp( beta*temperature + alpha) ) observed = mc.Bernoulli( "bernoulli_obs", p, \ value = D, observed=True) model = mc.Model( {"observed":observed, "beta":beta, "alpha":alpha} ) #mysterious code to be explained in Chapter 3 map_ = mc.MAP(model) map_.fit() mcmc = mc.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)