public
Created

Monte Carlo Elisp throwaways

  • Download Gist
coin-flip.el
Emacs Lisp
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
;; 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)
ibm-puzzle.el
Emacs Lisp
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
;; 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
Emacs Lisp
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
;;; 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)

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.