oh, i get it, duh. in both lua and javascript, math.random() with no arguments gives values in the range [0, 1) (non-inclusive upper bound.) so 1-math.rand() excludes zero, which is what you want. my bad.
ok @andrew, closing the loop on this dumb tangent, i found some things.
TL/DR:
i’d recommend using this “dice” method instead, which is a bit more performant and has predictable output bounds:
-- output has unit variance and is bounded to [-D, +D]
function randomGaussianDice()
local D = 3
x = 0
for i=1,D do
x = x + math.random()
end
return x*2 - D
end
- the choice of scaling in your original code isn’t arbitrary and is useful in a stats context: it gives a standard variance of 1. (that’s actually an interesting result of the approximation.) in a musical context, you are probably interested in both the bounds and the variance.
- but the downside is that the output bounds do indeed approach [+,-] infinity (i made a mistake earlier.) that also makes it potentially useful for research purposes but probably not for this context.
(here’s why:)
infs are avoided because of the non-inclusive upper bound of the random(), but you can still blow up arbitrarily large. one reason i don’t love this technique is that we rely on the random() implementation to be noninclusive and on whatever detail of that restriction to determine the output bounds range. as u approaches 0, log(u) approaches -inf (which causes the output to approach +/-inf with the other random term.) how big it actually gets is entirely implementation-defined! (max value is something like +/- e^(ε) where ε is the smallest representable floating point number. or maybe it’s 1 / 0x7ffffffff… if math.random() uses integer rand and scales it… who knows!)
-
i did some stats and benchmarks and found that with the “dice” method in lua, using only 3 dice, you get a about the same accuracy (didn’t really bother with confidence intervals but they are both within 0.1% of unit variance,) and about 20% speed boost, using math.random() to evaluate the dice. (if instead you prefer a better approximation and ~equal performance, increase dice count to 4.)
-
with the “dice” method, given normalized dice (in [0,1] inclusive or exclusive, doesn’t much matter), and taking N= number of dice, we get unit variance simply by summing the dice, multiplying by 2 and subtracting N, which also restricts our output bounds predictably to [-N, +N] (or i guess [-N, +N)?)
-
i also implemented the dice method with a custom, super simple PRNG (linear congruential.) this distribution is not statistically sound enough for, say, cryptography, but for music it’s fine and the final variance is just as good.
but somewhat weirdly, the performance was far worse in lua (in C it is far better; i guess this is about lua’s odd treatment of integer types under the hood and the fact that the library just calls down to <rand.h> in C.)
this is too bad, because 1) the dice performance boost scales with cheaper random(), 2) i think it is nice to be able to more deliberately manipulate the PRNG guts - re-seeding of course, but also being able to deliberately reduce the randomness in weird ways.
in any case, the not-very-performant lcg.lua:
lcg = {}
lcg.__index = lcg
lcg.new = function(seed, a, c)
local o = {
x = 1,
a = 1597334677,
c = 1289706101
}
if seed ~= nil then
o.x = seed
end
if a ~= nil then
o.a = a
end
if c ~= nil then
o.c = c
end
setmetatable(o, lcg)
return o
end
function lcg:seed(x) self.x = x end
function lcg:a(a) self.a = a end
function lcg:c(c) self.c = c end
function lcg:next()
self.x = (self.x * self.a + self.c) & 0x7fffffff
return self.x / 0x7fffffff
end
return lcg
i am wondering if it’s worth it to add this and some other related pseudorandom algos to norns via a C library.
if there is appetite for more random flavors/distros in the norns environment, i would be happy to add them (pink, brown, perlin, poisson, skewed gaussian, chaotic maps, etc?) and i would take some care to make them performant, accurate, usable etc. (not my first rodeo here.)
and the speed test:
socket = require 'socket' -- just for `gettime()`
function measure(name, func, n)
n = n or 10000000
print("measuring "..name.."...")
t = socket.gettime()
for i=1,n do
func()
end
elapsed = socket.gettime() - t
print("elapsed: ".. elapsed .. " seconds")
end
-- output has unit variance in (-inf, +inf)!
function randomGaussianBoxMuller()
u = 1 - math.random()
v = 1 - math.random()
return math.sqrt( -2.0 * math.log(u)) * math.cos(2.0 * math.pi * v)
end
-- output has unit variance in [-D, +D]
function randomGaussianDice()
local D = 3
x = 0
for i=1,D do
x = x + math.random()
end
return x*2 - D
end
x = 0
measure( "B-M", function() x = randomGaussianBoxMuller() end)
measure( "dice", function() x = randomGaussianDice() end)
and for completeness, the stats analysis in python, with result histograms:
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
def gauss_sqrt_log_cos(n):
u = 1 - np.random.rand(n)
v = 1 - np.random.rand(n)
x = np.sqrt(-2 * np.log(u)) * np.cos(2 * np.pi * v)
return x
# brute-force, average of uniform
def gauss_dice(n):
d = 3
x = np.random.rand(d, n)
print(x)
x = np.sum (x, axis=0)
x = x
return x* 2 - d
# for comparison
def gauss_builtin(n):
return np.random.normal(0, 1, n)
def evaluate(func, n=1000000, b=100):
name = func.__name__
x = func(n)
desc = stats.describe(x)
print("{}: {}".format(name, desc))
plt.figure()
plt.hist(x, b)
plt.title('{}, variance = {}'.format(name, desc[3]))
plt.show()
# plt.savefig('{}.png'.format(name))
print("")
evaluate(gauss_sqrt_log_cos)
#evaluate(gauss_builtin, 10000000)
evaluate(gauss_dice)
(note that the difference in the width of the bell curve is due to the fact that the dice method is predictably bounded and the direct B-M method is not.)
ok. now my brain will let me do something else.
