Things to explore

This looks very interesting:



As does this:

Towers of Hanoi

This famous mathematical game is described in wikipedia. Briefly, we have three pegs and N disks. We wish to move the whole stack of disks to a different peg, let's say # 3. The rules are:

• Only one disk may be moved at a time.
• Each move consists of taking the upper disk from one of the rods and sliding it onto another rod, on top of the other disks that may already be present on that rod.
• No disk may be placed on top of a smaller disk.

For example, here is an intermediate stage in the game with 4 disks. Can you see the three moves that brought us to this point?



A couple of interesting things: one is the recursive nature of Towers of Hanoi. If you already know how to solve the game with N disks, then how do you solve the game with N+1?

Easy, move the first N disks to peg # 2, then move disk N+1 to peg # 3, then move all the other N disks on top. This stereotyped pattern leads to the following visual aid. (I've forgotten where I saw it). It looks like a binary ruler.



This version of the ruler describes the series of moves (though not the target pegs) for the N=4 game. To extend it, add a new bar of the correct height for N = 5, then duplicate all the bars we already have.

The number of moves grows rapidly:

N moves 1 1 2 3 3 7 4 15 N 2N-1

According to wikipedia:

The puzzle was invented by the French mathematician Édouard Lucas in 1883. There is a legend about a Vietnamese temple which contains a large room with three time-worn posts in it surrounded by 64 golden disks. The monks of Hanoi, acting out the command of an ancient prophecy, have been moving these disks, in accordance with the rules of the puzzle, since that time. The puzzle is therefore also known as the Tower of Brahma puzzle. According to the legend, when the last move of the puzzle is completed, the world will end.

We can calculate, at one move per second, this game will take roughly 292 billion years, about 20 times the current age of the universe.

>>> x = 3600*24*365 >>> x 31536000 >>> 2**63/x 292471208677L

tic-tac-toe: game on!

This is my last post about this topic (see sideboard for more). What I have is a fairly long listing (250 lines) which uses the utilities from here, and analyzes the outcomes of playing tic-tac-toe. You can run it on any computer with Python installed (e.g. any Mac).

It has four major sections. (1) We build up game positions as shown previously, through move # 3. We make sure to check that a new position is not a permutation of one already known. Also we move from the center out, so that if there are two identical positions, we keep the one in which the center was X's first move. (2) We then do move # 4, continuing as before, but also checking to be sure that if X threatens to win, we block him. (3) For moves # 5-9 we check for several things in turn:

• is the position already won
• can the player win
• can the opponent win
• can the player make a "double threat"---and win on next turn.
• can the opponent make a double threat---then we block.


The last part organizes the output. All the results from both the intermediate stages and either won or filled (drawn) boards are stored in a dictionary labeled mD, keyed by the move list for that position (1-based indexing) For example, 27853149 means X took square 2, then O took square 7, and so on. We also store the logic for the game in another dictionary labeled rD, so for each position, we can say what the reason was for us to make the last move. We filter for boards which terminate a chain by looking for those with no "daughters"---that is, there are no keys in the dictionary which begin with this key.

The positions could be analyzed in a variety of ways. Here are two:
(1) If you've played you know that it is difficult to win if you do not go first. But it is not impossible! There are six winning games for O. Here is my favorite:

27314985 . X X O X X O X X O X X . . . . . . X . . X . . O . . O . . O . . O . O 273 2731 27314 273149 block X block center O X X O X X X . . X O . O X O O X O 2731498 27314985 block winner:O

We see that X took the side position, then O the corner, then X the opposing corner (273). It looks pretty good for X! O has to block on the next move. But now O threatens and X must block, and finally O has an unanswerable double threat.

The other question I looked at was: does X always win if he takes the center first? The answer is surprising. X wins only 4/10. The others are draws. Here is one of the wins:

5284793 . O . . O . . O . . O . . X . O X . O X . O X . . X . . X . X X . X X O 528 5284 52847 528479 DT for X block . O X O X . X X O 5284793 winner:X

Here is the code:

import sys from TTTUtils import * mD = dict() # master (as dict) # for moves 1-3 we keep all positions # as long as they are not permutations bL = [None] * 9 # "board" list # move #1 for i in [0,1,4]: L = bL[:] L[i] = 'X' mD[str(i+1)] = L def test(n): kL = [k for k in mD.keys() if len(k) == n] for k in sorted(kL,reverse=True): if len(k) == n: print k print boardRep(mD[k]) print '---' sys.exit() #test(1) # we do reverse sort to keep moves # where center(5) was first kL1 = [k for k in mD if len(k) == 1] kL1.sort(reverse=True) # move #2 added = list() for k in kL1: bL = mD[k] for i in openPos(bL): bL2 = bL[:] bL2[i] = 'O' good = True # now test them if bL2 in added: continue for p in allPermutations(bL2): if p in added: good = False break if good: added.append(bL2) mD[k + str(i+1)] = bL2 #test(2) # move #3 added = list() kL2 = [k for k in mD.keys() if len(k) == 2] kL2.sort(reverse=True) for k in kL2: bL = mD[k] for i in openPos(bL): bL2 = bL[:] bL2[i] = 'X' good = True # now test them if bL2 in added: continue for p in allPermutations(bL2): if p in added: good = False break if good: added.append(bL2) mD[k + str(i+1)] = bL2 #test(3) #========================================= # move #4, check for X threatens win # we are still checking permutations # from now on, we keep some info rD = dict() # what reason for move? added = list() kL3 = [k for k in mD if len(k) == 3] kL3.sort(reverse=True) for k in kL3: bL = mD[k] # a list of i,line for all winners iL = winners(bL,who='X') if iL: i = iL[0][0] bL2 = bL[:] bL2[i] = 'O' k2 = k + str(i+1) mD[k2] = bL2 rD[k2] = 'block X' added.append(bL2) continue for i in openPos(bL): bL2 = bL[:] bL2[i] = 'O' good = True # now test them if bL2 in added: continue for p in allPermutations(bL2): if p in added: good = False break if good: added.append(bL2) mD[k + str(i+1)] = bL2 #test(4) def testagain(): kL = [k for k in mD if len(k) == 4] print '4 move stage: N =', len(kL) sys.exit() #testagain() #========================================= # moves #5-9 def round(N): if N % 2: p1 = 'X'; p2 = 'O' else: p1 = 'O'; p2 = 'X' kL = [k for k in mD if len(k) == N-1] kL.sort(reverse=True) for k in kL: bL = mD[k] # check already have 3 in a row if wonPosition(bL): continue C = False # winners for player, then opponent for who in [p1,p2]: if C: continue # returns list of (i,line) iL = winners(bL,who) if iL: bL2 = bL[:] i = iL[0][0] bL2[i] = p1 k2 = k + str(i+1) mD[k2] = bL2 # save the reason if who == p1: rD[k2] = 'winner:' + p1 else: rD[k2] = 'block' C = True if C: continue # double threat for player,opponent for who in [p1,p2]: if C: continue # DT tries each move, returns i # if a move has two ways to win rL = doubleThreat(bL,who) if rL: i = rL[0] bL2 = bL[:] bL2[i] = p1 k2 = k + str(i+1) mD[k2] = bL2 # save the reason if who == p1: rD[k2] = 'DT for X' else: rD[k2] = 'DT for O' C = True # else, if the center is open, take it if not bL[4]: bL2 = bL[:] bL2[i] = p1 k2 = k + str(i+1) mD[k2] = bL2 rD[k2] = 'center' C = True if C: continue rL = promisingLines(bL,'X') for i,line in rL: if C: continue bL2 = bL[:] bL2[i] = p1 k2 = k + str(i+1) mD[k2] = bL2 rD[k2] = 'possible' C = True if C: continue if N == 9: i = bL.index(None) bL2 = bL[:] bL2[i] = p1 k2 = k + str(i+1) mD[k2] = bL2 rD[k2] = 'last' continue print boardRep(bL) print 'found nothing' sys.exit() for i in range(5,10): round(i) #----------------------------------------- def daughters(k,mD): rL = list() n = len(k) for k2 in mD: if len(k2) == n + 1: if k2[:n] == k: rL.append(k2) return rL #for k in mD: if len(k) == 9: print k def output(kL): pL = list() pL.append(multipleReps(mD,kL[:4]).strip()) s = '' for k in kL[:4]: if k in rD: s += rD[k].ljust(11) else: s += ' ' * 11 pL.append(s) pL.append('') pL.append(multipleReps(mD,kL[4:]).strip()) s = '' for k in kL[4:]: if k in rD: s += rD[k].ljust(11) else: s += ' ' * 11 pL.append(s) pL.append('-' * 30) return '\n'.join(pL) for k in sorted(mD.keys(),reverse=True): continue if len(k) < 4: continue if not daughters(k,mD): # look at wins for O if not 'winner' in rD[k]: continue if not 'O' == rD[k][-1]: continue print k kL = list() for i in range(3,len(k)+1): kL.append(k[:i]) print output(kL) for k in sorted(mD.keys(),reverse=True): if len(k) < 4: continue if not k[0] == '5': continue if not daughters(k,mD): print k #continue kL = list() for i in range(3,len(k)+1): kL.append(k[:i]) print output(kL)

tic-tac-toe: loop detail

In the program I'm writing, I have code structured like this:

for i in range(2): print 'i', i for p in 'xy': print 'p', p if True: continue

What I want is to continue back to the next i (outer loop), rather than the next p (inner loop). This code doesn't do that:

i 0 p x p y i 1 p x p y

It's not elegant, but what I did was add a flag labeled C. We check the value of C before starting the inner loop.

for i in range(2): print 'i', i C = False for p in 'xy': if C: continue print 'p', p if True: C = True

That does what I want:

i 0 p x i 1 p x

tic-tac-toe: building our repertoire

Now let's go through the first three rounds of building a library of all possible boards. A previous post did this by hand. Unfortunately, it shows. There are errors in the previous boards, situations where I did not recognize that two boards were the same under rotation or reflection.

The code is not very pretty at this stage, but it is typical of projects that I've done. While I'm still getting things to work, there are lots of extra print statements, and there are places where a similar operation is repeated, just because that seemed to be faster to write. We'll sort it out later. The main thing is to look carefully at the output to check that when we reject "duplicates" and "permutations" the decision is correct.

Here is a small part of the output at round 3:

123 X O X . . . . . . is accepted 124 X O . X . . . . . is accepted /snip/ 154 X . . X O . . . . is a perm of 152 X X . . O . . . .

A partial list of the mistakes from before (and the dups) is:

163 (127) 169 (129) 193 (137) 213 (132) 217 (136) 219 (196) 241 (243) ...

By my count now, we have 53 possible boards at stage 3. Next, going on to stage 4. But first we will have to analyze the boards for forced moves.

import sys from TTTUtils import * mD = dict() # master (as dict) # round 1 bL = [None] * 9 # "board" list for i in [0,1,4]: L = bL[:] L[i] = 'X' mD[str(i+1)] = L def test(n): for k in sorted(mD.keys()): if len(k) == n: print k print boardRep(mD[k]) print '---' sys.exit() #test(1) # round 2 added = list() kL = [k for k in mD.keys() if len(k) == 1] kL.sort() #kL.reverse() for k in kL: bL = mD[k] for i in openPos(bL): bL2 = bL[:] bL2[i] = 'O' good = True # now test them if bL2 in added: continue for p in allPermutations(bL2): if p in added: good = False break if good: added.append(bL2) mD[k + str(i+1)] = bL2 #test(2) def find_kforp(p): for k in mD: if mD[k] == p: return k return 'not found' # round 3 added = list() kL = [k for k in mD.keys() if len(k) == 2] kL.sort() #kL.reverse() for k in kL: bL = mD[k] for i in openPos(bL): bL2 = bL[:] bL2[i] = 'X' good = True # now test them if bL2 in added: print k + str(i+1) print boardRep(bL2) print 'is a dup' continue for p in allPermutations(bL2): if p in added: print k + str(i+1) print boardRep(bL2) print 'is a perm of' print find_kforp(p) print boardRep(p) good = False break if good: print k + str(i+1) print boardRep(bL2) print 'is accepted\n' added.append(bL2) mD[k + str(i+1)] = bL2 print len(mD.keys())

tic-tac-toe: utilities

Let's continue with what is admittedly a rather silly example: a simulation of tic-tac-toe boards with the ultimate goal of analyzing strategy. This post just gives some code for utility functions. For example, there are functions to give all rotations and reflections for a board.

There are other functions (not yet debugged properly) for assessing board status. If you want to follow along at home, save the code as TTTUtils.py on your Desktop so you can import it in the following steps.

Here is the output, along with comments to tell what is going on:

# start 1 2 3 4 5 6 7 8 9 # rotations (left) 3 6 9 2 5 8 1 4 7 9 8 7 6 5 4 3 2 1 7 4 1 8 5 2 9 6 3 # reflections 7 8 9 4 5 6 1 2 3 3 2 1 6 5 4 9 8 7 1 4 7 2 5 8 3 6 9 9 6 3 8 5 2 7 4 1

[Update: I added several functions (multipleReps, winners, doubleThreat, etc.)]
Here's the code:

rows = [ [0,1,2], [3,4,5], [6,7,8] ] cols = [ [0,3,6], [1,4,7], [2,5,8] ] diags = [ [0,4,8], [2,4,6] ] corners = [0,2,6,8] all = rows + cols + diags def boardRep(bL): rL = list() for iL in [[0,1,2],[3,4,5],[6,7,8]]: line = list() for i in iL: e = bL[i] if e: line.append(e) else: line.append('.') rL.append(' '.join(line)) return '\n'.join(rL) + '\n' def multipleReps(mD,kL): repL = list() for k in kL: repL.append(boardRep(mD[k])) pL = [[],[],[]] for rep in repL: L = rep.strip().split('\n') for i,line in enumerate(L): pL[i].append(line) L = [k.ljust(11) for k in kL] pL.append(''.join(L)) for i in range(3): pL[i] = (' '*4).join(pL[i]) return '\n'.join(pL) + '\n' R = range(1,10) def rotateLeft(bL): # easier to think in 1..9 D = {1:3,2:6,3:9,4:2,5:5, 6:8,7:1,8:4,9:7} # keep lists indexed in 0..8 return [bL[D[i]-1] for i in R] def threeRotations(bL): rL = list() for i in range(3): bL = rotateLeft(bL) rL.append(bL) return rL def switch(bL,i,j): bL[i],bL[j] = bL[j],bL[i] def fourReflections(L): rL = [ L[6:] + L[3:6] + L[:3] ] bL = L[:] for i,j in zip((0,3,6),(2,5,8)): switch(bL,i,j) rL.append(bL) bL = L[:] for i,j in zip((1,2,5),(3,6,7)): switch(bL,i,j) rL.append(bL) bL = L[:] for i,j in zip((1,0,3),(5,8,7)): switch(bL,i,j) rL.append(bL) return rL def allPermutations(bL): rL = threeRotations(bL) rL.extend(fourReflections(bL)) return rL # find a line that can win def winners(bL,who): #print 'test winners for', who #print boardRep(bL) rL = list() for line in all: vL = [bL[i] for i in line] vL.sort() if vL == [None, who, who]: for i in line: if not bL[i]: break # return a list of (i, line) rL.append((i,line)) #if rL: print 'returning', rL,'\n' return rL def doubleThreat(bL,p,v=False): if v: print 'doubleThreat DT' iL = openPos(bL) rL = list() for i in iL: if v: print 'testing', i bL2 = bL[:] bL2[i] = p if v: print boardRep(bL2) wL = winners(bL2,who=p) if wL and len(wL) > 1: if v: print 'DT saving', if v: print i rL.append(i) if v: print 'DT returning', rL return rL def promisingLines(bL,who): #print 'promisingLine', bL, who rL = list() for line in all: vL = [bL[i] for i in line] if vL.count(who) == 1: if vL.count(None) == 2: for i in line: if not i: break rL.append((i,line)) return rL def wonPosition(bL): for line in all: vL = [bL[i] for i in line] if vL.count(None) == 0: p = vL[0] if vL.count(p) == 3: return p def noPossibleWin(bL): for line in all: vL = [bL[i] for i in line] vL = [e for e in vL if e] if len(vL) == 1: return False if len(vL) == 2: if vL[0] == vL[1]: return False return True def openPos(bL): return [i for i in range(len(bL)) if not bL[i]] def closedPos(bL): return [i for i in range(len(bL)) if bL[i]] def findKforBoard(bL): for k in mD: if mD[k] == bL: return k if __name__ == '__main__': bL = list() for i in R: bL.append(str(i)) print boardRep(bL) for e in threeRotations(bL): print boardRep(e) for e in fourReflections(bL): print boardRep(e)

Enumerated theory of tic-tac-toe

This post will just set up the problem I have in mind.

