Bayesian credible interval for the difference between two proportions

bayesian_diff_prop.py

import numpy as np
import scipy.stats as stats

# Number of trials and successes for each group
n_1 = 100  
x_1 = 6   
n_2 = 125 
x_2 = 5    

# Parameters for the prior distributions
# Using non-informative priors with parameters alpha=1 and beta=1
alpha_prior = 1
beta_prior = 1

# Update the parameters for the posterior distributions
alpha_post1 = alpha_prior + x_1
beta_post1 = beta_prior + n_1 - x_1
alpha_post2 = alpha_prior + x_2
beta_post2 = beta_prior + n_2 - x_2

# Define the confidence level
confidence = 0.95

# Generate samples from the posterior distributions
samples1 = stats.beta.rvs(alpha_post1, beta_post1, size=10000)
samples2 = stats.beta.rvs(alpha_post2, beta_post2, size=10000)

# Calculate the difference between the samples
diff_samples = samples1 - samples2

# Calculate the credibile interval for the difference
ci_diff = np.percentile(diff_samples, [(1-confidence)/2*100, (1+confidence)/2*100])

print(f"The {confidence * 100}% credibile interval for the difference between the proportions is {ci_diff}")