解决模拟问题

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 上 post 因为我认为这更与代码相关而不是数学本身，但如果需要我会移动它，谢谢。

Answer 1

事实上，你的概率（几乎）永远不会达到 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 2

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

把你逻辑中的prob <- length(which(birthday == 365)) / people换成check = any(birthday == 365)，因为其中至少有一个人必须在12月31日出生。然后，你就能找到那个人数将至少有一个人在 12 月 31 日出生。

之后，您将不得不多次重新运行模拟以生成经验概率分布（有点像 Monte Carlo）。只有这样你才能检查概率。

模拟代码

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() 函数 returns 至少有一个人在 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 3

在这种情况下，基于概率的分析解决方案比尝试模拟更容易、更准确。我同意 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}"

产生：

With 439 people, the probability all miss is 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 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() 仅 return TRUE 或 FALSE。所以我运行 last_day() 每次迭代循环 10000 次以找出在 10000 次中有多少次有一个或多个人在最后一天出生（这将给出概率).然后我保持循环运行ning 直到概率为 70% 或更多，然后打印人数。

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

解决模拟问题

Having trouble solving simulation

simulation

r

probability

模拟代码

附录