--- a/.lispwords Wed Nov 07 20:39:43 2018 -0500
+++ b/.lispwords Thu Nov 08 19:23:33 2018 -0500
@@ -1,3 +1,5 @@
(1 codon-case)
(1 labels-memoized)
(4 define-problem)
+(1 do-sum Σ)
+(1 do-product Π)
--- a/rosalind.asd Wed Nov 07 20:39:43 2018 -0500
+++ b/rosalind.asd Thu Nov 08 19:23:33 2018 -0500
@@ -27,22 +27,27 @@
(:module "src" :serial t
:components ((:file "utils")
(:module "problems"
- :components ((:file "dna")
- (:file "rna")
- (:file "revc")
+ :components (
+
+ (:file "cons")
+ (:file "dna")
+ (:file "fib")
+ (:file "fibd")
(:file "gc")
+ (:file "grph")
(:file "hamm")
- (:file "prot")
- (:file "perm")
- (:file "fib")
- (:file "subs")
+ (:file "iev")
(:file "iprb")
- (:file "iev")
- (:file "fibd")
- (:file "cons")
- (:file "grph")
+ (:file "lcsm")
+ (:file "lia")
+ (:file "mrna")
+ (:file "perm")
+ (:file "prot")
(:file "prtm")
- (:file "mrna")
+ (:file "revc")
+ (:file "rna")
(:file "splc")
- (:file "lcsm")))))))
+ (:file "subs")
+ ))))))
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/problems/lia.lisp Thu Nov 08 19:23:33 2018 -0500
@@ -0,0 +1,109 @@
+(in-package :rosalind)
+
+;; When a heterozygous organism mates, its offspring have a 50% chance to be
+;; heterozygous themselves, *regardless of what the other mate happens to be*:
+;;
+;; A a A a A a
+;; A AA Aa A AA Aa a aA aa
+;; A AA Aa a Aa aa a aA aa
+;;
+;; Because everyone breeds with an Aa Bb mate, we don't need to worry about
+;; tracking populations along the way. Every child has a 1/2 chance of being
+;; Aa, and likewise a 1/2 chance of being Bb.
+;;
+;; Because A's and B's are independent, this means there's a 1/2 * 1/2 = 1/4
+;; chance of any given child being Aa Bb.
+
+(defmacro do-sum ((var from to) &body body)
+ "Sum `body` with `var` iterating over `[from, to]`.
+
+ It's just Σ:
+
+ to
+ ===
+ \
+ > body
+ /
+ ===
+ n = from
+
+ "
+ (once-only (to)
+ (with-gensyms (result)
+ `(do ((,var ,from (1+ ,var))
+ (,result 0))
+ ((> ,var ,to) ,result)
+ (incf ,result (progn ,@body))))))
+
+(defmacro do-product ((var from to) &body body)
+ "Multiply `body` with `var` iterating over `[from, to]`.
+
+ It's just Π:
+
+ to
+ =====
+ | |
+ | | body
+ | |
+ n = from
+
+ "
+ (once-only (to)
+ (with-gensyms (result)
+ `(do ((,var ,from (1+ ,var))
+ (,result 1))
+ ((> ,var ,to) ,result)
+ (setf ,result (* ,result (progn ,@body)))))))
+
+
+(defmacro Σ (bindings &body body) ;; lol
+ `(do-sum ,bindings ,@body))
+
+(defmacro Π (bindings &body body) ;; lol
+ `(do-product ,bindings ,@body))
+
+
+(defun binomial-coefficient (n k)
+ "Return `n` choose `k`."
+ ;; https://en.wikipedia.org/wiki/Binomial_coefficient#Multiplicative_formula
+ (Π (i 1 k)
+ (/ (- (1+ n) i) i)))
+
+(defun bernoulli-exactly (successes trials success-probability)
+ "Return the probability of exactly `successes` in `trials` Bernoulli trials.
+
+ Returns the probability of getting exactly `n` successes out of `trials`
+ Bernoulli trials with `success-probability`.
+
+ "
+ ;; For a group of N trials with success/failure probabilities s/f, any given
+ ;; ordering of exactly S successes and F failures will have probability:
+ ;;
+ ;; s * s * s * … * f * f * f
+ ;;
+ ;; There are N choose S possible orderings (or N choose F, it doesn't matter),
+ ;; so we just sum up the probabilities of all of them.
+ (let ((failures (- trials successes))
+ (failure-probability (- 1 success-probability)))
+ (* (binomial-coefficient trials successes)
+ (expt success-probability successes)
+ (expt failure-probability failures))))
+
+(defun bernoulli-at-least (successes trials success-probability)
+ "Return the probability of at least `successes` in `trials` Bernoulli trials.
+
+ Returns the probability of getting at least `n` successes out of `trials`
+ Bernoulli trials with `success-probability`.
+
+ "
+ ;; P(≥S) = P(=S) + P(=S+1) + … + P(N)
+ (Σ (n successes trials)
+ (bernoulli-exactly n trials success-probability)))
+
+(define-problem lia (data stream)
+ "2 1"
+ "0.684"
+ (let* ((generations (read data))
+ (target (read data))
+ (population (expt 2 generations)))
+ (format nil "~,3F" (bernoulli-at-least target population 1/4))))