monte-carlo-grouping.clj

(defn mean
  "The aritmetic mean (average) of xs"
  [xs] 
  (double
   (/ (reduce + xs) 
      (count xs))))

(defn gaussians
  "An infinite sequence of gaussians about 0 with variance v"
  [v]
  (let [gen (java.util.Random.)]
    (lazy-seq
     (cons
      (* v (.nextGaussian gen))
      (gaussians variance)))))

(defn abs
  "The absolute value of x"
  [x]
  (if (> 0 x)
    (* x -1)
    x))

(defn about?
  "Returns true if y is within r distance from x"
  [x d y]
  (let [d2 (abs d)
        h (+ x d2)
        l (- x d2)]
    (and
     (>= h y)
     (<= l y))))

(defn grab
  "The items within r distance from n; removes them from the-items ref"
  [n r the-items]
  (dosync
   (let [items @the-items
         in-range? #(about? n r %)]
     (alter the-items #(remove in-range? %))
     (filter in-range? items))))

(defn group
  "Groups of items of average size avg"
  [items avg]
  (let [variance (+ 1 (* 2 avg))
        the-items (ref items)
        r (first (gaussians variance))
        group-fns (for [p items] #(grab p r the-items))
        groups (apply pcalls group-fns)]
    (remove empty? groups)))

(def small-items
  [5 15 16 18 24])

(def big-items
  [392 396 400 406 422 436 446 448 450 454 462
   470 478 488 490 508 512 518 526 528 532 538
   548 562 570 580 592 594 596 598 600 604 608
   614 620 626 632 640])

(defn simulate [times]
  "Run a simulation, get back the mean of average groups"
  (mean
   (apply pcalls
          (for [n (range times)]
            (let [groups (group big-items 2)]
              #(mean (map count groups)))))))

(defn print-groups
  "Print groups n times"
  [items n]
  (dotimes [m n]
    (println
     (group items 2))))