diff --git a/spork/randgen.janet b/spork/randgen.janet
index 0c02130..6d0fc9a 100644
--- a/spork/randgen.janet
+++ b/spork/randgen.janet
@@ -28,6 +28,41 @@
   (def diff (- end start))
   (+ start (math/floor (* diff (rand-uniform)))))
 
+(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...
+  "
+  [&opt m sd]
+  (default m 0)
+  (default sd 1)
+  (defn scale [x] (+ m (* sd x)))
+
+  (def p (math/rng-uniform (get-rng)))
+  (def q (math/rng-uniform (get-rng)))
+
+  # We use the Box-Muller transform
+  (let [rho (math/sqrt (* -2 (math/log q)))
+        theta (* 2 math/pi p)
+        _box (* rho (math/cos theta))
+        _muller (* rho (math/sin theta))
+        box (scale _box)
+        muller (scale _muller)]
+
+    (yield box)
+    muller))
+
+(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)
+  (sample-n |(rand-gaussian 5 0.1) 4)
+  "
+  [f n]
+  (take n (generate [_ :iterate true]
+            (f))))
+
 (defn rand-index
   "Get a random numeric index of an indexed data structure"
   [xs]
@@ -85,4 +120,3 @@
   [weights & paths]
   ~(case (,rand-weights ,weights)
      ,;(array/concat @[] ;(map tuple (range (length paths)) paths))))
-
diff --git a/test/suite-randgen.janet b/test/suite-randgen.janet
new file mode 100644
index 0000000..b11e7f1
--- /dev/null
+++ b/test/suite-randgen.janet
@@ -0,0 +1,20 @@
+(use ../spork/test)
+(use ../spork/randgen)
+
+(start-suite)
+
+(assert-docs "/spork/randgen")
+
+(def delta 0.000000000000009)
+
+(assert (do
+          (set-seed 1)
+          (def expected @[5.00917440069385 
+                          5.0318756480189 
+                          4.97222830484943
+                          5.10069692026214])
+          (def actual
+            (sample-n |(rand-gaussian 5 0.1) 4))
+          (and (all |(> delta (math/abs (- $0 $1)))
+                    expected actual)))
+        "sample-rand-gaussian")