我有一個與概率論相關的問題，我試圖通過在 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 上發布這個，因為我認為這與代碼而不是數學本身更相關，但我會在必要時移動它，謝謝。
uj5u.com熱心網友回復：
您正在解決錯誤的問題。問題是，“需要多少人才能至少有 70% 的機會在 12 月的最后一天出生？”。您現在發現的是“需要多少人才能使 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% 的機會。
uj5u.com熱心網友回復：
實際上，您的概率（幾乎）永遠不會達到 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)
由于這是一個模擬，結果取決于運行。增加采樣率會使結果更穩定。
uj5u.com熱心網友回復：
在這種情況下，基于概率的決議解比嘗試模擬更容易、更準確。我同意 Harshvardhan 的觀點，即您的表述解決了錯誤的問題。
在 n 個池中至少有一個人在特定目標日期過生日的概率是1-P{all n miss the target date}。當 時，這個概率至少為 0.7 P{all n miss the target date} < 0.3。假設每個人未達到目標的概率為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 倍，更復雜，并引入了不確定性。為了提高精度，您需要更大的樣本量，因此需要更長的運行時間。鑒于這些考慮，我會重申我的建議，在這種情況下采用直接解決方案。
轉載請註明出處，本文鏈接：https://www.uj5u.com/net/373013.html
上一篇：使用具有功能真/假值的dplyr的if_else函式
下一篇：如何將正交多項式分配給函式？