Calculate posterior expected loss manually from bayesTest function of bayesAB package-CodePudding

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?

CodePudding user response：

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>