1

I got a question related to probability theory and I tried to solve it by simulating it in R. However, I ran into a problem as the while loop does not seem to break.

The question is asking: How many people are needed such that there is at least a 70% chance that one of them is born on the last day of December?

Here is my code:

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)

My guess is that it could never hit 70%, therefore the while loop never breaks, am I right? If so, did I interpret the question wrongly?

I did not want to post this on stats.stackexchange.com because I thought this is more related to code rather than math itself, but I will move it if necessary, thanks.

4 Answers 4

1

This is a case where an analytical solution based on probability is easier and more accurate than trying to simulate. I agree with Harshvardhan that your formulation is solving the wrong problem.

The probability of having at least one person in a pool of n have their birthday on a particular target date is 1-P{all n miss the target date}. This probability is at least 0.7 when P{all n miss the target date} < 0.3. The probability of each individual missing the target is assumed to be P{miss} = 1-1/365 (365 days per year, all birthdates equally likely). If the individual birthdays are independent, then P{all n miss the target date} = P{miss}^n.

I am not an R programmer, but the following Ruby should translate pretty easily:

# 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}"

which produces:

With 439 people, the probability all miss is 0.29987476838793214


Addendum

I got curious, since my answer differs from the accepted one, so I wrote a small simulation. Again, I think it's straightforward enough to understand even though it's not in 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)}"

Five runs produce outputs such as:

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

which is strongly consistent with the analytical result derived earlier and seems to differ significantly from other responder's answers.

Note that on my machine the simulation takes roughly 40 times longer than the analytical calculation, is more complex, and introduces uncertainty. To increase the precision you would need larger sample sizes, and thus longer run times. Given these considerations, I would reiterate my advice to go for the direct solution in this case.

1
  • I went over the answer I accepted and indeed it wasn't entirely correct. I couldn't entirely figure out what was the problem but it might have been the second part starting from the lapply that might have gone wrong. Anyway, while your answer is not provided in R, I think it's fairly easy to understand with the help of your comments on each line of code. I think its important for people to understand the question first before trying to code and your post did explain the question well, so I have decided to accept your answer. I also provided an answer in R code, but not as well written.
    – Wei
    Dec 6, 2021 at 5:51
1

Indeed, your probability will (almost) never reach 0.7, because you hardly will hit the point where exactly 1 person has got birthday = 365. When people gets larger, there will be more people having a birthday = 365, and the probability for exactly 1 person will decrease.

Furthermore, to calculate a probability for a given number of persons, you should draw many samples and then calculate the probability. Here is a way to achieve that:

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)

As this is a simulation, the result depends on the run. Increasing the sample rate would make the result more stable.

1

You are solving the wrong problem. The question is, "How many people are needed such that there is at least a 70% chance that one of them is born on the last day of December?". What you are finding now is "How many people are needed such that 70% have their birthdays on the last day of December?". The answer to the second question is close to zero. But the first one is much simpler.

Replace prob <- length(which(birthday == 365)) / people with check = any(birthday == 365) in your logic because at least one of them has to be born on Dec 31. Then, you will be able to find if that number of people will have at least one person born on Dec 31.

After that, you will have to rerun the simulation multiple times to generate empirical probability distribution (kind of Monte Carlo). Only then you can check for probability.

Simulation Code

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() function returns the number of people required to have so that at least one of them was born on Dec 31. Then I rerun the simulation 10,000 times.

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

Let's see the distribution of the number of people required.

histogram of empirical distribution

To find actual probability, I'll use cumsum().

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

So, on average, you would need 292 people to have more than a 70% chance.

7
  • 1
    Yup. In case anyone is looking for answer and still unclear, I found this website: bandolier.org.uk/booth/Risk/birthday.html and clears up my confusion about the question. Was about to comment on a solution but thanks!
    – Wei
    Dec 4, 2021 at 17:15
  • 1
    @Wei I'm not so sure about the site you linked, they got the probability of two people sharing the same birthday wrong (they said 1/370, when it's obviously 1/365) which makes me doubt the accuracy of their table. I'd recommend the Wikipedia page on the birthday problem.
    – pjs
    Dec 4, 2021 at 17:30
  • 1
    Since this is a simulation, you should be reporting the margin of error. I came up with a different answer using analytical probability rather than sampling, and it's not clear whether your result is the "same" to within sampling error.
    – pjs
    Dec 4, 2021 at 17:37
  • @pjs if you look at roughly two paragraphs above the table it did specify that the probability of two people sharing a particular birthday is lower than the probabilities in the birthday paradox. I agree that the table isn't accurate and should probably use the wiki page instead, but ultimately this was the page that helped me solve the problem so I thought I'd share.
    – Wei
    Dec 6, 2021 at 5:14
  • 1
    @wei I went straight to the table and saw it was wrong, didn’t bother to read the text at that point. I’m a fan of the birthday problem, it helped me get my very first refereed publication 35 years ago.
    – pjs
    Dec 6, 2021 at 5:58
1

In addition to @pjs answer, I would like to provide one myself, written in R. I attempted to solve this question by simulation rather than an analytical approach, and I am sharing it in case it is helpful for someone else who also has the same problem. Its not that well written but the idea is there:

# 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() only return TRUE or FALSE. So I run last_day() 10000 times in the loop for every iteration to find out, out of 10000 times, how many times does it have one or more people born on the last day (This will give the probability). I then keep the loop running until the probability is 70% or more, then print the number of people.

The answer I get from running the loop once is 440 which is quite close to the answer provided by @pjs.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.