gauss.nim

import math, random

proc gauss(mu = 0.0, sigma = 1.0): float =
   var
      s = 0.0
      u = 0.0
      v = 0.0
   while s > 1 or s == 0:
      u = 2.0 * rand(1.0) - 1.0
      v = 2.0 * rand(1.0) - 1.0
      s = (u * u) + (v * v)
   let x = u * sqrt(-2.0 * ln(s) / s)
   result = mu + (sigma * x)

proc gauss1(mu = 0.0, sigma = 1.0): float =
   var
      x = 0.0
   while true:
      let u = 1.0 - rand(1.0)
      let v = rand(1.0)
      # Note: u > 0.0
      x = sqrt(8 / E) * (v - 0.5) / u
      if x * x <= -4.0 * ln(u): break
   result = mu + x * sigma

proc gauss2(mu = 0.0, sigma = 1.0): float =
   const
      s = 0.449871 # Constants from Leva
      t = -0.386595
      a = 0.19600
      b = 0.25472
      r1 = 0.27597
      r2 = 0.27846

   var
      u = 0.0
      v = 0.0
   while true:
      u = 1 - rand(1.0)
      v = rand(1.0) - 0.5
      # Constant 1.7156 > sqrt(8/e)
      v *= 1.7156
      # Compute Leva's quadratic form Q
      let x = u - s
      let y = abs(v) - t
      let q = x * x + y * (a * y - b * x)
      # Accept P if Q < r1 (Leva)
      # Reject P if Q > r2 (Leva)
      # Accept if v^2 <= -4 u^2 log(u) (K+M)
      if q < r1 or (q <= r2 and v * v <= -4 * u * u * ln(u)): break

   result = mu + sigma * (v / u) # Return slope

proc gaus(x, mu, sigma: float, norm = false): float =
   # Calculate a gaussian function with mean and sigma.
   # If norm=kTRUE (default is kFALSE) the result is divided
   # by sqrt(2*Pi)*sigma.

   if sigma == 0: return 1e30
   let arg = (x - mu) / sigma
   # for |arg| > 39  result is zero in double precision
   if arg < -39.0 or arg > 39.0: return
   result = exp(-0.5 * arg * arg)
   if not norm: return
   result = result / (sqrt(Tau) * sigma)

when isMainModule:
   import std / [times, stats, strformat, sums]

   #randomize()
   var rs: RunningStat
   for j in 1..5:
      for i in 1 .. 1_000:
         rs.push(gauss1())
      echo("mean: ", rs.mean,
         " stdDev: ", rs.standardDeviation())
      rs.clear()

   proc warmup() =
      # Warmup - make sure cpu is on max perf
      let start = cpuTime()
      var a = 123
      for i in 0 ..< 300_000_000:
         a += i * i mod 456
         a = a mod 789
      let dur = cpuTime() - start
      echo &"Warmup: {dur:>4.4f} s ", a

   proc printStats(name: string, stats: RunningStat, dur: float) =
      echo &"""{name}:
      Collected {stats.n} samples in {dur:>4.4f} s
      Average time: {stats.mean * 1000:>4.4f} ms
      Stddev  time: {stats.standardDeviationS * 1000:>4.4f} ms
      Min     time: {stats.min * 1000:>4.4f} ms
      Max     time: {stats.max * 1000:>4.4f} ms"""

   template bench(name, samples, code: untyped) =
      proc runBench() {.gensym, nimcall.} =
         var stats: RunningStat
         let globalStart = cpuTime()
         for i {.inject.} in 0 ..< samples:
            let start = cpuTime()
            code
            let duration = cpuTime() - start
            stats.push duration
         let globalDuration = cpuTime() - globalStart
         printStats(name, stats, globalDuration)
      runBench()

   proc main =
      const n = 10_000
      warmup()

      var s = newSeq[float](n)
      bench("gauss polar", n):
         s[i] = gauss()
      echo s.sumKbn

      bench("gauss py", n):
         s[i] = gauss1()
      echo s.sumKbn

      bench("gauss gsl", n):
         s[i] = gauss2()
      echo s.sumKbn

   main()

mean: -0.01044065606097552 stdDev: 1.001754054755673
mean: -0.01984169728775705 stdDev: 1.007737011761979
mean: 0.01861571323384367 stdDev: 1.012283632285083
mean: 0.04522488626175587 stdDev: 1.008427782476539
mean: 0.001956187534907565 stdDev: 1.05766215589672
Warmup: 0.9025 s 224
gauss polar:
      Collected 10000 samples in 0.0070 s
      Average time: 0.0004 ms
      Stddev  time: 0.0000 ms
      Min     time: 0.0002 ms
      Max     time: 0.0012 ms
134.6366307323438
gauss py:
      Collected 10000 samples in 0.0069 s
      Average time: 0.0004 ms
      Stddev  time: 0.0000 ms
      Min     time: 0.0002 ms
      Max     time: 0.0023 ms
-37.30513539942112
gauss gsl:
      Collected 10000 samples in 0.0068 s
      Average time: 0.0003 ms
      Stddev  time: 0.0000 ms
      Min     time: 0.0000 ms
      Max     time: 0.0012 ms
-59.05391651987865

With 1M samples

mean: -0.01044065606097552 stdDev: 1.001754054755673
mean: -0.01984169728775705 stdDev: 1.007737011761979
mean: 0.01861571323384367 stdDev: 1.012283632285083
mean: 0.04522488626175587 stdDev: 1.008427782476539
mean: 0.001956187534907565 stdDev: 1.05766215589672
Warmup: 0.9025 s 224
gauss polar:
      Collected 1000000 samples in 0.6946 s
      Average time: 0.0004 ms
      Stddev  time: 0.0000 ms
      Min     time: 0.0000 ms
      Max     time: 0.0019 ms
-1816.760808309902
gauss py:
      Collected 1000000 samples in 0.6837 s
      Average time: 0.0003 ms
      Stddev  time: 0.0000 ms
      Min     time: 0.0000 ms
      Max     time: 0.0018 ms
1744.328667953836
gauss gsl:
      Collected 1000000 samples in 0.6790 s
      Average time: 0.0003 ms
      Stddev  time: 0.0000 ms
      Min     time: 0.0000 ms
      Max     time: 0.0021 ms
-2288.719327933872