Instantly share code, notes, and snippets.

skeeto/coin-flip.el Created Jul 30, 2012

What would you like to do?
Monte Carlo Elisp throwaways
 ;; Minimum Number Of Coin Tosses Such That Probability Of Getting 3 ;; Consecutive Heads Is Greater Than 1/2 ;; http://blog.eduflix.tv/2012/05/contest-answer-minimum-number-of-coin-tosses-such-that-probability-of-getting-3-consecutive-heads-is-greater-than-12/ (require 'cl) (defun flip () "Flip a coin." (if (< 0 (random)) 'H 'T)) (defun flipn (n) "Generate a list of N flips." (if (= n 0) () (cons (flip) (flipn (- n 1))))) (defun* count-run (flips &optional (count 0) (best 0)) "Count the longest length run of heads in FLIPS." (if (null flips) best (let ((ncount (if (eq 'H (car flips)) (1+ count) 0))) (count-run (cdr flips) ncount (max best ncount))))) (defun trial (ntrials nflips) "Run a Monte Carlo for NFLIPS coin flips." (let ((match 0)) (dotimes (i ntrials) (if (>= (count-run (flipn nflips)) 3) (incf match))) (/ match 1.0 ntrials))) ;; Create a results table (let ((n 1000000)) (map 'vector (apply-partially 'trial n) (number-sequence 1 20))) ;;; 1,000,000 trials ;; 1 0.0 ;; 2 0.0 ;; 3 0.124649 ;; 4 0.186849 ;; 5 0.250329 ;; 6 0.312954 ;; 7 0.36751 ;; 8 0.417768 ;; 9 0.465094 ;; 10 0.508078 ;; 11 0.547962 ;; 12 0.583849 ;; 13 0.617091 ;; 14 0.647718 ;; 15 0.676534 ;; 16 0.702703 ;; 17 0.725677 ;; 18 0.748897 ;; 19 0.768097 ;; 20 0.787552 ;;; Actual solution (defun T (n) (if (<= n 1) 1 (+ (H (- n 1)) (T (- n 1))))) (defun H (n) (if (<= n 1) 1 (+ (T (- n 1)) (T (- n 2))))) (defun p (n) (- 1 (/ (+ (T n) (H n)) (expt 2.0 n)))) (memoize 'T) (memoize 'H) (map 'vector 'p (number-sequence 1 60)) ;;; Junk (setq eval-expression-print-length 1000) (setq eval-expression-print-level nil) (mapatoms 'byte-compile)
 ;; Problem: Alice and Bob are playing two games: they start ;; simultaneously and the first one to win his game is the winner. ;; ;; Alice is given an urn with N balls, colored in N different colors ;; and in every second she randomly picks two balls, colors the first ;; ball in the color of the second ball and returns both to the urn. ;; ;; Her game is over once all the balls in the urn have the same color. ;; ;; Bob is given an urn with M balls, all colorless, and in every ;; second he picks a random ball, color it, and puts it back to the ;; urn. Bob’s game is over once all his balls are colored. ;; ;; Our question is: what are the possible values of M for which (on ;; average) Bob is expected to lose for N=80 and win for N=81? (require 'cl) ;; Support (defun uniformp (bowl color) (catch 'fail (dotimes (i (length bowl) t) (if (not (eq (aref bowl i) color)) (throw 'fail nil))))) ;; Bob (defun bob-init (size) (make-vector size nil)) (defun bob (bowl) (let ((count 0)) (while (not (uniformp bowl 't)) (aset bowl (random (length bowl)) 't) (incf count)) count)) ;; Alice (defun alice-init (size) (apply 'vector (number-sequence 1 size))) (defun alice (bowl) (let ((count 0)) (while (not (uniformp bowl (aref bowl 0))) (let ((a (random (length bowl))) (b (random (length bowl)))) (aset bowl a (aref bowl b)) (incf count))) count)) ;; Monte Carlo (defun sim (playerf initf n runs) (let ((mean 0)) (dotimes (i runs mean) (setq mean (+ mean (/ (funcall playerf (funcall initf n)) 1.0 runs)))))) (defun table (playerf initf min max) (let ((runs 4000) (table '())) (dotimes (i (- max min) (reverse table)) (let* ((n (+ i min)) (result (sim playerf initf n runs))) (setq table (cons (cons n result) table)))))) (defun tables (max) (list (cons 'alice (table 'alice 'alice-init 1 max)) (cons 'bob (table 'bob 'bob-init 1 max)))) ;; Run ;(insert "\n" (pp (tables 101))) (setq results '((alice (1 . 0.0) (2 . 2.002457) (3 . 6.105002) (4 . 11.825988) (5 . 20.008257) (6 . 29.953535) (7 . 42.47124) (8 . 56.0001) (9 . 70.63585) (10 . 91.07466) (11 . 109.917145) (12 . 131.69872) (13 . 155.38876) (14 . 182.4785) (15 . 212.73148) (16 . 239.16151) (17 . 270.24463) (18 . 310.75928) (19 . 347.11017) (20 . 377.62512) (21 . 418.27908) (22 . 464.14706) (23 . 506.0562) (24 . 552.4154) (25 . 596.2141) (26 . 642.0864) (27 . 703.79) (28 . 749.7222) (29 . 808.4281) (30 . 871.599) (31 . 936.3256) (32 . 999.0706) (33 . 1054.7343) (34 . 1111.2745) (35 . 1190.465) (36 . 1255.0137) (37 . 1331.1805) (38 . 1407.3064) (39 . 1487.6403) (40 . 1550.6205) (41 . 1618.6235) (42 . 1727.6727) (43 . 1813.0775) (44 . 1893.2147) (45 . 1992.5043) (46 . 2077.9153) (47 . 2144.497) (48 . 2277.776) (49 . 2382.2654) (50 . 2490.0784) (51 . 2493.3884) (52 . 2656.8643) (53 . 2769.7937) (54 . 2875.5928) (55 . 2981.537) (56 . 3019.7056) (57 . 3226.736) (58 . 3291.034) (59 . 3454.302) (60 . 3504.081) (61 . 3694.1597) (62 . 3802.1477) (63 . 3883.796) (64 . 4025.8975) (65 . 4263.0576) (66 . 4326.8696) (67 . 4345.3965) (68 . 4550.426) (69 . 4658.369) (70 . 4880.2505) (71 . 4971.99) (72 . 5086.295) (73 . 5238.4614) (74 . 5443.6943) (75 . 5543.5015) (76 . 5749.335) (77 . 5872.9556) (78 . 5913.406) (79 . 6128.497) (80 . 6365.5366) (81 . 6527.9507) (82 . 6561.5415) (83 . 6824.895) (84 . 6898.22) (85 . 7107.244) (86 . 7269.1426) (87 . 7417.496) (88 . 7704.0415) (89 . 7870.894) (90 . 8041.8257) (91 . 8271.157) (92 . 8341.924) (93 . 8488.021) (94 . 8780.235) (95 . 8891.476) (96 . 9049.977) (97 . 9291.529) (98 . 9524.37) (99 . 9734.146) (100 . 9812.312)) (bob (1 . 0.9999731) (2 . 3.0097127) (3 . 5.490466) (4 . 8.321182) (5 . 11.495453) (6 . 14.828194) (7 . 18.069202) (8 . 21.86045) (9 . 25.529926) (10 . 29.414255) (11 . 33.083313) (12 . 37.56758) (13 . 41.730797) (14 . 45.249268) (15 . 49.84833) (16 . 54.137608) (17 . 57.98371) (18 . 63.100243) (19 . 67.44546) (20 . 71.922554) (21 . 76.699974) (22 . 80.76666) (23 . 86.63025) (24 . 90.386055) (25 . 95.33459) (26 . 99.70857) (27 . 104.59158) (28 . 110.2528) (29 . 114.441154) (30 . 120.62691) (31 . 124.60795) (32 . 131.11934) (33 . 134.76894) (34 . 140.57611) (35 . 144.69154) (36 . 150.1664) (37 . 155.35785) (38 . 160.1864) (39 . 166.40111) (40 . 172.37294) (41 . 175.43832) (42 . 183.1514) (43 . 186.60425) (44 . 191.17455) (45 . 197.35544) (46 . 204.3226) (47 . 209.11603) (48 . 212.70107) (49 . 218.54588) (50 . 224.89168) (51 . 231.84097) (52 . 237.5957) (53 . 240.1086) (54 . 247.65561) (55 . 252.85236) (56 . 258.91718) (57 . 264.9629) (58 . 268.86707) (59 . 275.50586) (60 . 281.0956) (61 . 286.4899) (62 . 293.90005) (63 . 298.26697) (64 . 302.80377) (65 . 309.73608) (66 . 315.5479) (67 . 321.16452) (68 . 329.4876) (69 . 334.26077) (70 . 338.07364) (71 . 347.04492) (72 . 352.44537) (73 . 356.4692) (74 . 362.56384) (75 . 368.52573) (76 . 374.26178) (77 . 381.79062) (78 . 384.69492) (79 . 394.21405) (80 . 396.02576) (81 . 400.15216) (82 . 409.19498) (83 . 418.5294) (84 . 421.8215) (85 . 428.92725) (86 . 430.54443) (87 . 437.92624) (88 . 447.21964) (89 . 449.65555) (90 . 459.1126) (91 . 461.70978) (92 . 467.96185) (93 . 476.41946) (94 . 481.19977) (95 . 490.34436) (96 . 496.3833) (97 . 500.05154) (98 . 508.84027) (99 . 510.18802) (100 . 516.7757)))) (dolist (row (cdr (assoc 'bob results))) (insert (pp (car row)) " " (pp (cdr row)) "\n")) ;(table 'alice 'alice-init 1 90) ;(sim 'alice 'alice-init 4 4000) ;(alice (alice-init 80)) ;(bob (bob-init 20)) ;; (mapatoms 'byte-compile)