Skip to content

Instantly share code, notes, and snippets.

@skeeto
Created July 30, 2012 14:55
Show Gist options
  • Save skeeto/3207572 to your computer and use it in GitHub Desktop.
Save skeeto/3207572 to your computer and use it in GitHub Desktop.
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)
;;; marbles.el --- solve the marble problem
;; Problem: You have a bag with one white marble and two black
;; marbles. You remove a marble, note whether it was white or black,
;; then put it back in the bag. Repeat this process indefinitely. What
;; is the chance the number of white marbles procured – at some point
;; – surpasses the number of black marbles procured?
;;
;; More generally: As above, with m white marbles, n black marbles,
;; m <= n, gcd(m, n) = 1.
;;
;; http://redd.it/x91i3/
(require 'cl)
(require 'eieio)
;; Define a generic marble bag. Only add-marble, marble-types, draw,
;; and count-marbles know the internal structure of the object.
(defclass bag ()
((total :initform 0 :reader size)
(counts :initform ()))
:documentation "A bag of marbles.")
(defmethod add-marble ((bag bag) marble &optional count)
"Add a marble to the bag."
(if (null count) (setf count 1))
(let ((key (assoc marble (slot-value bag 'counts))))
(if (null key)
(setf (slot-value bag 'counts)
(cons (cons marble count) (slot-value bag 'counts)))
(incf (cdr key) count))
(incf (slot-value bag 'total) count)))
(defmethod count-marbles ((bag bag) marble)
"Return the number of MARBLE in the bag."
(cdr (assoc marble (slot-value bag 'counts))))
(defmethod marble-types ((bag bag))
"Return the types of marbles this bag contains."
(mapcar (lambda (key) (car key)) (slot-value bag 'counts)))
(defmethod draw ((bag bag))
"Draw a marble from the bag and replace it."
(let ((n (random (size bag))))
(catch 'found
(dolist (key (slot-value bag 'counts))
(if (< n (cdr key))
(throw 'found (car key))
(decf n (cdr key)))))))
(defmethod combine ((a bag) &rest bags)
"Create a new bag containing the contents of the given bags."
(let ((bag (make-instance 'bag)))
(dolist (old (cons a bags) bag)
(dolist (type (marble-types old))
(add-marble bag type (count-marbles old type))))))
;; Problem-specific functions
(defmethod initialize-instance ((bag bag) &rest args)
(if (= 2 (length (car args)))
(destructuring-bind ((white black)) args
(add-marble bag 'white white)
(add-marble bag 'black black))))
(defun simulate (bag timeout)
"Simulate a game using BAG, stopping at TIMEOUT."
(let ((white 0))
(catch 'overtake
(dotimes (n timeout nil)
(if (eq 'white (draw bag))
(incf white))
(if (> white (/ (1+ n) 2))
(throw 'overtake t))
(if (<= (+ white (- timeout n 1)) (/ timeout 2))
(throw 'overtake nil))))))
(defun monte-carlo (bag timeout trials)
"Use the Monte Carlo method compute the probability that white
overtakes black for this BAG."
(let ((overtake 0))
(dotimes (n trials (/ overtake 1.0 trials))
(if (simulate bag timeout)
(incf overtake)))))
(defun analyze (bag)
"Analytically compute the probability that white overtakes
black for this BAG."
(/ (count-marbles bag 'white) 1.0 (count-marbles bag 'black)))
(defun compare (white black)
"Compare the results of the Monte Carlo method to the
analytical, returning the percent error."
(let ((bag (make-instance 'bag white black))
(limit 4000))
(abs (* 100 (- (analyze bag) (monte-carlo bag limit limit))))))
;; Test
;(compare 3 10)
;=> 0.050000000000000044
;(setq a (make-instance 'bag 10 10))
;(add-marble a 'black)
;(add-marble a 'white)
;(add-marble a 'red)
;(add-marble a 'green)
;(size a)
;(draw a)
;(combine a (make-instance 'bag 10 10))
;(marble-types a)
;(monte-carlo a 1000 1000)
;(analyze a)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment