无法解决模拟问题

Question

我遇到了一个与概率论有关的问题，我试图通过在 R 中模拟它来解决它。但是，我遇到了一个问题，因为 while 循环似乎没有中断。

问题是：需要多少人才能使其中一个人在 12 月的最后一天出生的概率至少为 70%？

这是我的代码：

prob <- 0 
people <- 1 

while (prob <= 0.7) {
  people <- people + 1 #start the iteration with 2 people in the room and increase 1 for every iteration
  birthday <- sample(365, size = people, replace = TRUE) 
  prob <- length(which(birthday == 365)) / people
}
return(prob)

我的猜测是它永远不会达到 70%，因此while 循环永远不会中断，对吗？如果是这样，我是否错误地解释了这个问题？

我不想在 stats.stackexchange.com 上发布此内容，因为我认为这与代码而不是数学本身更相关，但如有必要，我会移动它，谢谢。

Answer 1

在这种情况下，基于概率的分析解决方案比尝试模拟更容易、更准确。我同意 Harshvardhan 的观点，即您的表述解决了错误的问题。

n 池中至少有一个人在特定目标日期生日的概率是1-P{all n miss the target date}。当P{all n miss the target date} < 0.3 时，这个概率至少为 0.7。假设每个人错过目标的概率为P{miss} = 1-1/365（每年 365 天，所有生日的可能性相同）。如果个人生日是独立的，那么P{all n miss the target date} = P{miss}^n。

我不是 R 程序员，但以下 Ruby 应该很容易翻译：

# Use rationals to avoid cumulative float errors.
# Makes it slower but accurate.
P_MISS_TARGET = 1 - 1/365r   
p_all_miss = P_MISS_TARGET
threshold = 3r / 10   # seeking P{all miss target} < 0.3
n = 1
while p_all_miss > threshold
  p_all_miss *= P_MISS_TARGET
  n += 1
end
puts "With #{n} people, the probability all miss is #{p_all_miss.to_f}"

产生：

439 人，所有未命中的概率为 0.29987476838793214

---

附录

我很好奇，因为我的答案与公认的不同，所以我写了一个小模拟。同样，我认为即使它不在 R 中，也很容易理解：

require 'quickstats'  # Stats "gem" available from rubygems.org

def trial
  n = 1
  # Keep adding people to the count until one of them hits the target
  n += 1 while rand(1..365) != 365
  return n
end

def quantile(percentile = 0.7, number_of_trials = 1_000)
  # Create an array containing results from specified number of trials.
  # Defaults to 1000 trials
  counts = Array.new(number_of_trials) { trial }
  # Sort the array and determine the empirical target percentile.
  # Defaults to 70th percentile
  return counts.sort[(percentile * number_of_trials).to_i]
end

# Tally the statistics of 100 quantiles and report results,
# including margin of error, formatted to 3 decimal places.
stats = QuickStats.new
100.times { stats.new_obs(quantile) }
puts "#{"%.3f" % stats.avg}+/-#{"%.3f" % (1.96*stats.std_err)}"

五次运行产生如下输出：

440.120+/-3.336
440.650+/-3.495
435.820+/-3.558
439.500+/-3.738
442.290+/-3.909

这与之前得出的分析结果非常一致，并且似乎与其他响应者的答案存在显着差异。

请注意，在我的机器上，模拟所需的时间大约是解析计算的 40 倍，更复杂，并且引入了不确定性。为了提高精度，您需要更大的样本量，因此需要更长的运行时间。考虑到这些因素，我会重申我的建议，即在这种情况下采用直接解决方案。

Answer 2

确实，您的概率（几乎）永远不会达到 0.7，因为您几乎不会达到恰好 1 人的生日 = 365。当人变大时，会有更多人的生日 = 365，而恰好 1 人的概率会降低。

此外，要计算给定人数的概率，您应该抽取许多样本，然后计算概率。这是实现这一目标的一种方法：

N = 450  # max. number of peoples being tried
probs = array(numeric(), N)  # empty array to store found probabilities

# try for all people numbers in range 1:N
for(people in 1:N){
  # do 200 samples to calculate prop
  samples = 200
  successes = 0
  for(i in 1:samples){
    birthday <- sample(365, size = people, replace = TRUE)
    total_last_day <- sum(birthday == 365)
    if(total_last_day >= 1){
      successes <- successes + 1
    }
  }
  # store found prop in array
  probs[people] = successes/samples
}

# output of those people numbers that achieved a probability of > 0.7
which(probs>0.7)

由于这是一个模拟，结果取决于运行。提高采样率会使结果更稳定。

Answer 3

您正在解决错误的问题。问题是，“需要多少人才能使其中一个人在 12 月的最后一天出生的可能性至少为 70%？”。您现在发现的是“需要多少人才能使 70% 的人在 12 月的最后一天过生日？”。第二个问题的答案接近于零。但是第一个要简单得多。

在您的逻辑中将 prob <- length(which(birthday == 365)) / people 替换为 check = any(birthday == 365)，因为其中至少有一个必须在 12 月 31 日出生。然后，您将能够找到 的人数是否会有至少有一个人出生于 12 月 31 日。

之后，您将不得不多次重新运行模拟以生成经验概率分布（类似于蒙特卡洛）。只有这样你才能检查概率。

模拟代码

people_count = function(i)
{
  set.seed(i)
  for (people in 1:10000)
  {
    birthday = sample(365, size = people, replace = TRUE)
    check = any(birthday == 365)
    if(check == TRUE)
    {
      pf = people
      break
    }
  }
  return(pf)
}

people_count() 函数返回所需的人数，以便其中至少有一个出生于 12 月 31 日。然后我重新运行模拟 10,000 次。

# Number of simulations
nsim = 10000
l = lapply(1:nsim, people_count) %>%
  unlist()

我们来看看所需人数的分布情况。

要找到实际概率，我将使用cumsum()。

> cdf = cumsum(l/nsim)
> which(cdf>0.7)[1]
[1] 292

因此，平均而言，您需要 292 人才能获得超过 70% 的机会。

Answer 4

除了@pjs 答案之外，我想自己提供一个，用 R 编写。我试图通过模拟而不是分析方法来解决这个问题，我分享它以防它对其他人有帮助也有同样的问题。写得不是很好，但想法就在那里：

# create a function which will find if anyone is born on last day
last_day <- function(x){
  birthdays <- sample(365, size = x, replace = TRUE) #randomly get everyone's birthdays
  if(length(which(birthdays == 365)) >= 1) { 
    TRUE #find amount of people born on last day and return true if >1  
  } else {
    FALSE
  }
}

# find out how many people needed to get 70%
people <- 0 #set number of people to zero
prob <- 0 #set prob to zero

while (prob <= 0.7) { #loop does not stop until it hits 70%
  people <- people + 1 #increase the number of people every iteration
  prob <- mean(replicate(10000, last_day(people))) #run last_day 10000 times to find the mean of probability
}
print(no_of_people)

last_day() 只返回TRUE 或FALSE。所以我在每次迭代的循环中运行last_day() 10000 次，以找出在 10000 次中，有多少次有一个或多个人在最后一天出生（这将给出概率）。然后我保持循环运行，直到概率达到 70% 或更多，然后打印人数。

我从运行循环一次得到的答案是440，这与@pjs 提供的答案非常接近。