Monte Carlo Elisp throwaways

coin-flip.el

;; 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

;; 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

;;; 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)