Вычислить байесовский фактор набора данных теста A / B в r

Question

Я пытаюсь вычислить байесовский коэффициент набора тестовых данных A / B, который можно найти здесь . Тем не менее, я получаю NaN, потому что коэффициент бета оценивается как ноль. При расчете вероятностей я предполагаю, что он следует биномиальному распределению. Следовательно, я следую этой формуле:

правдоподобие = выбор (n, k) * бета (k + 1, n-k + 1)

Код можно найти ниже

data <- read.csv(file="ab_data.csv", header=TRUE, sep=",")

control <- data[which(data$group == "control"),]
treatment <- data[which(data$group == "treatment"),]

#compute bayes factor 
n1 = nrow(control)
r1 = sum(control$converted)
n2 = nrow(treatment)
r2 = sum(treatment$converted)

likelihood_control <- choose(n1,r1) * beta(r1+1, n1-r1+1)
likelihood_treatment <- choose(n2,r2) * beta(r2+1, n2-r2+1)
bayes_factor <- likelihood_control/ likelihood_treatment
beta(r1+1, n1+r1+1)
beta(r2+1, n2-r2+1)
bayes_factor

Answer 1

Как вы заметили, проблема в том, что бета-функция возвращает 0, но это не потому, что вероятность на самом деле равна 0, просто вероятность того, что компьютер настолько мал, что компьютер хранит его как 0. Вторая проблема заключается в том, что выбрать возвращает инф. Опять же, это не потому, что значение на самом деле бесконечно, просто R не может внутренне хранить такие большие значения. Решение состоит в том, чтобы использовать логарифмы, которые растут намного медленнее, а затем возводятся в степень в конце. Ниже должно работать (я тестировал функцию logchoose, и, похоже, она работает)

logchoose <- function(n, k){
  num <- sum(log(seq(n - k  + 1, n)))
  denom <- sum(log(1:k))
  return(num - denom)
}

likelihood_control <- logchoose(n1,r1) + lbeta(r1+1, n1-r1+1)
likelihood_treatment <- logchoose(n2,r2) + lbeta(r2+1, n2-r2+1)
bayes_factor <- exp(likelihood_control - likelihood_treatment)
bayes_factor