59 lines
1.6 KiB
Text
59 lines
1.6 KiB
Text
require "table2"
|
|
require "io2"
|
|
local fmt = require "fmt"
|
|
|
|
-- Box-Muller method from Wikipedia.
|
|
local function normal(mu, sigma)
|
|
local u1 = math.random()
|
|
local u2 = math.random()
|
|
local mag = sigma * math.sqrt(-2 * math.log(u1))
|
|
local z0 = mag * math.cos(2 * math.pi * u2) + mu
|
|
local z1 = mag * math.sin(2 * math.pi * u2) + mu
|
|
return {z0, z1}
|
|
end
|
|
|
|
local N = 100_000
|
|
local NUM_BINS = 12
|
|
|
|
local bins = table.rep(NUM_BINS, 0)
|
|
local binSize = 0.1
|
|
local samples = table.rep(N, 0)
|
|
local mu = 0.5
|
|
local sigma = 0.25
|
|
for i = 0, N / 2 - 1 do
|
|
local rns = normal(mu, sigma)
|
|
for j = 0, 1 do
|
|
local rn = rns[j + 1]
|
|
local bn
|
|
if rn < 0 then
|
|
bn = 0
|
|
elseif rn >= 1 then
|
|
bn = 11
|
|
else
|
|
bn = rn // binSize + 1
|
|
end
|
|
bins[bn + 1] += 1
|
|
samples[i * 2 + j + 1] = rn
|
|
end
|
|
end
|
|
|
|
local labels = table.rep(NUM_BINS, "")
|
|
for i = 1, NUM_BINS do
|
|
if i == 1 then
|
|
labels[i] = "-inf ..< 0.00 "
|
|
elseif i < NUM_BINS then
|
|
labels[i] = string.format("%4.2f ..< %4.2f ", binSize * (i-1), binSize * i)
|
|
else
|
|
labels[i] = "1.00 ... +inf "
|
|
end
|
|
end
|
|
|
|
fmt.print("Normal distribution with mean %0.2f and S/D %0.2f for %,s samples:\n", mu, sigma, N)
|
|
local title = " Range Number of samples within that range"
|
|
io.barchart(title, 85, labels, bins, true, true)
|
|
|
|
local m = samples:mean()
|
|
fmt.print("\nActual mean for these samples : %0.5f", m)
|
|
local c = #samples
|
|
local v = (samples:sqsum() - m * m * c) / (c - 1)
|
|
fmt.print("Actual S/D for these samples : %0.5f", math.sqrt(v))
|