82 lines
2.8 KiB
Text
82 lines
2.8 KiB
Text
Rebol [
|
|
title: "Rosetta code: Statistics/Normal_distribution"
|
|
file: %Statistics-Normal_distribution.r3
|
|
url: https://rosettacode.org/wiki/Statistics/Normal_distribution
|
|
author: @ldci
|
|
needs: 3.19.0
|
|
]
|
|
|
|
random/seed 1
|
|
randMax: 2147483647 ;; max integer value
|
|
nMax: 500000 ;; number of random values, can be modified
|
|
|
|
;; Normal random numbers generator using Marsaglia algorithm.
|
|
;; Generates 2 independent series of random values in the range [-1.0, 1.0].
|
|
generate: function [
|
|
n [integer!] ;; number of pairs to generate
|
|
][
|
|
m: n * 2
|
|
values: make vector! [f64! :m] ;; vector to hold generated values
|
|
for i 1 m 2 [
|
|
rsq: 0.0
|
|
;; Repeat while radius squared is outside the unit circle or zero
|
|
while [any [(rsq >= 1.0) (rsq == 0.0)]] [
|
|
x: (2.0 * random randMax) / randMax - 1.0
|
|
y: (2.0 * random randMax) / randMax - 1.0
|
|
rsq: (x * x) + (y * y)
|
|
]
|
|
f: sqrt ((-2.0 * log-e rsq) / rsq)
|
|
values/(i): x * f
|
|
values/(i + 1): y * f
|
|
]
|
|
values
|
|
]
|
|
|
|
;; Show histogram of the values vector
|
|
printHistogram: function [
|
|
values [vector!]
|
|
][
|
|
width: 50.0 ;; width of histogram bars
|
|
low: -3.0 ;; lower bound of histogram
|
|
high: 3.0 ;; upper bound of histogram
|
|
delta: 0.1 ;; bin width
|
|
n: values/length ;; length of data
|
|
nbins: to integer! ((high - low) / delta) ;; number of bins (60)
|
|
bins: make vector! [i32! :nbins] ;; initialize bins vector
|
|
repeat i n [
|
|
j: to integer! ((values/:i - low) / delta)
|
|
if all [(j >= 1) (j <= nbins)] [
|
|
bins/:j: bins/:j + 1 ;; increment bin counter
|
|
]
|
|
]
|
|
maxi: bins/maximum ;; max count in any bin
|
|
repeat j nbins [
|
|
lbin: round/to (low + j * delta - high + 0.25) 0.01 ;; low limit for bin
|
|
hbin: round/to (low + j + 1 * delta - high + 0.25) 0.01 ;; high limit for bin
|
|
s: ajoin ["[" lbin " " hbin "] "] ;; bin label string
|
|
pad s -15 ;; pad string left for alignment
|
|
k: width * bins/:j / maxi ;; number of block characters to print
|
|
while [k > 0] [
|
|
append s to-char 9609 ;; append block character (unicode 9609)
|
|
k: k - 1
|
|
]
|
|
append s ajoin [" " round/to (bins/:j * 100 / n) 0.01 "%"] ;; append percentage
|
|
print s
|
|
]
|
|
]
|
|
|
|
;;********************** Main ***********************
|
|
|
|
print "Be patient! Generating Data and Gaussian Histogram..."
|
|
print-horizontal-line
|
|
|
|
time: dt [
|
|
values: generate nMax ;; generate nMax pairs of random values
|
|
printHistogram values ;; print histogram of generated values
|
|
]
|
|
|
|
print-horizontal-line
|
|
print [nMax * 2 "Values processed in:" round/to third time 0.01 "sec"]
|
|
print ["Mean: " values/mean]
|
|
print ["STD : " values/sample-deviation]
|
|
print-horizontal-line
|