Saturday 11 January 2020

Birthday Problem and other R snippets

Birthday problem solver code for R. I need a place to store this for future. I am working through second year statistics and will post other interesting things in this page as time goes on.
Here is a link to the full explanation.
In summary
You need 23 people to have at least a 50/50 chance of having a clashing birthday. See the graph below to see how that dynamic works as n increases.
#[1]: 10 people - number of possible birthday combinations
365^10

#Chances of a single combination of interest occurring
1/365^10

#number of combinations within sample space [1] having a unique birthday
(365)*(364)*(363)*(362)*(361)*(360)*(359)*(358)*(357)*(356)

#probability of no one in the room having the same birthday
distinctTuples <- function(elements,sampleSize) 
{
  n = sampleSize
  n_a <- elements
  
  for (x in 1:(n-1)) {
    n_a <- n_a*(elements-x)
  }
  
  return(n_a)
}

distinctTuples(365,10)

#Wrap the distinct tuples up and subtract from one to figure out the birthday problem
birthdayProblem <- function(numberOfPeople) {
  return (1-(distinctTuples(365,numberOfPeople)/365^numberOfPeople))
}

#Print the results
birthdayResults <- lapply(2:50,birthdayProblem)
plot(x=2:50,y=birthdayResults, type="l",main="Birthday Problem", 
     xlab="Number of people", ylab="Chance")

#Binomial distribution function where n=sample size, p=success chance, y=number of 
#successfull outcomes
calc_binomial <- function(p=1,n=1,y=1) {
  q <- 1-p
  s <- factorial(n)/(factorial(y)*factorial(n-y))
  p <- s * (p^y) * (q^(n-y))
  return (p)
}


plot(x=1:20,y=calc_binomial(.5,20,y=1:20), type="l",
     main="20 Coin flips and count of heads per trial", xlab="Number of flips", 
     ylab="Chance")

#quick factorial
easyFact <- function(n,y) {return (factorial(n)/(factorial(y)*factorial(n-y))) }

Output


No comments:

Post a Comment