fast-continent-33282
08/16/2023, 2:48 PMhelpful-application-7107
08/16/2023, 3:59 PMfast-continent-33282
08/16/2023, 4:09 PMimport numpy as np
from scipy.stats import norm, beta
from scipy.special import digamma, polygamma, roots_hermitenorm # , roots_sh_jacobi
from scipy import linalg
from scipy.special import betaln, eval_jacobi, roots_sh_jacobi
def update_beta_prior(x_a, n_a, x_b, n_b,alpha_0 = 1, beta_0=1):
alpha_a = alpha_0 + x_a
beta_a = beta_0 + n_a - x_a
alpha_b = alpha_0 + x_b
beta_b = beta_0 + n_b - x_b
return alpha_a, beta_a, alpha_b, beta_b
def calculate_beta_binomial(alpha_a, beta_a, alpha_b, beta_b):
log_beta_mean = lambda a, b: digamma(a) - digamma(a + b)
var_beta_mean = lambda a, b: polygamma(1, a) - polygamma(1, a + b)
d1_beta = norm(loc=beta.mean(alpha_b, beta_b) - beta.mean(alpha_a, beta_a),
scale=np.sqrt(beta.var(alpha_b, beta_b) + beta.var(alpha_a, beta_a)))
d2_beta = norm(loc=log_beta_mean(alpha_b, beta_b) - log_beta_mean(alpha_a, beta_a),
scale=np.sqrt(var_beta_mean(alpha_b, beta_b) + var_beta_mean(alpha_a, beta_a)))
return d1_beta, d2_beta
x_a, n_a = 1287, 2384 # converted & total users in A
x_b, n_b = 1271, 2234 # converted & total users in B
alpha_0, beta_0 = 1, 1 # Beta prior
alpha_a, beta_a, alpha_b, beta_b = update_beta_prior(x_a, n_a,x_b, n_b, alpha_0, beta_0)
d1_beta, d2_beta = calculate_beta_binomial(alpha_a, beta_a, alpha_b, beta_b)
print(f'The probability the conversion in B is higher is {round(d1_beta.sf(0) * 100, 2)}%')
print(f'The 95% credibility interval of (p_b/p_a-1) is {(np.exp(d2_beta.ppf((.025, .975))) - 1)}')
Itamar Farhan
code from medium that i use to double check these numbers.helpful-application-7107
08/16/2023, 4:39 PMfast-continent-33282
08/16/2023, 5:02 PMhelpful-application-7107
08/16/2023, 5:39 PMfast-continent-33282
08/16/2023, 5:41 PMhelpful-application-7107
08/16/2023, 5:43 PMfast-continent-33282
08/16/2023, 5:44 PMhelpful-application-7107
08/16/2023, 5:45 PM