有没有一种方法可以用R从一个大样本蒙特卡罗模拟中确定一个精确的分数？

Question

在过去的一个周末，我参加了一年一度的麻省理工学院迷思狩猎。这个谜题有12个我们必须解决的概率问题。虽然我们没有及时完成它，但我决定自己尝试解决它(所以没有破坏者)。

问题8目前我被绊倒了。它说:爱丽丝从小费中得到16个硬币，鲍勃得到3个硬币。他们决定玩下面的游戏。每个玩家都会选择自己的硬币。他们掷硬币，直到一个是正面，另一个是尾巴。有头的人把两枚硬币都留着。这是重复的，直到一个玩家得到所有的硬币。爱丽丝把小费丢给鲍勃的可能性有多大？

我不擅长统计学，我相信有一种方法可以用公式来计算这一点(这里的任何指导都会受到赞赏)，但我意识到我可以通过模拟R中的结果来计算出一个答案。

下面是我为模拟程序运行的代码(警告，对R来说是新的，所以很抱歉，jankiness和糟糕的格式化)：

coin_battle <- function(alice_points,bob_points){
Alice <- alice_points
Bob <- bob_points
    while(Alice !=0 & Bob !=0){
        result <- sample(1:2,1)
        if(result == 1){
            Alice <- Alice + 1
            Bob <- Bob - 1
        }
        else{
            Alice <- Alice - 1
            Bob <- Bob + 1
        }
    }
if(Alice==0){
    return("Bob")
}
if(Bob==0){
    return("Alice")
}
}
simulate <- function(sample_size){
n <- sample_size
res <- rep(0,n)
    for(i in 1:n){
        res[i] <- coin_battle(16,3)
    }
    return(res)
}
vector1 <- simulate(10000)
table(vector1)

我刚刚运行了这个样本大小为10,000,000的代码，结果是，Alice赢了8,420,661次，Bob赢了1,579,339次，所以Bob赢了15.79339%。我让一个朋友确认这个十进制答案是正确的。然而，为了这个谜题的目的，答案必须是一个分数的形式。I怀疑简化形式的分子和分母都是26或更小，因为这个谜题可能会让我把分数转换成字母。我可能可以用1到26之间的数字来猜测和检查潜在的选项，但是我想知道是否有一个更通用的解决方案。

谢谢!

Answer 1

我想不出一种直接的方法来估计准确的分数。但是你有一个很大的模拟，所以你可以在真实概率的周围设置一个置信区间。在这里，你可以假设分母小于某个任意值，检查所有符合这个标准的分数，看看它们是否落在置信区间内。我写了一些代码来做这件事，因为这个问题吸引了我。请注意-这是非常糟糕的代码，因为我有点不舒服，思维也不太清楚。

ps <- prop.test(1579339,10000000)

# arbitrarily assume the denominator is <= 100
denominators <- 1:100
# numerators
num <- lapply(denominators, function(x){sequence(x)})
# the decimal values
vals <- lapply(denominators, function(x){sequence(x)/x})

# helper function to pluck the values that fall between the interval
between <- function(x, a, b){
  a <= x & x <= b
}

# find the values that fall in the interval
btw <- lapply(vals, between, a = ps$conf.int[1], b = ps$conf.int[2])
# not too proud of this.  it's ugly and there is a better way
output <- data.frame(percs = unlist(vals),
                     num = unlist(num), 
                     den = rep(denominators, denominators), 
                     between = unlist(btw))
# only grab the output that meets our criteria
possible <- output[output$between,]
possible

这给了我们

> possible
         percs num den between
174  0.1578947   3  19    TRUE
709  0.1578947   6  38    TRUE
1605 0.1578947   9  57    TRUE
2862 0.1578947  12  76    TRUE
4480 0.1578947  15  95    TRUE

所有这些都减少到了3/19。