def plot_ci(p, post, num_samples, lower_q, upper_q):
 # Compute a large sample by resampling with replacement
 samples = np.random.choice(p, size=num_samples, replace=True, p=post)
 ci = scipy.percentile(samples, [lower_q*100, upper_q*100]) # compute the quantiles

 interval = upper_q - lower_q
 plt.title('Posterior density with %.3f credible interval' % interval)
 plt.plot(p, post, color='blue')
 plt.xlabel('Parameter value')
 plt.ylabel('Density')
 plt.axvline(x=ci[0], color='red')
 plt.axvline(x=ci[1], color='red')
 print('The %.3f credible interval is %.3f to %.3f' 
 % (interval, lower_q, upper_q))

plot_ci(p, post, num_samples, lower_q, upper_q)