summaryrefslogtreecommitdiff
path: root/statistics.lua
diff options
context:
space:
mode:
authorTenPlus1 <kinsellaja@yahoo.com>2015-12-08 14:47:22 +0000
committerTenPlus1 <kinsellaja@yahoo.com>2015-12-08 14:47:22 +0000
commit171a67ebf2be68ce08c706c436cb579bd2fa0ad6 (patch)
tree6db4fd6b3cdc58048f08ee303c73839c3d8136fc /statistics.lua
parentdc894725c71fa2fbb6d33cd4d76c8f78a009d241 (diff)
Code tweak and tidy, fixed rhubarb spawning
Diffstat (limited to 'statistics.lua')
-rw-r--r--statistics.lua36
1 files changed, 30 insertions, 6 deletions
diff --git a/statistics.lua b/statistics.lua
index 265ccb3..c8defa2 100644
--- a/statistics.lua
+++ b/statistics.lua
@@ -1,9 +1,9 @@
local statistics = {}
-
local ROOT_2 = math.sqrt(2.0)
-- Approximations for erf(x) and erfInv(x) from
--- https://en.wikipedia.org/wiki/Error_function
+-- https://en.wikipedia.org/wiki/Error_function
+
local erf
local erf_inv
@@ -13,19 +13,26 @@ local C = 2.0/(math.pi * A)
local D = 1.0 / A
erf = function(x)
+
if x == 0 then return 0; end
+
local xSq = x * x
local aXSq = A * xSq
local v = math.sqrt(1.0 - math.exp(-xSq * (B + aXSq) / (1.0 + aXSq)))
+
return (x > 0 and v) or -v
end
erf_inv = function(x)
+
if x == 0 then return 0; end
+
if x <= -1 or x >= 1 then return nil; end
+
local y = math.log(1 - x * x)
local u = C + 0.5 * y
local v = math.sqrt(math.sqrt(u * u - D * y) - u)
+
return (x > 0 and v) or -v
end
@@ -37,6 +44,7 @@ local poisson
local cdf_table = {}
local function generate_cdf(lambda_index, lambda)
+
local max = math.ceil(4 * lambda)
local pdf = math.exp(-lambda)
local cdf = pdf
@@ -56,6 +64,7 @@ for li = 1, 100 do
end
poisson = function(lambda, max)
+
if max < 2 then
return (math.random() < math.exp(-lambda) and 0) or 1
elseif lambda >= 2 * max then
@@ -67,6 +76,7 @@ poisson = function(lambda, max)
local cdfs = cdf_table[lambda_index]
if cdfs then
+
lambda = 0.25 * lambda_index
if u < cdfs[0] then return 0; end
@@ -74,9 +84,13 @@ poisson = function(lambda, max)
if u >= cdfs[max - 1] then return max; end
if max > 4 then -- Binary search
+
local s = 0
+
while s + 1 < max do
+
local m = math.floor(0.5 * (s + max))
+
if u < cdfs[m] then max = m; else s = m; end
end
else
@@ -88,6 +102,7 @@ poisson = function(lambda, max)
return max
else
local x = lambda + math.sqrt(lambda) * std_normal(u)
+
return (x < 0.5 and 0) or (x >= max - 0.5 and max) or math.floor(x + 0.5)
end
end
@@ -102,14 +117,17 @@ statistics.erf_inv = erf_inv
--
-- @return
-- Any real number (actually between -3.0 and 3.0).
- --
+
statistics.std_normal = function()
+
local u = math.random()
+
if u < 0.001 then
return -3.0
elseif u > 0.999 then
return 3.0
end
+
return std_normal(u)
end
@@ -121,14 +139,17 @@ end
-- The distribution standard deviation.
-- @return
-- Any real number (actually between -3*sigma and 3*sigma).
- --
+
statistics.normal = function(mu, sigma)
+
local u = math.random()
+
if u < 0.001 then
return mu - 3.0 * sigma
elseif u > 0.999 then
return mu + 3.0 * sigma
end
+
return mu + sigma * std_normal(u)
end
@@ -140,11 +161,14 @@ end
-- The distribution maximum.
-- @return
-- An integer between 0 and max (both inclusive).
- --
+
statistics.poisson = function(lambda, max)
+
lambda, max = tonumber(lambda), tonumber(max)
+
if not lambda or not max or lambda <= 0 or max < 1 then return 0; end
+
return poisson(lambda, max)
end
-return statistics \ No newline at end of file
+return statistics