Symmetry plays an important role in tac-tac-toe. For example, there are only 3 first moves: center, corner, and side. (Let's say X plays first). We indicate the move by the square number (1-based indexing, for a change).

. . . X . . . X . . X . . . . . . . . . . . . . . . . 5 1 2

Starting from the center, there are 2 second moves, while from the corner or side, there are 5 second moves. That's because 16 = 18 and 12=14 and so on. A total of 12 positions:

. . . O . . . O . . X . . X . . X . . . . . . . . . . 5 51 52 ------------------------------------------- X . . X O . X . O X . . . . . . . . . . . . O . . . . . . . . . . . . . 1 12 13 15 X . . X . . . . O . . . . . . . . O 16 19 ------------------------------------------- . X . O X . . X . . X . . . . . . . O . . . O . . . . . . . . . . . . . 2 21 24 25 . X . . X . . . . . . . O . . . O . 27 28

I see that 12 and 21 look identical except that X and O are switched, but I will keep them distinct, to remember the difference in order of play between player 1 and player 2.

For the third move, it is still feasible to enumerate all the moves (although I'll admit it was quite tedious). Those starting positions with left-right (reflection) symmetry (51, 52, 15, 19, 25, 28) or rotational symmetry( ) have 4 additional possibilities. The rest have 7. We generate a total of 6 * 4 + 6 * 7 = 66 boards with 3 plays (but see below).

O . . O X . O . X O . . O . . . X . . X . . X . . X X . X . . . . . . . . . . . . . . . X 51 512 513 516 519 ------------------------------------------------------- . O . X O . . O . . O . . O . . X . . X . X X . . X . . X . . . . . . . . . . X . . . X . 52 521 524 527 528 ------------------------------------------------------- X O . X O X X O . X O . X O . . . . . . . X . . . X . . . X . . . . . . . . . . . . . . . 12 123 124 125 126 X O . X O . X O . . . . . . . . . . X . . . X . . . X 127 128 129 ------------------------------------------------------- X . O X X O X . O X . O X . O . . . . . . X . . . X . . . X . . . . . . . . . . . . . . . 13 132 134 135 136 X . O X . O X . O . . . . . . . . . X . . . X . . . X 137 138 139 ------------------------------------------------------- X . . X X . X . X X . . X . . . O . . O . . O . . O X . O . . . . . . . . . . . . . . . X 15 152 153 156 159 ------------------------------------------------------- X . . X X . X . X X . . X . . . . O . . O . . O X . O . X O . . . . . . . . . . . . . . . 16 162 163 164 165 X . . X . . X . . . . O . . O . . O X . . . X . . . X 167 168 169 ------------------------------------------------------- X . . X X . X . X X . . X . . . . . . . . . . . . X . . . X . . O . . O . . O . . O . . O 19 192 193 195 196 ------------------------------------------------------- O X . O X X O X . O X . O X . . . . . . . X . . . X . . . X . . . . . . . . . . . . . . . 21 213 214 215 216 O X . O X . O X . . . . . . . . . . X . . . X . . . X 217 218 219 ------------------------------------------------------- . X . X X . . X X . X . . X . O . . O . . O . . O X . O . X . . . . . . . . . . . . . . . 24 241 243 245 246 . X . . X . . X . O . . O . . O . . X . . . X . . . X 247 248 249 ------------------------------------------------------- . X . X X . . X . . X . . X . . O . . O . X O . . O . . O . . . . . . . . . . X . . . X . 25 251 254 257 258 ------------------------------------------------------- . X . X X . . X X . X . . X . . . . . . . . . . X . . . X . O . . O . . O . . O . . O . . 27 271 273 274 275 . X . . X . . X . . . X . . . . . . O . . O X . O . X 276 278 279 ------------------------------------------------------- . X . X X . . X . . X . . X . . . . . . . X . . . X . . . . . O . . 0 . . 0 . . 0 . X 0 . 28 281 284 285 287

As you probably have noticed, some of these are the same, just reflected or rotated. Just considering the 5.. series, I found six duplicates:

512 = 215 513 = 135 519 = 195 521 = 125 524 = 245 528 = 285

And notice that the other permutation is not the same, for example:

O X . O X . X X . . X . . X . . O . . . . . . . . . . 512 215 125

It only works if we switch X's, the first and third positions. We would have expected 132 = 231, but 231 is not in our list because 23 is the same as 21. I found two more:

152 = 251 192 = 273

So it looks like we have 66 - 8 = 58 unique positions.

From this point, there are 6! positions possible for each, neglecting the chance that some of these may be rotational isomers. So there is at most a total of 58 * 720 = 41760 possible positions, which is almost an order of magnitude less than 9! = 362880.

Actually, there are significantly fewer than that, since a number of the 3-mers force the next move by 'O' in order to block the win by 'X'. Also, some have symmetry, allowing only 3 possibilities instead of 6.

From the summary below (D = duplicate, S = symmetric, F = forced), I count 8 duplicates (removed), 23 forced moves (1 each), 9 with symmetry admitting 3 moves, and 1 with special symmetry (159) admitting only 2 moves. That leaves 66 - (8 + 23 + 9 + 1) = 66 - 41 = 25 with 6 possible moves. That's a total of 23 + 9 * 3 + 2 + 25 * 6 = 23 + 27 + 2 + 150 = 202 positions at the 4-mer stage. And some of these may be duplicates.

It's about now that I start begging for a simulation. I'll see what I can do.

Duplicate, Symmetric, Forced 123 S 124 F 125 F 126 127 F 128 129 F 132 134 F 135 F 136 137 F 138 139 F 152 153 S 156 159 S* 162 163 164 165 167 S 168 169 192 F 193 F 195 S 196 213 214 S 215 F 216 217 218 F 219 D 241 F 243 F 245 F 246 247 248 F 249 251 D 254 S 257 258 S 271 F 273 F 274 275 F 276 278 F 279 S 281 F 284 285 S 287 512 D 513 D 516 F 519 D 521 D 524 D 527 F 528 D

Wits

I recently read Stan (Stanislaw) Ulam's autobiography titled "Adventures of a Mathematician." Ulam was a Polish mathematician who worked on the Manhattan project and stayed at Los Alamos labs for much of his career. He was involved in the invention of the Monte Carlo method, supposedly named because he had an uncle who borrowed money to gamble at Monte Carlo.

Ulam tells some funny stories, many of them involving his friends Johnny von Neumann and Enrico Fermi. Here is a taste of Fermi's wit (p.164):

Once Segrè, who was very fond of fishing on weekends in the streams of the Los Alamos mountains, was expounding on the subtleties of the art, saying that it was not easy to catch trout. Enrico, who was not a fisherman, said with a smile, "Oh, I see Emilio, it is a battle of wits."

There is a wonderful picture of Ulam on the web (here).

Duly quoted

“It’s like they’re coming in and saying to you, ‘I’m going to drive my car off a cliff. Should I or should I not wear a seatbelt?’ And you say, ‘I don’t think you should drive your car off the cliff.’ And they say, ‘No, no, that bit’s already been decided – the question is whether to wear a seatbelt.’ And you say, ‘Well, you might as well wear a seatbelt.’ And then they say, ‘We’ve consulted with policy expert Rory Stewart and he says ...’”

Read it here. h/t Kevin Drum

I have to read the whole interview, Rory Stewart sounds like a very interesting guy.

The birthday problem

Everyone knows somebody who has the same birthday as they do. And if not, you can go to wikipedia, like I did, to find that both Anton van Leeuwenhoek and Kevin Kline were also born on October 24. In this post, I want to explore this famous problem. If we ignore Feb 29 and assume that births are evenly distributed over the days of the year (which may not be true), then the probability that two individuals chosen at random share the same birthday is 1/365.

So let's think about this famous gathering, and ask the question: if the individual chance that Albert and Max share the same birthday is 1/365, what is the probability that there are at least two people in this group that do share the same special day?



The key, as you probably guess, is that this is a combinations problem. Consider a group of 5 individuals (A-E), if F walks up to the group and introduces herself, there are 5 introductions to make.



If we think about building up a group in steps, then for n individuals the number of introductions that have been made is:

Σ 1 + 2 ... + n-1

This is the problem supposedly solved by Gauss as a young boy.

Another way to think about it is to consider that if each person in the group of n people shakes hands with every other person, there are n * (n-1) hands involved in handshakes, but then we must divide by two to get the unique interactions, since we've counted one hand for both of the interacting partners.

In general, we have the formula for combinations: C(n,k) = n! / (n-k)! k!, where k = 2.

We can solve the birthday problem in a couple of ways. We may say that we have the probability for each pair that they do not share a birthday, which is 364/365. The probability that all the independent combinations do not share a birthday is (364/365)**C(n,2). The probability of the complementary event, that at least one pair does share a birthday, is 1 - (364/365)**C(n,2).

The second approach is to consider the group with 2 people and P = 364/365 that they do not share a birthday. If a new person walks up to the group, there are 363 birthdays which would preserve the "no shared birthday" criterion. The probability of the desired event is then 1 - 364/365 * 363/365..., extended for n-1 steps.

This is easy to program, and I won't bore you with the details, but I will show this pretty plot I made using R:



Now, to the point of the post. I thought about finding a group near the critical size and testing it for the birthday criterion. I'm not such a big sports fan anymore, unless you consider politics a sport. How about the Presidents of the United States? Barack Obama is #44. You can get their vitals from wikipedia, but I found a text version on the web here.

The data needs just a bit of cleanup. One date lacks the comma, one date is listed as April 28th. And one entry has two tabs separating the name and the birthday. And why, exactly, are we presented with Obama's middle name? ("I got my middle name from someone who obviously, never thought I'd be running for President")-video.

Here are the results:

James K. Polk November 2 Warren G. Harding November 2 Millard Fillmore January 7 Richard M. Nixon January 9 William McKinley January 29 Franklin D. Roosevelt January 30 Ronald Reagan February 6 William Henry Harrison February 9 Abraham Lincoln February 12 George Washington February 22 Andrew Jackson March 15 James Madison March 16 Grover Cleveland March 18 John Tyler March 29 Thomas Jefferson April 13 James Buchanan April 23 Ulysses S. Grant April 27 James Monroe April 28 Harry S Truman May 8 John Kennedy May 29 George H. W. Bush June 12 Calvin Coolidge July 4 George W. Bush July 6 John Quincy Adams July 11 Gerald R. Ford July 14 Barack Hussein Obama August 4 Herbert Hoover August 10 William J. Clinton August 19 Benjamin Harrison August 20 Lyndon B. Johnson August 27 William Howard Taft September 15 Jimmy Carter October 1 Rutherford B. Hayes October 4 Chester A. Arthur October 5 Dwight D. Eisenhower October 14 Theodore Roosevelt October 27 John Adams October 30 Warren G. Harding November 2 James K. Polk November 2 James A. Garfield November 19 Franklin Pierce November 23 Zachary Taylor November 24 Martin Van Buren December 5 Woodrow Wilson December 28 Andrew Johnson December 29

And here is the code:

fn = 'presidents.txt' FH = open(fn,'r') data = FH.read().strip() FH.close() L = data.split('\n') def parseDate(s): s = s.strip() y = int(s.split()[-1]) s = s.split(',')[0] m,d = s.split() return {'year':y,'month':m, 'day':int(d) } D = dict() for e in L: t = e.split('\t') name,date = t[0],t[-1] # extra tabs in some D[name] = parseDate(date) for i,k1 in enumerate(D.keys()): for j in range(i): k2 = D.keys()[j] if D[k1]['month'] == D[k2]['month']: if D[k1]['day'] == D[k2]['day']: print k1.ljust(25), print D[k1]['month'],D[k1]['day'] print k2.ljust(25), print D[k2]['month'],D[k2]['day'] print mL = ['January','February','March','April', 'May','June','July','August','September', 'October','November','December'] def f(k): return mL.index(D[k]['month']),D[k]['day'] for k in sorted(D.keys(),key=f): print k.ljust(25), print D[k]['month'],D[k]['day']

the full Monty

There is one crucial point that I didn't make clear in the last post about the Monty Hall problem. If the probability that the other unopened door is the door with the prize changes after the host's action, information must have been received somehow. It comes from the fact that Monty is a "knowledgeable host"---he knows which door hides the prize, and he always opens a door that reveals a goat.

This contingency is made explicit by considering other potential host behaviors as described in the Wikipedia entry:

• Monty from Hell
• Angelic Monty
• Ignorant Monty
• Monty only offers sometimes

One good way to see that the result is correct is to extend the problem to a deck of cards. Suppose you are to choose among 52 cards, hoping to get the Ace of Spades.



You choose one card, which remains hidden, and now I turn over 50 cards, none of which turns out to be the Ace of Spades. It is pretty clear the probability that your first choice was correct is 1 in 52 and now the odds for the one card remaining are obviously much improved.

Monty, Monty...

By now, I'm sure you know about the "Monty Hall problem." It is a wonderful problem because many people, even those knowledgeable about statistics, find it difficult to believe the correct answer. If you don't know the story, here is wikipedia.

The short version:

There are 3 doors, behind one is a prize and behind the other two are goats. You first choose a door, which remains closed. The host must now open one of the other two doors. He does so, and behind this door is a goat. At this point, the host offers you the possibility of changing your choice to the third door. Should you switch?

The intuitive answer is that since there are two unopened doors, and ostensibly no information, they are equally logical choices. But this is not correct. For a detailed discussion, see Grinstead and Snell (example 4.6) or Krauss and Wang (2003 J. Exp. Psychol.: General. 132:3; pdf available for both on Wikipedia).

I wrote a Python simulation for the problem. Here is the output:

p = A c = A m = B p = B c = B m = A p = C c = B m = C p = A c = C m = A p = B c = C m = B p = B c = A m = B p = A c = C m = A p = C c = B m = C p = B c = B m = C p = C c = B m = C p = A c = C m = A p = A c = B m = A stay: 3332 switch: 6668

Here is a nice Java applet with a simulation.

And here is a syntax-colored screenshot of my Python code:

Simpson's paradox

In Dennis Lindley's book, Understanding Uncertainty, I came across a striking statistical paradox I had never heard about but which seems to be well known to statisticians of all ages. It's called Simpson's Paradox. The wikipedia entry has such a good example that I will shamelessly appropriate it:

Here are the batting averages for two major league baseball players for consecutive seasons in the 90's:

1995 1996 A .250 (12/48) .314 (183/582) B .253 (104/411) .321 (45/140) Combined A .310 (195/630) B .270 (149/551)

Player A's average was much higher considering the two seasons together, but in each individual year, player B had better numbers.

The key to understanding how this works is to notice that the number of attempts in individual years varies markedly, with player A having just 48 attempts in 1995 and player B having only 140 in 1996. This makes the comparisons for individual years depend on results that do not contribute so much to the combined total.

Remarkably, the situation continues even in the third season, although the combined totals are now pretty close:

1997 Combined A .291 (190/654) .300 (385/1284) B .329 (163/495) .298 (312/1046) Player A: Derek Jeter Player B: David Justice

In a medical situation, one could get a similar result with A and B being different treatments, and the years corresponding to some "lurking variable" (sex is a common example). The moral of the story: if you wish to break down results by some factor, you need to have similar sample sizes from the two groups.

Here is another striking example involving vectors, which I generated using R code from wikipedia.



Individual trials 1 and 2 show that red has a higher "success rate", but when the two trials are combined blue is better. There are other great examples in the article.

Code simplicity

This is hilarious:

If somebody comes up to you and says something like, “How do I make this pony fly to the moon?”, the question you need to ask is, “What problem are you trying to solve?” You’ll find out that they really need to collect gray rocks. Why they thought they had to fly to the moon, and use a pony to do it, only they know. People do get confused like this.

Tic Tac Toe

This project is a rewrite of one I posted here. I don't want to show all the code, but I put the zipped files for the project (including the nib file) on the .mac server. The buttons have an image assigned depending on which player chose that square. This is done with bindings, but the buttons are bound individually to variables in the DisplayController, so there is a variable for each button. I've been trying to figure out a way to do this more elegantly, but no luck yet. I used one-based indexing for the board, which accounts for stuff like L = range(1,10). Here is the ugly code for that class:

class DisplayController(NSObject): myBoard = objc.IBOutlet() sq1 = objc.ivar('sq1') sq2 = objc.ivar('sq2') sq3 = objc.ivar('sq3') sq4 = objc.ivar('sq4') sq5 = objc.ivar('sq5') sq6 = objc.ivar('sq6') sq7 = objc.ivar('sq7') sq8 = objc.ivar('sq8') sq9 = objc.ivar('sq9') def setsq1_(self,value): self.sq1 = value def setsq2_(self,value): self.sq2 = value def setsq3_(self,value): self.sq3 = value def setsq4_(self,value): self.sq4 = value def setsq5_(self,value): self.sq5 = value def setsq6_(self,value): self.sq6 = value def setsq7_(self,value): self.sq7 = value def setsq8_(self,value): self.sq8 = value def setsq9_(self,value): self.sq9 = value def init(self): self.bush = NSImage.imageNamed_('bush') self.chimp = NSImage.imageNamed_('chimp') self.blank = NSImage.imageNamed_('blank') return self def update_(self,sender,which): fL = [None, self.setsq1_,self.setsq2_,self.setsq3_, self.setsq4_,self.setsq5_,self.setsq6_, self.setsq7_,self.setsq8_,self.setsq9_] if which is not 'all': L = [which] else: L = range(1,10) for i in L: marker = self.myBoard.b[i] if not marker: fL[i](self.blank) elif b marker == 'P': fL[i](self.chimp) else: fL[i](self.bush)

Fun with Sudoku

In a previous life, I was interested in Sudoku for a while. I never got into doing the puzzles, but I wanted to write a Sudoku solver in Python. It seemed like a fun problem, and it is not too hard, until you get to puzzles which have no two-way decision points where one decision leads to an invalid puzzle, forcing the other path. Then I stumbled across this. Naturally, Peter Norvig is way smarter than me, that's one reason Google pays him so much. If you don't believe me, watch this. Norvig says: "as computer security expert Ben Laurie has stated, Sudoku is "a denial of service attack on human intellect." That sounds about right to me.

Anyway, while I was thinking about Sudoku, I read the Wikipedia entry. It says: "The numerals...are used for convenience...any set of distinct symbols will do, such as letters, shapes or colours." And that got me thinking... I ended up writing a Cocoa application for Color Sudoku. At the time (late 2005), nobody else had implemented one that I could find on the web. However, a quick Google search shows that the situation has changed.

Still, I don't think anybody has done it like I did. It is so long ago that I really don't feel like reading through the code. Like many projects, it grew without apparent direction. But the app still works under Leopard and it is fun to play. It is in my public folder at .Mac. There are other goodies there as well. Here is a screenshot:



Each number is represented by a different color. If a given position in the 9 x 9 grid is completely determined, it is filled in with a solid color. Otherwise, the possible values for the square are given as smaller squares. You click on squares to make them go away, when you deduce that they are not possible values (or command-click to select one). I find it very easy and intuitive to play, and it beats the heck out of crossing out numbers.

I notice that the version in the folder does not bring the window back with command-N. I only learned how to do that later.

One more point. There is a lot of Perl used for Bioinformatics. As far as I can tell, there is no reason to write any new project in Perl. I know it is partly a matter of taste. Still, this guy is very proud that his Perl Sudoku solver is only four lines of completely obscure code (actually, I see he has made it three lines now). Not all Perl is like this, but I found that after a week or two, I simply could not understand how my programs worked without a lot of effort. It's not like that in Python. Here are a couple of citations on point:

Programs must be written for people to read, and only incidentally for machines to execute.

- Abelson & Sussman, SICP

How do we convince people that in programming simplicity and clarity—in short: what mathematicians call "elegance"— are not a dispensable luxury, but a crucial matter that decides between success and failure?

- E. W. Dijkstra