Rebol3 Code Examplex
Statistics/Normal_distribution
Work with the normal distribution and its properties.
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