nimbus-eth2/research/fork_choice_rule/ethresearch/distributions.nim

77 lines
2.5 KiB
Nim
Raw Normal View History

# beacon_chain
# Copyright (c) 2018 Status Research & Development GmbH
# Licensed and distributed under either of
# * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT).
# * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0).
# at your option. This file may not be copied, modified, or distributed except according to those terms.
# A port of https://github.com/ethereum/research/blob/master/clock_disparity/ghost_node.py
# Specs: https://ethresear.ch/t/beacon-chain-casper-ffg-rpj-mini-spec/2760
# Part of Casper+Sharding chain v2.1: https://notes.ethereum.org/SCIg8AH5SA-O4C1G1LYZHQ#
# Note that implementation is not updated to the latest v2.1 yet
import math, random
proc normal_distribution*(mean = 0, std = 1): int =
## Return an integer sampled from a normal distribution (gaussian)
2018-09-06 16:16:11 +00:00
## ⚠ This is not thread-safe
# Implementation via the Box-Muller method
# See https://en.wikipedia.org/wiki/BoxMuller_transform
let
mean = mean.float
std = std.float
var
z1 {.global.}: float
generate {.global.}: bool
generate = not generate
if not generate:
return int(z1 * std + mean)
let
u1 = rand(1.0)
u2 = rand(1.0)
R = sqrt(-2.0 * ln(u1))
z0 = R * cos(2 * PI * u2)
z1 = R * sin(2 * PI * u2)
return int(z0 * std + mean)
when isMainModule:
import sequtils, stats, strformat
func absolute_error(y_true, y: float): float =
## Absolute error: |y_true - y|
abs(y_true - y)
func relative_error(y_true, y: float): float =
## Relative error: |y_true - y|/|y_true|
abs(y_true - y)/abs(y_true)
let
mu = 1000
sigma = 12
a = newSeqWith(10000000, normal_distribution(mean = mu, std = sigma))
var statistics: RunningStat
for val in a:
statistics.push val
# Note: we use the sample standard deviation, not population
# See Bessel's correction and standard deviation estimation.
proc report(stat: string, value, expected: float) =
echo &"{stat:<20} {value:>9.4f} | Expected: {expected:>9.4f}"
echo &"Statistics on {a.len} samples"
report "Mean: ", statistics.mean, mu.float
report "Standard deviation: ", statistics.standardDeviationS, sigma.float
# Absolute error
doAssert absolute_error(mu.float, statistics.mean) < 0.6
doAssert absolute_error(sigma.float, statistics.standardDeviationS) < 0.01
# Relative error
doAssert relative_error(mu.float, statistics.mean) < 0.01
doAssert relative_error(sigma.float, statistics.standardDeviationS) < 0.01