A Monte Carlo analysis of a problem by Euler.

euler-population.clj

(ns euler-population.core)

;; "After the Flood, all people are descended from six people. If we
;; suppose that after 200 years the population had grown to 1,000,000,
;; then by what part must the population grow each year?"
;; — Leonhard Euler's "Introductio in Analysin Infinitorum", Ch. 6, example 3

(def initial-conditions {:men 3
                         :women 3})

(def ^:dynamic *reproduction-rate* 0.90)
(def ^:dynamic *polygamy* false)

(defn compute-reproduction-rate 
  ([iter-limit] (compute-reproduction-rate 1e6 iter-limit))
  ([desired-pop iter-limit]
    (let [C (* 2 (Math/pow (double (/ desired-pop 6)) (double (/ 1 iter-limit))))]
      (/ C (+ 1 C)))))

(defn- avg [coll]
  (when (seq coll)
    (double (/ (reduce + 0 coll) (count coll)))))

(defn- family-size [{men :men women :women}] (+ men women))

(defn number-of-couples [state]
  (range (if *polygamy*
           (:women state)
           (min (:men state) (:women state)))))

(defn- have-child [] (if (> (rand) 0.5) :men :women))

(defn generate-family []
  (loop [family {:men 0 :women 0}]
    (if (>= (rand) *reproduction-rate*)
      family
      (recur (update-in family [(have-child)] inc)))))

(defn next-generation [state]
  (let [families (number-of-couples state)]
    (if (empty? families)
      (throw (Exception. "No more families..."))
      (apply (partial merge-with +)
             (for [couple families]
               (generate-family))))))

(defn run-for [N-generations]
  (try
    (loop [state initial-conditions
           k 0]
      (if (>= k N-generations)
        (family-size state)
        (recur (next-generation state) (inc k))))
    (catch Exception e 0)))