When I use the function on some simulated data
library(bayesAB)
set.seed(255)
# simulate data
A <- rbinom(250, 1, .25)
B <- rbinom(250, 1, .2)
# apply the function for AB group comparison
AB_test <- bayesTest(A,
B,
priors = c('alpha' = 65, 'beta' = 200),
n_samples = 1e5,
distribution = 'bernoulli')
# obtain the output
summary(AB_test)
I get the following output
# Quantiles of posteriors for A and B:
#
# $Probability
# $Probability$A
# 0% 25% 50% 75% 100%
# 0.1775006 0.2469845 0.2598399 0.2730324 0.3506919
#
# $Probability$B
# 0% 25% 50% 75% 100%
# 0.1510354 0.2146442 0.2268472 0.2394675 0.3182802
#
#
# --------------------------------------------
#
# P(A > B) by (0)%:
#
# $Probability
# [1] 0.89305
#
# --------------------------------------------
#
# Credible Interval on (A - B) / B for interval length(s) (0.9) :
#
# $Probability
# 5% 95%
# -0.04278263 0.37454069
#
# --------------------------------------------
#
# Posterior Expected Loss for choosing A over B:
#
# $Probability
# [1] 0.00587424
I know how to manually obtain the first 3 sections, using the posterior data
quantile(AB_test$posteriors$Probability$A, c(0, 0.25, 0.50, 0.75, 1))
# 0% 25% 50% 75% 100%
# 0.1775006 0.2469845 0.2598399 0.2730324 0.3506919
quantile(AB_test$posteriors$Probability$B, c(0, 0.25, 0.50, 0.75, 1))
# 0% 25% 50% 75% 100%
# 0.1510354 0.2146442 0.2268472 0.2394675 0.3182802
mean(AB_test$posteriors$Probability$A > AB_test$posteriors$Probability$B)
# [1] 0.89305
quantile(AB_test$posteriors$Probability$A / AB_test$posteriors$Probability$B - 1, c(0.05, 0.95))
# 5% 95%
# -0.04278263 0.37454069
But I'm not sure how to calculate the posterior expected loss, shown in the last section of the output. Is that possible?
You are able to figure this out by looking at how the summary is created, i.e. running bayesAB:::summary.bayesTest in the console. Doing this myself I found:
> lapply(
+ AB_test$posteriors,
+ function(x) do.call(bayesAB:::getPostError, unname(x))
+ )
$Probability
[1] 0.00587424
Running bayesAB:::getPostError in the console shows exactly how the calculation works:
> bayesAB:::getPostError
function (A_samples, B_samples, f = function(a, b) (a - b)/b)
{
BoverA <- B_samples > A_samples
loss <- f(B_samples, A_samples)
coalesce(mean(BoverA) * mean(loss[BoverA]))
}
<bytecode: 0x0000017fe9329040>
<environment: namespace:bayesAB>