Birthday Paradox

October 12, 2012

We first write a function to determine if two or more people in a group of n people have the same birthday; it is sufficient to keep a list of random numbers as they arrive, compare each new number to the list, and report #t if there is a match and #f if you reach the desired size without a match:

(define (same? n)
  (let loop ((n n) (rs (list)))
    (if (zero? n) #f
      (let ((r (randint 365)))
        (if (member r rs) #t
          (loop (- n 1) (cons r rs)))))))

To verify the percentages, we simulate a number of trials:

(define (sim size trials)
  (let loop ((trials trials) (same 0))
    (if (zero? trials) same
      (loop (- trials 1)
        (+ same (if (same? size) 1 0))))))

And here are our results:

> (sim 23 1000)
500
> (sim 57 1000)
989

I was shocked when the first simulation reported exactly 500 matches; I expected something close but not exact. You can run the simulation at http://programmingpraxis.codepad.org/mK7SkBOG, where you will also see the random number generator from the Standard Prelude.

About these ads

Pages: 1 2

16 Responses to “Birthday Paradox”

  1. [...] today’s Programming Praxis exercise, our goal is to simulate the well-known birthday paradox. Let’s [...]

  2. My Haskell solution (see http://bonsaicode.wordpress.com/2012/10/12/programming-praxis-birthday-paradox/ for a version with comments):

    import Control.Monad
    import Data.List
    import System.Random
    
    paradox :: Int -> IO Bool
    paradox n = fmap (\g -> or [elem h t | (h:t) <- tails . take n $
                randomRs (1, 365 :: Int) g]) newStdGen
    
    trial :: Int -> IO Float
    trial = fmap ((/ 100) . genericLength . filter id) . replicateM 10000 . paradox
    
    main :: IO ()
    main = do print =<< trial 23
              print =<< trial 57
    
  3. Steve Massey said

    Should one not allow for leap years?

  4. Paul said

    A Python version

    import random as RA
    import itertools as IT

    def analytic(N):
    prob = 1
    for i in range(0, N):
    prob *= (1.0 – i / 365.0)
    return 1 – prob

    def gen_cases(N):
    while 1:
    yield [RA.randint(1, 365) for i in range(N)]

    def monte_carlo(N, mcsamples = 1000):
    “””calculate chance that 2 or more people have same birthday”””
    cnt = sum(len(c) != len(set(c)) for c in IT.islice(gen_cases(N), mcsamples))
    return cnt / float(mcsamples)
    #[/sourcecode]

  5. Paul said

    Repost. A python version.

    import random as RA
    import itertools as IT
    
    def analytic(N):
        prob = 1
        for i in range(0, N):
            prob *= (1.0 - i / 365.0)
        return 1 - prob
    
    def gen_cases(N):
        while 1:
            yield [RA.randint(1, 365) for i in range(N)]
    
    def monte_carlo(N, mcsamples = 1000):
        """calculate chance that 2 or more people have same birthday"""
        cnt = sum(len(c) != len(set(c)) for c in IT.islice(gen_cases(N), mcsamples))
        return cnt / float(mcsamples)
    
  6. tompko said

    A Python version which takes into account leap years, 0-364 are the non-leap days of the year and 365 is Feb 29th:

    import random
    
    def gen_day():
        return random.randrange(366) if random.randrange(4) == 0 else random.randrange(365)
    
    def trial(people):
        birthdays = [gen_day() for p in range(people)]
        return len(birthdays) != len(set(birthdays))
    
    def simulate(people, trials):
        return len([x for x in range(trials) if trial(people)])
    
  7. JP said

    Hah, beat you to one for once. :)

    Here’s a post that I wrote a week and a half ago on my own birthday that basically works out the math behind the problem and then lets you run a simulation build into the web page (Javascript w/ JQuery). You can choose any number of birthdays to simulate and it will generate that many for you over the next year, keeping statistics for each number, along with the expected value. I think it’s pretty neat at least. :)

    Check it out here: The Birthday Paradox

  8. David said

    Clojure solution.

    (defn birthday? [n]
       (let [sample  (map rand-int (repeat n 365))]
          (not= (count (distinct sample))  (count sample))))
    
    (defn paradox [n, trials]
       (float (/ (count (filter birthday? (repeat trials n))) trials)))
    
    user=> (paradox 23 100000)
    0.50874
    user=> (paradox 33 100000)
    0.7745
    user=> (paradox 57 100000)
    0.99029
    
  9. Christian Siegert said

    Solution written in Go (source code https://gist.github.com/3885112):

    package main

    import (
    "fmt"
    "math/rand"
    )

    func main() {
    fmt.Printf("Chance of same birthday in group with 23 people: %.0f%%\n", getChanceOfSameBirthdayInGroup(23)*100)
    fmt.Printf("Chance of same birthday in group with 57 people: %.0f%%\n", getChanceOfSameBirthdayInGroup(57)*100)
    }

    func getChanceOfSameBirthdayInGroup(groupSize uint) float32 {
    n := 10000
    numberOfSameBirthdays := 0

    for i := 0; i < n; i++ {
    if isSameBirthdayInGroup(groupSize) {
    numberOfSameBirthdays++
    }
    }

    return float32(numberOfSameBirthdays) / float32(n)
    }

    func isSameBirthdayInGroup(groupSize uint) bool {
    birthdays := make(map[int]uint8)

    for i := uint(0); i < groupSize; i++ {
    birthday := rand.Intn(365)
    birthdays[birthday]++

    if birthdays[birthday] == 2 {
    return true
    }
    }

    return false
    }

  10. treeowl said

    A year whose number is divisible by 100 is not a leap year, unless it is also divisible by 400. So the year 2000 was a leap year, but the year 2100 won’t be. A program trying to calculate birthday statistics for people currently alive can ignore this issue, but if it’s supposed to work across history it needs to take it into account. Of course, if one goes back before the Gregorian calendar was adopted, matters get extremely complicated, and figuring out how to deal with the eventual breakdown (at some not-yet-predictable time) of the Gregorian calendar is even trickier.

  11. A Java Solution:

    import java.util.HashSet;
    import java.util.Random;
    import java.util.Set;
    
    public class BirthdayParadox {
    
        private static final Random randGen = new Random();
    
        static float simulate(int nPeople, int nTrials) {
            int n = 0;
            for(int i = 0; i < nTrials; i++) {
                if(trial(nPeople)) n++;
            }
            return (float)n/nTrials;
        }
    
        static boolean trial(int nPeople) {
            Set<Integer> birthdays = new HashSet<Integer>();
            for (int i = 0; i < nPeople; i++) {
                int day = randGen.nextInt(365);
                if (birthdays.contains(day)) return true;
                birthdays.add(day);
            }
            return false;
        }
    
        public static void main(String[] args) {
            System.out.printf("%.4f\n", simulate(23, 1000)); # => 0.5220
            System.out.printf("%.4f\n", simulate(57, 1000)); # => 0.9900
        }
    }
    
  12. drew said

    Leap years? Not sure what leap years have to do with the original problem stated in the main blog homepage. Here is my PHP solution for the first part.

    for($i = 1; $i < 1000; $i++)
    {

    $arr = array();

    for($q = 0; $q 1) { $total++; }

    }

    echo $total / 1000;

  13. drew said

    Not sure what’s wrong with your blog but it just can’t off half my code. Forget this. Moving on.

  14. PBC said
    (defun sim-bday-paradox (size trials)
      (loop repeat trials
         counting (loop repeat size
    		 for r = (random 365) then (random 365)
    		 when (member r bday-lst) do (return t)
    		 collect r into bday-lst)))
    (sim-bday-paradox 23 10000)
    (sim-bday-paradox 57 10000)
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 574 other followers

%d bloggers like this: