lin-reg-ppc.py

model = Model()
with model:
  beta = Uniform('beta', lower=-10, upper=10, shape=(1, n_predictors))
  
  p = tinvlogit(sum(beta * predictors, 1))
  Y_obs = Bernoulli('Y_obs', p=p, observed=Y)
  
  start = find_MAP()
  step = NUTS(scaling=start)
  trace = sample(mcmc_training_sample_size, step, start)

ppc_sample_size = 500   
beta_rnd = trace[(-1 * ppc_sample_size):]['beta']

# get point to test and calculate posterior distribution 
pred_test = get_point_to_test() 
ps_np = tinvlogit(sum(beta_rnd * pred_test), 1) 
posterior_mean = ps_np.mean()
posterior_sigma = ps_np.std() 
posterior_interval = hpd(ps_np)