0

我正在尝试计算可以在此处找到的 A/B 测试数据集的贝叶斯因子。但是，我最终得到一个 NaN，因为 beta 系数的计算结果为零。在计算可能性时，我假设它遵循二项分布。因此，我遵循这个公式：

可能性 = 选择(n,k) * Beta(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
4

1 回答 1

0

正如您所观察到的，问题是 beta 函数返回 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
于 2019-01-06T14:25:27.017 回答