diff options
author | TenPlus1 <kinsellaja@yahoo.com> | 2015-12-08 14:47:22 +0000 |
---|---|---|
committer | TenPlus1 <kinsellaja@yahoo.com> | 2015-12-08 14:47:22 +0000 |
commit | 171a67ebf2be68ce08c706c436cb579bd2fa0ad6 (patch) | |
tree | 6db4fd6b3cdc58048f08ee303c73839c3d8136fc /statistics.lua | |
parent | dc894725c71fa2fbb6d33cd4d76c8f78a009d241 (diff) |
Code tweak and tidy, fixed rhubarb spawning
Diffstat (limited to 'statistics.lua')
-rw-r--r-- | statistics.lua | 36 |
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 |