prob.lisp (4921B)
1 ;; Editor's note: Helper script written for a Nature of Reason (80-150) probability assignment. 2 ;; An underrated quality of s-expressions is that you get to have simple rational literals... 3 ;; Written around the spring of 2024. 4 5 (load "~/quicklisp/setup.lisp") 6 7 (defmacro defalias (to fn) `(setf (fdefinition ',to) #',fn)) 8 (defalias filter remove-if-not) 9 (defmacro fn (&body body) `(lambda () ,@body)) 10 (defmacro fnx (&body body) `(lambda (x) ,@body)) 11 12 (defun fact (x &optional (acc 1)) 13 (if (zerop x) acc 14 (fact (- x 1) (* acc x)))) 15 16 (defun choose (n k) 17 (/ (fact n) (* (fact k) (fact (- n k))))) 18 19 (defun generate (f n) 20 (loop for i from 0 to (- n 1) collect (funcall f))) 21 22 (defun generate-cond (f p n &optional (acc nil)) 23 (if (eq (length acc) n) acc 24 (let ((res (funcall f))) 25 (generate-cond f p n (if (funcall p res) (cons res acc) acc))))) 26 27 (defun exp/prob (sample pred &optional (samples 10000)) 28 "Estimates probability of P being true by repeatedly sampling space via SAMPLE." 29 (/ (length (filter pred (generate sample samples))) 30 (float samples))) 31 32 (defun exp/cond-prob (sample cond pred &optional (samples 10000)) 33 "Estimates probability of PRED being true given COND by repeatedly sampling space via SAMPLE." 34 (let ((cands (generate-cond sample cond samples))) 35 (/ (length (filter pred cands)) 36 (float (length cands))))) 37 38 (defun exp/compare (exp theory) 39 (abs (- exp theory))) 40 41 (defun thm/mutex-or (&rest args) (reduce '+ args)) 42 (defun thm/or (a b both) (- (+ a b) both)) 43 (defun thm/indp-and (a b) (* a b)) 44 (defun thm/indp-trials (trials successes psuccess) 45 (abs (* (choose trials successes) 46 (expt psuccess successes) 47 (expt (1- psuccess) (- trials successes))))) 48 (defun thm/cond (a given both) (/ both given)) 49 (defun thm/bayes (a given flipped) (/ (* flipped a) given)) 50 51 (defun thm/not (a) (- 1 a)) 52 53 (defun is-prime (n &optional (d (- n 1))) 54 (or (= d 1) 55 (and (/= (rem n d) 0) 56 (is-prime n (- d 1))))) 57 58 (defun flip-coin (&key (times 1) (heads-bias nil)) 59 (let ((f (if heads-bias (lambda () (if (< (random 1.0) heads-bias) 'heads 'tails)) 60 (lambda () (nth (random 2) '(heads tails)))))) 61 (if (eq times 1) (funcall f) (generate f times)))) 62 63 ;; PROBLEM 1 64 ;; a) 65 (exp/compare 66 (exp/prob (fn (1+ (random 6))) (fnx (or (oddp x) (is-prime x))) 1000000) 67 (thm/or 3/6 3/6 2/6)) 68 69 ;; b) 70 71 (exp/compare 72 (exp/cond-prob 73 (lambda () (flip-coin :times 11)) 74 (lambda (flips) (every (lambda (x) (eq x 'tails)) (subseq flips 0 9))) 75 (lambda (flips) (eq (car (last flips)) 'heads)) 1000) 76 (thm/cond 1/2 (thm/indp-trials 10 10 1/2) (thm/indp-trials 11 11 1/2))) 77 78 ;; c) 79 80 (exp/compare 81 (exp/prob 82 (fn (flip-coin :heads-bias 2/3 :times 2)) 83 (fnx (equal x '(heads heads)))) 84 (thm/indp-trials 2 2 2/3)) 85 86 (exp/compare 87 (exp/prob 88 (fn (flip-coin :heads-bias 2/3 :times 2)) 89 (fnx (equal x '(heads tails)))) 90 (thm/indp-and 2/3 1/3)) 91 92 ;; d) 93 (exp/compare 94 (exp/prob 95 (fn (flip-coin :times 3)) 96 (fnx (not (equal x '(heads heads heads))))) 97 (thm/not (thm/indp-trials 3 3 1/2))) 98 99 (exp/compare 100 (exp/prob 101 (fn (flip-coin :times 3)) 102 (fnx (some (fnx (eq x 'heads)) x))) 103 (thm/mutex-or (thm/indp-trials 3 1 1/2) 104 (thm/indp-trials 3 2 1/2) 105 (thm/indp-trials 3 3 1/2))) 106 107 108 (defun exp/frac (p) 109 (< (random 1.0) p)) 110 111 ;; Shitty imperative code. 112 (defun exp/weighted-choose (choices) 113 (let ((x (random 1.0)) (cmp 0) (out nil)) 114 (dolist (pair choices) 115 (incf cmp (cadr pair)) 116 (if (< x cmp) 117 (return-from exp/weighted-choose (car pair)))))) 118 119 ;; problem 2 120 ;; 121 ;; P(A|S) = 3/5 122 ;; P(S) = 4/5 123 ;; P(A|~S) = 1/9 124 ;; P(A)? 125 ;; P(A) = P(A|S)P(S) + P(A|~S)P(~S) by Bayes' 126 (exp/compare 127 (exp/prob (fn (if (exp/frac 4/5) (exp/frac 3/5) (exp/frac 1/9))) 128 (fnx (eq x t)) 1000000) 129 (+ (* 3/5 4/5) (* 1/9 1/5))) 130 131 ;; 2b 132 ;; 133 ;; P(B) = 0.2 134 ;; P(L|B) = 0.07 135 ;; P(L) = P(L|W)P(W) + P(L|B)P(B) + P(L|C)P(C) = 0.1*0.5 + 0.03*0.3 + 0.07*0.2 = 0.2 136 ;; P(B|L) = (P(L|B)P(B))/P(L) by Bayes' 137 138 (exp/compare 139 (exp/cond-prob 140 (fn (let ((transit (exp/weighted-choose '((walk 0.5) (car 0.3) (bus 0.2))))) 141 (list transit (case transit 142 (walk (exp/frac 0.1)) 143 (car (exp/frac 0.03)) 144 (bus (exp/frac 0.07)))))) 145 (fnx (eq (cadr x) t)) 146 (fnx (eq (car x) 'bus)) 10000) 147 (thm/bayes 2/10 148 (thm/mutex-or (thm/indp-and 1/10 5/10) 149 (thm/indp-and 3/100 3/10) 150 (thm/indp-and 7/100 2/10)) 151 7/100)) 152 153 ;; problem 3 154 155 ;; example: numbers from 0 to 50 (exclusive) 156 ;; event B: x is even 157 ;; event A: x is divisible by 10 158 159 (exp/prob 160 (fn (random 50)) 161 (fnx (or (not (zerop (mod x 2))) 162 (zerop (mod x 10))))) 163 164 (remove-if-not (fnx (zerop (mod x 10))) (loop for i from 0 to 49 collect i)) 165 166 (exp/cond-prob 167 (fn (random 50)) 168 (fnx (zerop (mod x 2))) 169 (fnx (zerop (mod x 10)))) 170