standard-sampling.clj

(defn log-add 
 "log (sum xs) from seq of log-x"
 [log-xs]
 (let [max-log-x (apply max log-xs)]
   (+ max-log-x
      (Math/log (sum          
          (for [log-x log-xs 
                 :let [diff (- log-x max-log-x)]
                 :when (> diff -30)]
             (Math/exp diff)))))))

(defn log-normalize [log-xs]
 (let [log-sum (log-add log-xs)]
   (map (fn [log-x] (Math/exp (- log-x log-sum))) log-xs)))
   
(defn sample-from-scores [log-scores]
 (let [trg (.nextDouble +rand+)]
  (loop [so-far 0.0
         posts (indexed (log-normalize log-scores))]
   (if-let [[i p] (first posts)]
     (cond
       (< trg (+ so-far p)) i
       :default (recur (+ so-far p) (rest posts)))
     (throw (RuntimeException. "Impossible"))))))