Skip to content

Instantly share code, notes, and snippets.

@hyiltiz
Created January 13, 2024 04:56
Show Gist options
  • Save hyiltiz/3a64f6bbe4e797890a008dfe806cd75f to your computer and use it in GitHub Desktop.
Save hyiltiz/3a64f6bbe4e797890a008dfe806cd75f to your computer and use it in GitHub Desktop.
Lazy generator for Gaussian random samples using efficient box-muller algorithm
(ns statistics.lazy.rand
(:require [clojure.math :as math]))
(defn rand-gaussian
"Get a random sample from the standard Gaussian distribution.
Optionall specify the mean m and the standard deviation sd. E.g.:
(rand-gaussian) # => 0.1324...
(rand-gaussian 5 0.1) # => 5.3397...
"
([] (rand-gaussian 0 1))
([m sd]
(defn scale [x] (+ m (* sd x)))
(def p (rand))
(def q (rand))
; We use the Box-Muller transform
(let [rho (clojure.math/sqrt (* -2 (clojure.math/log q)))
theta (* 2 clojure.math/PI p)
_box (* rho (clojure.math/cos theta))
_muller (* rho (clojure.math/sin theta))
box (scale _box)
muller (scale _muller)]
; (clojure.pprint/pprint {:mean m :sd sd
; :_box _box :_muller _muller
; :box box :muller muller})
; tail call
(lazy-seq (cons box
(cons muller
(rand-gaussian m sd)))))))
(defn sample-n
"Generate n samples based on the random sampler f. E.g.
(sample-n |(rand-int 0 3) 4) # => @[0 1 2 0]
(sample-n |(rand-uniform) 4)
"
[f n]
(take n ; (generate [_ :iterate true])
(f)))
(take 4 (rand-gaussian)) ; =>
; (-1.100044284767874
; -0.8147643792371908
; 2.09931093798589
; -0.8900201976621036)
(rand-gaussian) ; generates infinitely
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment