Green Eyes

September 29, 2009

We’ll look first at the standard closed-form solution. The (choose n r) function calculates the number of ways to choose r items from a group of n items:

(define (choose n r)
  (let loop ((n n) (r r) (c 1))
    (if (zero? r) c
      (loop (- n 1) (- r 1) (* c n (/ r))))))

There are 27C3 ways to choose three people from a group of twenty-seven. (That 27C3 is the notation used in my daughter’s high-school textbook. When I was in high school, the notation was \left( \begin{array}{c} 27 \\ 3 \end{array} \right)_C.) Of those, 24C2 will have one person with green eyes (that is, choose two of the twenty-four people who don’t have green eyes). Since any of the three can have green eyes, the result must be multiplied by three:

> (* (/ (choose 24 2) (choose 27 3)) 3)
92/325

There are 27C3 = 2925 possible sets of three people, of which 3 × 24C2 = 828 have one person with green eyes. The fraction 828/2925 reduces to 92/325, or 28.3%.

We could instead enumerate all 2925 collections of three people and count those that have one person with green eyes. Here we define a vector of length twenty-seven with the appropriate mix of blue, brown and green:

(define eyes (make-vector 27 'blue))
(for (i 0 3) (vector-set! eyes i 'green))
(for (i 3 16) (vector-set! eyes i 'brown))

We also need a function that checks the color of a person’s eyes:

(define (green x)
  (if (eq? (vector-ref eyes x) 'green) 1 0))

Now we form the 2925 possibilities and count how many have one person with green eyes:

> (length
    (list-of (list a b c)
      (a range 0 27)
      (b range (+ a 1) 27)
      (c range (+ b 1) 27)
      (= (+ (green a) (green b) (green c)) 1)))
828

A third approach is to simulate a large number of selections from appropriate groups of people and count those that have one person with green eyes. We need to select three persons from a population of twenty-seven, or more generally m items from a set of n, without replacement. We can make the selection by running through the items x in order, selecting each with probability s/r, where s is the remaining number of items to be selected and r is the remaining number of items in the selection pool (this is Knuth’s algorithm 3.4.2S):

(define (sample n m)
  (let loop ((s m) (r n) (x 0) (xs '()))
    (cond ((= x n) xs)
          ((< (rand) (/ s r))
            (loop (- s 1) (- r 1) (+ x 1) (cons x xs)))
          (else (loop s (- r 1) (+ x 1) xs)))))

Then we can simulate n tests:

(define (sim n)
  (define (s?)
    (= (sum (map green (sample 27 3))) 1))
  (let loop ((n n) (g 0))
    (if (zero? n) g
      (loop (- n 1) (+ g (if (s?) 1 0))))))

Here’s a sample:

> (sim 2925)
827

That our sample run randomly gave an answer so close to the true answer is an example of sheer dumb luck. Math problems like this let us show off the Standard Prelude; we used sum, for, list-of, and rand. You can run the program at http://codepad.org/gSUgP0pn.

About these ads

Pages: 1 2

5 Responses to “Green Eyes”

  1. P. Riva said

    I’d say we have

    P = 3/27 * 24/27 * 24/27 * 3 =~ 26.34%

  2. P. Riva said

    OPS!!! I was a bit in a hurry:

    3/27 * 24/26 * 23/25 * 3 =~ 28.31%

  3. Peter said

    It may be easier to conceptualize as choosing one child w/ green eyes and two children with eyes of any other color. Your numerator is P(choosing one green of the 3 greens)*P(choosing 2 not greens of the 24 not greens). The denominator is the total number of combinations of choosing three children out of the 27.

    Number of ways of:
    -Choosing 1 green out of the 3 greens: 3C1 = 3
    -Choosing 2 non-greens out of the 24 non-greens: 24C2 = 276
    -Choosing 3 kids out of 27 kids: 27C3 = 2925

    So, the probability of choosing only one green eyed kid out of the three selected is:
    (3C1*24C2)/(27C3) = (3*276)/2925 = 0.283076923

  4. Graham said
    from operator import mul
    
    def falling_factorial(n, k):
        """(n)_k = n! / (n - k)! = n * (n - 1) * ... * (n - k + 1)"""
        return reduce(mul, xrange(n - k + 1, n + 1), 1)
    
    
    def binom(n, k):
        """n! / (k! * (n-k)!) = (n)_k / k!; k! = (k)_k"""
        return falling_factorial(n, k) / falling_factorial(k, k)
    
    
    if __name__ == "__main__":
        print 100 * binom(24, 2) * binom(3, 1) / float(binom(27, 3))
    

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 621 other followers

%d bloggers like this: