Statistics/Normal distribution
You are encouraged to solve this task according to the task description, using any language you may know.
The Normal (or Gaussian) distribution is a frequently used distribution in statistics. While most programming languages provide a uniformly distributed random number generator, one can derive normally distributed random numbers from a uniform generator.
- The task
- Take a uniform random number generator and create a large (you decide how large) set of numbers that follow a normal (Gaussian) distribution. Calculate the dataset's mean and stddev, and show the histogram here.
- Mention any native language support for the generation of normally distributed random numbers.
- Reference
- You may refer to code in Statistics/Basic if available.
C#
<lang csharp>using System; using MathNet.Numerics.Distributions; using MathNet.Numerics.Statistics;
class Program {
static void RunNormal(int sampleSize)
{
double[] X = new double[sampleSize];
var norm = new Normal(new Random());
norm.Samples(X);
const int numBuckets = 10;
var histogram = new Histogram(X, numBuckets);
Console.WriteLine("Sample size: {0:N0}", sampleSize);
for (int i = 0; i < numBuckets; i++)
{
string bar = new String('#', (int)(histogram[i].Count * 360 / sampleSize));
Console.WriteLine(" {0:0.00} : {1}", histogram[i].LowerBound, bar);
}
var statistics = new DescriptiveStatistics(X);
Console.WriteLine(" Mean: " + statistics.Mean);
Console.WriteLine("StdDev: " + statistics.StandardDeviation);
Console.WriteLine();
}
static void Main(string[] args)
{
RunNormal(100);
RunNormal(1000);
RunNormal(10000);
}
}</lang>
- Output:
Sample size: 100 -2.12 : ####### -1.66 : ############################ -1.19 : ####################################### -0.72 : ############################################## -0.26 : ############################################################################### 0.21 : ###################################################################################### 0.68 : ################################ 1.14 : ######################### 1.61 : ### 2.07 : ########## Mean: 0.0394411345297757 StdDev: 0.925286665513647 Sample size: 1,000 -2.98 : ## -2.34 : ########## -1.69 : ############################## -1.05 : ################################################################ -0.40 : ########################################################################################### 0.24 : ######################################################################################## 0.88 : ############################################## 1.53 : ################## 2.17 : ##### 2.82 : ## Mean: 0.0868718238400114 StdDev: 0.989120264661867 Sample size: 10,000 -3.88 : -3.12 : ## -2.35 : ################# -1.59 : #################################################### -0.82 : ################################################################################################ -0.06 : #################################################################################################### 0.71 : ############################################################### 1.47 : ##################### 2.23 : #### 3.00 : Mean: 0.0208920122989818 StdDev: 1.00046328880424
C++
showing features of C++11 here <lang cpp>#include <random>
- include <map>
- include <string>
- include <iostream>
- include <cmath>
- include <iomanip>
int main( ) {
std::random_device myseed ;
std::mt19937 engine ( myseed( ) ) ;
std::normal_distribution<> normDistri ( 2 , 3 ) ;
std::map<int , int> normalFreq ;
int sum = 0 ; //holds the sum of the randomly created numbers
double mean = 0.0 ;
double stddev = 0.0 ;
for ( int i = 1 ; i < 10001 ; i++ )
++normalFreq[ normDistri ( engine ) ] ;
for ( auto MapIt : normalFreq ) {
sum += MapIt.first * MapIt.second ;
}
mean = sum / 10000 ;
stddev = sqrt( sum / 10000 ) ;
std::cout << "The mean of the distribution is " << mean << " , the " ;
std::cout << "standard deviation " << stddev << " !\n" ;
std::cout << "And now the histogram:\n" ;
for ( auto MapIt : normalFreq ) {
std::cout << std::left << std::setw( 4 ) << MapIt.first <<
std::string( MapIt.second / 100 , '*' ) << std::endl ;
} return 0 ;
}</lang> Output:
The mean of the distribution is 1 , the standard deviation 1 ! And now the histogram: -10 -9 -8 -7 -6 -5 -4 * -3 ** -2 **** -1 ****** 0 ********************* 1 ************ 2 ************ 3 *********** 4 ********* 5 ****** 6 **** 7 ** 8 * 9 10 11 12 13
D
This uses the Box-Muller method as in the Go entry, and the module from the Statistics/Basic. A ziggurat-based normal generator for the Phobos standard library is in the works. <lang d>import std.stdio, std.random, std.math, std.range, std.algorithm,
statistics_basic;
struct Normals {
double mu, sigma; double[2] state; size_t index = state.length; enum empty = false;
void popFront() pure nothrow { index++; }
@property double front() {
if (index >= state.length) {
immutable r = sqrt(-2 * uniform!"]["(0., 1.0).log) * sigma;
immutable x = 2 * PI * uniform01;
state = [mu + r * x.sin, mu + r * x.cos];
index = 0;
}
return state[index];
}
}
void main() {
const data = Normals(0.0, 0.5).take(100_000).array;
writefln("Mean: %8.6f, SD: %8.6f\n", data.meanStdDev[]);
data.map!q{ max(0.0, min(0.9999, a / 3 + 0.5)) }.showHistogram01;
}</lang>
- Output:
Mean: 0.000528, SD: 0.502245 0.0: * 0.1: ****** 0.2: ***************** 0.3: *********************************** 0.4: ************************************************* 0.5: ************************************************** 0.6: ********************************** 0.7: ***************** 0.8: ****** 0.9: *
Elixir
<lang elixir>defmodule Statistics do
def normal_distribution(n, w\\5) do
{sum, sum2, hist} = generate(n, w)
mean = sum / n
stddev = :math.sqrt(sum2 / n - mean*mean)
IO.puts "size: #{n}"
IO.puts "mean: #{mean}"
IO.puts "stddev: #{stddev}"
{min, max} = Map.to_list(hist)
|> Enum.filter_map(fn {_k,v} -> v >= n/120/w end, fn {k,_v} -> k end)
|> Enum.min_max
Enum.each(min..max, fn i ->
bar = String.duplicate("=", trunc(120 * w * Map.get(hist, i, 0) / n))
:io.fwrite "~4.1f: ~s~n", [i/w, bar]
end)
IO.puts ""
end
defp generate(n, w) do
Enum.reduce(1..n, {0, 0, %{}}, fn _,{sum, sum2, hist} ->
z = :rand.normal
{sum+z, sum2+z*z, Map.update(hist, round(w*z), 1, &(&1+1))}
end)
end
end
Enum.each([100,1000,10000], fn n ->
Statistics.normal_distribution(n)
end)</lang>
- Output:
size: 100 mean: 0.027742416007234007 stddev: 1.0209597927405403 -2.6: ============ -2.4: -2.2: ============ -2.0: ====== -1.8: -1.6: -1.4: ============================== -1.2: ====== -1.0: ============================== -0.8: ========================================== -0.6: ========================================== -0.4: ================================================ -0.2: ================================================ 0.0: ============================== 0.2: ==================================== 0.4: ========================================== 0.6: ====================================================== 0.8: ========================================== 1.0: ================================================ 1.2: ============================== 1.4: ====== 1.6: ============ 1.8: ============ 2.0: 2.2: 2.4: ====== 2.6: ====== size: 1000 mean: -0.025562168667763084 stddev: 1.0338288521306742 -3.2: = -3.0: -2.8: = -2.6: === -2.4: == -2.2: ====== -2.0: == -1.8: ============= -1.6: =============== -1.4: ================= -1.2: ================= -1.0: ==================================== -0.8: =================================== -0.6: ============================================ -0.4: ============================================ -0.2: =============================================== 0.0: ========================================= 0.2: =========================================== 0.4: ============================================= 0.6: ======================================= 0.8: ================================ 1.0: ============================ 1.2: ======================== 1.4: ================== 1.6: ========== 1.8: ===== 2.0: ======== 2.2: ==== 2.4: ===== 2.6: = 2.8: = size: 10000 mean: -0.009132420943007152 stddev: 0.9979508347451509 -2.6: = -2.4: === -2.2: ==== -2.0: ===== -1.8: ========= -1.6: ============== -1.4: ================ -1.2: ======================= -1.0: ============================ -0.8: ================================= -0.6: ============================================ -0.4: =========================================== -0.2: ============================================== 0.0: ================================================== 0.2: ============================================ 0.4: =========================================== 0.6: ======================================= 0.8: ===================================== 1.0: ============================ 1.2: ======================= 1.4: ================ 1.6: ============== 1.8: ========= 2.0: ====== 2.2: === 2.4: == 2.6: =
Fortran
Using the Marsaglia polar method <lang fortran>program Normal_Distribution
implicit none
integer, parameter :: i64 = selected_int_kind(18)
integer, parameter :: r64 = selected_real_kind(15)
integer(i64), parameter :: samples = 1000000_i64
real(r64) :: mean, stddev
real(r64) :: sumn = 0, sumnsq = 0
integer(i64) :: n = 0
integer(i64) :: bin(-50:50) = 0
integer :: i, ind
real(r64) :: ur1, ur2, nr1, nr2, s
n = 0
do while(n <= samples)
call random_number(ur1)
call random_number(ur2)
ur1 = ur1 * 2.0 - 1.0
ur2 = ur2 * 2.0 - 1.0
s = ur1*ur1 + ur2*ur2
if(s >= 1.0_r64) cycle
nr1 = ur1 * sqrt(-2.0*log(s)/s)
ind = floor(5.0*nr1)
bin(ind) = bin(ind) + 1_i64
sumn = sumn + nr1
sumnsq = sumnsq + nr1*nr1
nr2 = ur2 * sqrt(-2.0*log(s)/s)
ind = floor(5.0*nr2)
bin(ind) = bin(ind) + 1_i64
sumn = sumn + nr2
sumnsq = sumnsq + nr2*nr2
n = n + 2_i64
end do
mean = sumn / n
stddev = sqrt(sumnsq/n - mean*mean)
write(*, "(a, i0)") "sample size = ", samples
write(*, "(a, f17.15)") "Mean : ", mean,
write(*, "(a, f17.15)") "Stddev : ", stddev
do i = -15, 15
write(*, "(f4.1, a, a)") real(i)/5.0, ": ", repeat("=", int(bin(i)*500/samples))
end do
end program</lang>
- Output:
sample size = 1000 Mean : 0.043096320705032 Stddev : 0.981532585231540 -3.0: -2.8: -2.6: == -2.4: == -2.2: ==== -2.0: ====== -1.8: ======= -1.6: ============ -1.4: ================ -1.2: ===================== -1.0: =========================== -0.8: ======================= -0.6: ================================== -0.4: ===================================== -0.2: ========================================== 0.0: =============================================== 0.2: ==================================== 0.4: ================================= 0.6: ================================== 0.8: ============================= 1.0: ==================== 1.2: ========================== 1.4: =========== 1.6: ========= 1.8: ==== 2.0: ====== 2.2: === 2.4: 2.6: 2.8: = 3.0: sample size = 1000000 Mean : 0.000166653231289 Stddev : 1.000025612171690 -3.0: -2.8: = -2.6: = -2.4: == -2.2: ==== -2.0: ====== -1.8: ========= -1.6: ============ -1.4: ================= -1.2: ===================== -1.0: ========================== -0.8: =============================== -0.6: =================================== -0.4: ====================================== -0.2: ======================================= 0.0: ======================================= 0.2: ====================================== 0.4: ================================== 0.6: =============================== 0.8: ========================== 1.0: ===================== 1.2: ================= 1.4: ============ 1.6: ========= 1.8: ====== 2.0: ==== 2.2: == 2.4: = 2.6: = 2.8: 3.0:
FreeBASIC
<lang freebasic>' FB 1.05.0 Win64
Const pi As Double = 3.141592653589793 Randomize
' Generates normally distributed random numbers with mean 0 and standard deviation 1 Function randomNormal() As Double
Return Cos(2.0 * pi * Rnd) * Sqr(-2.0 * Log(Rnd))
End Function
Sub normalStats(sampleSize As Integer)
If sampleSize < 1 Then Return Dim r(1 To sampleSize) As Double Dim h(-1 To 10) As Integer all zero by default Dim sum As Double = 0.0 Dim hSum As Integer = 0
' Generate 'sampleSize' normally distributed random numbers with mean 0.5 and standard deviation 0.25
' calculate their sum
' and in which box they will fall when drawing the histogram
For i As Integer = 1 To sampleSize
r(i) = 0.5 + randomNormal / 4.0
sum += r(i)
If r(i) < 0.0 Then
h(-1) += 1
ElseIf r(i) >= 1.0 Then
h(10) += 1
Else
h(Int(r(i) * 10)) += 1
End If
Next
For i As Integer = -1 To 10 : hSum += h(i) : Next
' adjust one of the h() values if necessary to ensure hSum = sampleSize
Dim adj As Integer = sampleSize - hSum
If adj <> 0 Then
For i As Integer = -1 To 10
h(i) += adj
If h(i) >= 0 Then Exit For
h(i) -= adj
Next
End If
Dim mean As Double = sum / sampleSize
Dim sd As Double sum = 0.0 ' Now calculate their standard deviation For i As Integer = 1 To sampleSize sum += (r(i) - mean) ^ 2.0 Next sd = Sqr(sum/sampleSize)
' Draw a histogram of the data with interval 0.1
Dim numStars As Integer
' If sample size > 300 then normalize histogram to 300
Dim scale As Double = 1.0
If sampleSize > 300 Then scale = 300.0 / sampleSize
Print "Sample size "; sampleSize
Print
Print Using " Mean #.######"; mean;
Print Using " SD #.######"; sd
Print
For i As Integer = -1 To 10
If i = -1 Then
Print Using "< 0.00 : ";
ElseIf i = 10 Then
Print Using ">=1.00 : ";
Else
Print Using " #.## : "; i/10.0;
End If
Print Using "##### " ; h(i);
numStars = Int(h(i) * scale + 0.5)
Print String(numStars, "*")
Next
End Sub
normalStats 100 Print normalStats 1000 Print normalStats 10000 Print normalStats 100000 Print Print "Press any key to quit" Sleep</lang> Sample output:
- Output:
Sample size 100 Mean 0.486977 SD 0.244147 < 0.00 : 2 ** 0.00 : 5 ***** 0.10 : 4 **** 0.20 : 14 ************** 0.30 : 12 ************ 0.40 : 15 *************** 0.50 : 17 ***************** 0.60 : 11 *********** 0.70 : 9 ********* 0.80 : 7 ******* 0.90 : 1 * >=1.00 : 3 *** Sample size 1000 Mean 0.489234 SD 0.247606 < 0.00 : 18 ***** 0.00 : 32 ********** 0.10 : 73 ********************** 0.20 : 111 ********************************* 0.30 : 138 ***************************************** 0.40 : 151 ********************************************* 0.50 : 153 ********************************************** 0.60 : 114 ********************************** 0.70 : 101 ****************************** 0.80 : 51 *************** 0.90 : 38 *********** >=1.00 : 20 ****** Sample size 10000 Mean 0.498239 SD 0.249235 < 0.00 : 225 ******* 0.00 : 333 ********** 0.10 : 589 ****************** 0.20 : 999 ****************************** 0.30 : 1320 **************************************** 0.40 : 1542 ********************************************** 0.50 : 1581 *********************************************** 0.60 : 1323 **************************************** 0.70 : 963 ***************************** 0.80 : 591 ****************** 0.90 : 314 ********* >=1.00 : 220 ******* Sample size 100000 Mean 0.500925 SD 0.248910 < 0.00 : 2173 ******* 0.00 : 3192 ********** 0.10 : 5938 ****************** 0.20 : 9715 ***************************** 0.30 : 13351 **************************************** 0.40 : 15399 ********************************************** 0.50 : 15680 *********************************************** 0.60 : 13422 **************************************** 0.70 : 9633 ***************************** 0.80 : 5993 ****************** 0.90 : 3207 ********** >=1.00 : 2297 *******
Go
Box-Muller method shown here. Go has a normally distributed random function in the standard library, as shown in the Go Random numbers solution. It uses the ziggurat method. <lang go>package main
import (
"fmt" "math" "math/rand" "strings"
)
// Box-Muller func norm2() (s, c float64) {
r := math.Sqrt(-2 * math.Log(rand.Float64())) s, c = math.Sincos(2 * math.Pi * rand.Float64()) return s * r, c * r
}
func main() {
const (
n = 10000
bins = 12
sig = 3
scale = 100
)
var sum, sumSq float64
h := make([]int, bins)
for i, accum := 0, func(v float64) {
sum += v
sumSq += v * v
b := int((v + sig) * bins / sig / 2)
if b >= 0 && b < bins {
h[b]++
}
}; i < n/2; i++ {
v1, v2 := norm2()
accum(v1)
accum(v2)
}
m := sum / n
fmt.Println("mean:", m)
fmt.Println("stddev:", math.Sqrt(sumSq/float64(n)-m*m))
for _, p := range h {
fmt.Println(strings.Repeat("*", p/scale))
}
}</lang> Output:
mean: -0.0034970888831523488 stddev: 1.0040682925006286 * **** ********* *************** ******************* ****************** ************** ********* **** *
Haskell
<lang haskell>import Data.Map (Map, empty, insert, findWithDefault, toList) import Data.Maybe (fromMaybe) import Text.Printf (printf) import Data.Function (on) import Data.List (sort, maximumBy, minimumBy) import Control.Monad.Random (RandomGen, Rand, evalRandIO, getRandomR) import Control.Monad (replicateM)
-- Box-Muller getNorm :: RandomGen g => Rand g Double getNorm = do
u0 <- getRandomR (0.0, 1.0)
u1 <- getRandomR (0.0, 1.0)
let r = sqrt $ (-2.0) * log u0
theta = 2.0 * pi * u1
return $ r * sin theta
putInBin :: Double -> Map Int Int -> Double -> Map Int Int putInBin binWidth t v =
let bin = round (v / binWidth)
count = findWithDefault 0 bin t
in insert bin (count+1) t
runTest :: Int -> IO () runTest n = do
rs <- evalRandIO $ replicateM n getNorm let binWidth = 0.1
tally v (sv, sv2, t) = (sv+v, sv2 + v*v, putInBin binWidth t v)
(sum, sum2, tallies) = foldr tally (0.0, 0.0, empty) rs
tallyList = sort $ toList tallies
printStars tallies binWidth maxCount selection =
let count = findWithDefault 0 selection tallies
bin = binWidth * fromIntegral selection
maxStars = 100
starCount = if maxCount <= maxStars
then count
else maxStars * count `div` maxCount
stars = replicate starCount '*'
in printf "%5.2f: %s %d\n" bin stars count
mean = sum / fromIntegral n
stddev = sqrt (sum2/fromIntegral n - mean*mean)
printf "\n" printf "sample count: %d\n" n printf "mean: %9.7f\n" mean printf "stddev: %9.7f\n" stddev
let maxCount = snd $ maximumBy (compare `on` snd) tallyList
maxBin = fst $ maximumBy (compare `on` fst) tallyList
minBin = fst $ minimumBy (compare `on` fst) tallyList
mapM_ (printStars tallies binWidth maxCount) [minBin..maxBin]
main = do
runTest 1000 runTest 2000000</lang>
- Output:
sample count: 1000 mean: -0.0269949 stddev: 0.9795285 -3.10: ** 2 -3.00: 0 -2.90: 0 -2.80: ** 2 -2.70: * 1 -2.60: **** 4 -2.50: ** 2 -2.40: ** 2 -2.30: 0 -2.20: *** 3 -2.10: ***** 5 -2.00: ****** 6 -1.90: ****** 6 -1.80: *********** 11 -1.70: ************ 12 -1.60: ******* 7 -1.50: ************* 13 -1.40: ***************** 17 -1.30: ******************** 20 -1.20: **************** 16 -1.10: ***************** 17 -1.00: ********************** 22 -0.90: *************************** 27 -0.80: ********************** 22 -0.70: ******************************** 32 -0.60: ********************************* 33 -0.50: ****************************************** 42 -0.40: *********************************************** 47 -0.30: ************************************************ 48 -0.20: *************************** 27 -0.10: ***************************** 29 0.00: *********************************************** 47 0.10: *************************************************** 51 0.20: ****************************************** 42 0.30: ******************************** 32 0.40: ********************************* 33 0.50: ***************************************** 41 0.60: ****************************************** 42 0.70: **************************** 28 0.80: ********************** 22 0.90: *************************** 27 1.00: ******************* 19 1.10: ********************** 22 1.20: ************************ 24 1.30: ******************** 20 1.40: ***************** 17 1.50: ********** 10 1.60: ************* 13 1.70: **** 4 1.80: *** 3 1.90: ******* 7 2.00: ****** 6 2.10: * 1 2.20: * 1 2.30: ******* 7 2.40: *** 3 2.50: 0 2.60: * 1 2.70: 0 2.80: 0 2.90: 0 3.00: * 1 3.10: 0 3.20: 0 3.30: * 1 sample count: 2000000 mean: 0.0001017 stddev: 0.9994329 -4.60: 3 -4.50: 2 -4.40: 3 -4.30: 9 -4.20: 18 -4.10: 19 -4.00: 20 -3.90: 41 -3.80: 77 -3.70: 84 -3.60: 105 -3.50: 189 -3.40: 245 -3.30: 350 -3.20: 460 -3.10: 619 -3.00: * 838 -2.90: * 1234 -2.80: * 1586 -2.70: ** 2063 -2.60: *** 2716 -2.50: **** 3503 -2.40: ***** 4345 -2.30: ******* 5678 -2.20: ******** 7160 -2.10: *********** 8856 -2.00: ************* 10915 -1.90: **************** 13299 -1.80: ******************* 15599 -1.70: *********************** 19004 -1.60: *************************** 22321 -1.50: ******************************** 25940 -1.40: ************************************* 29622 -1.30: ****************************************** 34213 -1.20: ************************************************ 38922 -1.10: ****************************************************** 43415 -1.00: ************************************************************ 48250 -0.90: ****************************************************************** 53210 -0.80: ************************************************************************ 58127 -0.70: ****************************************************************************** 62438 -0.60: *********************************************************************************** 66650 -0.50: **************************************************************************************** 70298 -0.40: ******************************************************************************************** 73739 -0.30: *********************************************************************************************** 75831 -0.20: ************************************************************************************************** 78222 -0.10: *************************************************************************************************** 79412 0.00: **************************************************************************************************** 79801 0.10: *************************************************************************************************** 79255 0.20: ************************************************************************************************* 78163 0.30: ************************************************************************************************ 76667 0.40: ******************************************************************************************** 73554 0.50: **************************************************************************************** 70391 0.60: *********************************************************************************** 66566 0.70: ****************************************************************************** 62857 0.80: ************************************************************************ 57962 0.90: ****************************************************************** 53407 1.00: ************************************************************ 48565 1.10: ****************************************************** 43496 1.20: ************************************************ 38799 1.30: ****************************************** 34156 1.40: ************************************* 29713 1.50: ******************************** 25946 1.60: *************************** 22264 1.70: *********************** 18843 1.80: ******************* 15780 1.90: **************** 13151 2.00: ************* 10905 2.10: ********** 8690 2.20: ******** 7102 2.30: ******* 5693 2.40: ***** 4459 2.50: **** 3550 2.60: *** 2603 2.70: ** 2155 2.80: ** 1619 2.90: * 1121 3.00: * 914 3.10: 607 3.20: 478 3.30: 349 3.40: 216 3.50: 170 3.60: 113 3.70: 79 3.80: 58 3.90: 48 4.00: 33 4.10: 20 4.20: 9 4.30: 8 4.40: 7 4.50: 3 4.60: 3 4.70: 0 4.80: 1 4.90: 1
J
Solution <lang j>runif01=: ?@$ 0: NB. random uniform number generator rnorm01=. (2 o. 2p1 * runif01) * [: %: _2 * ^.@runif01 NB. random normal number generator (Box-Muller)
mean=: +/ % # NB. mean stddev=: (<:@# %~ +/)&.:*:@(- mean) NB. standard deviation histogram=: <:@(#/.~)@(i.@#@[ , I.)</lang> Example Usage <lang j> DataSet=: rnorm01 1e5
(mean , stddev) DataSet
0.000781667 1.00154
require 'plot' plot (5 %~ i: 25) ([;histogram) DataSet</lang>
Java
<lang java>import static java.lang.Math.*; import static java.util.Arrays.stream; import java.util.Locale; import java.util.function.DoubleSupplier; import static java.util.stream.Collectors.joining; import java.util.stream.DoubleStream; import static java.util.stream.IntStream.range;
public class Test implements DoubleSupplier {
private double mu, sigma; private double[] state = new double[2]; private int index = state.length;
Test(double m, double s) {
mu = m;
sigma = s;
}
static double[] meanStdDev(double[] numbers) {
if (numbers.length == 0)
return new double[]{0.0, 0.0};
double sx = 0.0, sxx = 0.0;
long n = 0;
for (double x : numbers) {
sx += x;
sxx += pow(x, 2);
n++;
}
return new double[]{sx / n, pow((n * sxx - pow(sx, 2)), 0.5) / n};
}
static String replicate(int n, String s) {
return range(0, n + 1).mapToObj(i -> s).collect(joining());
}
static void showHistogram01(double[] numbers) {
final int maxWidth = 50;
long[] bins = new long[10];
for (double x : numbers)
bins[(int) (x * bins.length)]++;
double maxFreq = stream(bins).max().getAsLong();
for (int i = 0; i < bins.length; i++)
System.out.printf(" %3.1f: %s%n", i / (double) bins.length,
replicate((int) (bins[i] / maxFreq * maxWidth), "*"));
System.out.println();
}
@Override
public double getAsDouble() {
index++;
if (index >= state.length) {
double r = sqrt(-2 * log(random())) * sigma;
double x = 2 * PI * random();
state = new double[]{mu + r * sin(x), mu + r * cos(x)};
index = 0;
}
return state[index];
}
public static void main(String[] args) {
Locale.setDefault(Locale.US);
double[] data = DoubleStream.generate(new Test(0.0, 0.5)).limit(100_000)
.toArray();
double[] res = meanStdDev(data);
System.out.printf("Mean: %8.6f, SD: %8.6f%n", res[0], res[1]);
showHistogram01(stream(data).map(a -> max(0.0, min(0.9999, a / 3 + 0.5)))
.toArray());
}
}</lang>
Mean: -0.001870, SD: 0.500539 0.0: ** 0.1: ******* 0.2: ****************** 0.3: ************************************ 0.4: *************************************************** 0.5: ************************************************** 0.6: *********************************** 0.7: ****************** 0.8: ******* 0.9: **
Lasso
<lang Lasso>define stat1(a) => { if(#a->size) => { local(mean = (with n in #a sum #n) / #a->size) local(sdev = math_pow(((with n in #a sum Math_Pow((#n - #mean),2)) / #a->size),0.5)) return (:#sdev, #mean) else return (:0,0) } } define stat2(a) => { if(#a->size) => { local(sx = 0, sxx = 0) with x in #a do => { #sx += #x #sxx += #x*#x } local(sdev = math_pow((#a->size * #sxx - #sx * #sx),0.5) / #a->size) return (:#sdev, #sx / #a->size) else return (:0,0) } } define histogram(a) => { local( out = '\r', h = array(0,0,0,0,0,0,0,0,0,0,0), maxwidth = 50, sc = 0 ) with n in #a do => { if((#n * 10) <= 0) => { #h->get(1) += 1 else((#n * 10) >= 10) #h->get(#h->size) += 1 else #h->get(integer(decimal(#n)*10)+1) += 1 }
} local(mx = decimal(with n in #h max #n)) with i in #h do => { #out->append((#sc/10.0)->asString(-precision=1)+': '+('+' * integer(#i / #mx * #maxwidth))+'\r') #sc++ } return #out } define normalDist(mean,sdev) => { // Uses Box-Muller transform return ((-2 * decimal_random->log)->sqrt * (2 * pi * decimal_random)->cos) * #sdev + #mean }
with scale in array(100,1000,10000) do => {^ local(n = array) loop(#scale) => { #n->insert(normalDist(0.5, 0.2)) } local(sdev1,mean1) = stat1(#n) local(sdev2,mean2) = stat2(#n) #scale' numbers:\r'
'Naive method: sd: '+#sdev1+', mean: '+#mean1+'\r' 'Second method: sd: '+#sdev2+', mean: '+#mean2+'\r' histogram(#n) '\r\r'
^}</lang>
- Output:
100 numbers: Naive method: sd: 0.199518, mean: 0.506059 Second method: sd: 0.199518, mean: 0.506059 0.0: ++ 0.1: ++++ 0.2: +++++++++++++++++ 0.3: ++++++++++++++++++++++ 0.4: ++++++++++++++++++++++++++++++++++++++++++++++++++ 0.5: +++++++++++++++++++++++++++++++++++++++ 0.6: +++++++++++++++++++++++++++++++++ 0.7: ++++++++++++++++++++++++ 0.8: ++++++++++++++++++++ 0.9: ++++ 1.0: ++ 1000 numbers: Naive method: sd: 0.199653, mean: 0.504046 Second method: sd: 0.199653, mean: 0.504046 0.0: +++ 0.1: ++++++ 0.2: ++++++++++++++++ 0.3: ++++++++++++++++++++++++++++++ 0.4: +++++++++++++++++++++++++++++++++++++++++++++++ 0.5: ++++++++++++++++++++++++++++++++++++++++++++++++++ 0.6: ++++++++++++++++++++++++++++++++++++++++++++++ 0.7: +++++++++++++++++++++++++ 0.8: +++++++++++++++++++ 0.9: +++++++ 1.0: ++++ 10000 numbers: Naive method: sd: 0.202354, mean: 0.502519 Second method: sd: 0.202354, mean: 0.502519 0.0: +++ 0.1: +++++++ 0.2: +++++++++++++++ 0.3: +++++++++++++++++++++++++++++ 0.4: ++++++++++++++++++++++++++++++++++++++++++ 0.5: ++++++++++++++++++++++++++++++++++++++++++++++++++ 0.6: +++++++++++++++++++++++++++++++++++++++++++ 0.7: ++++++++++++++++++++++++++++++ 0.8: ++++++++++++++++ 0.9: +++++++ 1.0: ++++
Liberty BASIC
Uses LB Statistics/Basic <lang lb>call sample 100000
end
sub sample n
dim dat( n)
for i =1 to n
dat( i) =normalDist( 1, 0.2)
next i
'// show mean, standard deviation. Find max, min.
mx =-1000
mn = 1000
sum =0
sSq =0
for i =1 to n
d =dat( i)
mx =max( mx, d)
mn =min( mn, d)
sum =sum +d
sSq =sSq +d^2
next i
print n; " data terms used."
mean =sum / n print "Largest term was "; mx; " & smallest was "; mn range =mx -mn print "Mean ="; mean
print "Stddev ="; ( sSq /n -mean^2)^0.5
'// show histogram
nBins =50
dim bins( nBins)
for i =1 to n
z =int( ( dat( i) -mn) /range *nBins)
bins( z) =bins( z) +1
next i
for b =0 to nBins -1
for j =1 to int( nBins *bins( b)) /n *30)
print "#";
next j
print
next b
print
end sub
function normalDist( m, s) ' Box Muller method
u =rnd( 1) v =rnd( 1) normalDist =( -2 *log( u))^0.5 *cos( 2 *3.14159265 *v)
end function</lang>
100000 data terms used. Largest term was 4.12950792 & smallest was -4.37934139 Mean =-0.26785425e-2 Stddev =1.00097669
# ## ### ##### ######## ############ ################ ######################## ############################## ###################################### ############################################## ######################################################## ################################################################### ############################################################################## ####################################################################################### ################################################################################################ #################################################################################################### ######################################################################################################## ##################################################################################################### ############################################################################################## ######################################################################################### ################################################################################## ######################################################################### ############################################################## #################################################### ########################################## ################################# ########################## ################## ############# ######### ###### #### ## # #
Lua
Lua provides math.random() to generate uniformly distributed random numbers. The function gaussian() shown here uses math.random() to generate normally distributed random numbers with given mean and variance. <lang Lua>function gaussian (mean, variance)
return math.sqrt(-2 * variance * math.log(math.random())) *
math.cos(2 * variance * math.pi * math.random()) + mean
end
function mean (t)
local sum = 0
for k, v in pairs(t) do
sum = sum + v
end
return sum / #t
end
function std (t)
local squares, avg = 0, mean(t)
for k, v in pairs(t) do
squares = squares + ((avg - v) ^ 2)
end
local variance = squares / #t
return math.sqrt(variance)
end
function showHistogram (t)
local lo = math.ceil(math.min(unpack(t)))
local hi = math.floor(math.max(unpack(t)))
local hist, barScale = {}, 200
for i = lo, hi do
hist[i] = 0
for k, v in pairs(t) do
if math.ceil(v - 0.5) == i then
hist[i] = hist[i] + 1
end
end
io.write(i .. "\t" .. string.rep('=', hist[i] / #t * barScale))
print(" " .. hist[i])
end
end
math.randomseed(os.time()) local t, average, variance = {}, 50, 10 for i = 1, 1000 do
table.insert(t, gaussian(average, variance))
end print("Mean:", mean(t) .. ", expected " .. average) print("StdDev:", std(t) .. ", expected " .. math.sqrt(variance) .. "\n") showHistogram(t)</lang>
- Output:
Mean: 50.008328894275, expected 50 StdDev: 3.2374717425824, expected 3.1622776601684 41 3 42 = 8 43 == 11 44 ==== 22 45 ======= 38 46 ============ 60 47 ============== 73 48 ================== 92 49 ======================= 118 50 =========================== 136 51 ========================= 128 52 ================= 89 53 ================= 89 54 =========== 56 55 ======= 37 56 === 18 57 = 7 58 = 5 59 = 6 60 2
Maple
Maple has a built-in for sampling directly from Normal distributions: <lang maple>with(Statistics): n := 100000: X := Sample( Normal(0,1), n ); Mean( X ); StandardDeviation( X ); Histogram( X );</lang>
Mathematica
<lang Mathematica>x:= RandomReal[1] SampleNormal[n_] := (Print[#//Length, " numbers, Mean : ", #//Mean, ", StandardDeviation : ", #//StandardDeviation];
Histogram[#, BarOrigin -> Left,Axes -> False])& [(Table[(-2*Log[x])^0.5*Cos[2*Pi*x], {n} ]]
Invocation:
SampleNormal[ 10000 ]
->10000 numbers, Mean : -0.0122308, StandardDeviation : 1.00646
</lang>
MATLAB / Octave
<lang Matlab> N = 100000;
x = randn(N,1); mean(x) std(x) [nn,xx] = hist(x,100); bar(xx,nn);</lang>
PARI/GP
<lang parigp>rnormal()={ my(u1=random(1.),u2=random(1.); sqrt(-2*log(u1))*cos(2*Pi*u1) \\ Could easily be extended with a second normal at very little cost. }; mean(v)={
sum(i=1,#v,v[i])/#v
}; stdev(v,mu="")={
if(mu=="",mu=mean(v)); sqrt(sum(i=1,#v,(v[i]-mu)^2))/#v
}; histogram(v,bins=16,low=0,high=1)={
my(u=vector(bins),width=(high-low)/bins); for(i=1,#v,u[(v[i]-low)\width+1]++); u
}; show(n)={
my(v=vector(n,i,rnormal()),m=mean(v),s=stdev(v,m),h,sz=ceil(n/300));
h=histogram(v,,vecmin(v)-.1,vecmax(v)+.1);
for(i=1,#h,for(j=1,h[i]\sz,print1("#"));print());
}; show(10^4)</lang>
For versions before 2.4.3, define <lang parigp>rreal()={
my(pr=32*ceil(default(realprecision)*log(10)/log(4294967296))); \\ Current precision random(2^pr)*1.>>pr
};</lang>
and use rreal() in place of random(1.).
A PARI implementation: <lang C>GEN rnormal(long prec) { pari_sp ltop = avma; GEN u1, u2, left, right, ret; u1 = randomr(prec); u2 = randomr(prec); left = sqrtr_abs(shiftr(mplog(u1), 1)); right = mpcos(mulrr(shiftr(mppi(prec), 1), u2));
ret = mulrr(left, right);
ret = gerepileupto(ltop, ret);
return ret;
}</lang>
Use mpsincos and caching to generate two values at nearly the same cost.
Pascal
//not neccessary include unit math if using function rnorm
got some trouble with math.randg needs this call randg(cMean,cMean*cStdDiv), whereas randg(cMean,cStdDiv) to get the same results??
From Free Pascal Docs unit math <lang pascal>Program Example40; {$IFDEF FPC}
{$MOde objFPC}
{$ENDIF} { Program to demonstrate the randg function. } Uses Math;
type
tTestData = extended;//because of math.randg
ttstfunc = function (mean, sd: tTestData): tTestData;
tExArray = Array of tTestData;
tSolution = record
SolExArr : tExArray;
SollowVal,
SolHighVal,
SolMean,
SolStdDiv : tTestData;
SolSmpCnt : LongInt;
end;
function getSol(genFunc:ttstfunc;Mean,StdDiv: tTestData;smpCnt: LongInt): tSolution; var
GenValue, sumValue, sumsqrVal : extended;
Begin
with result do
Begin
SolSmpCnt := smpCnt;
SolMean := 0;
SolStdDiv := 0;
SolLowVal := Mean+50* StdDiv;
SolHighVal := Mean-50* StdDiv;
setlength(SolExArr,smpCnt);
if smpCnt <= 0 then
EXIT;
sumValue := 0;
sumsqrVal := 0;
repeat
GenValue := genFunc(Mean,StdDiv);
sumValue := sumvalue+GenValue;
sumsqrVal := sumsqrVal+sqr(GenValue);
IF GenValue < SollowVal then
SollowVal:= GenValue
else
IF GenValue > SolHighVal then
SolHighVal := GenValue;
dec(smpCnt);
SolExArr[smpCnt] := GenValue;
until smpCnt<= 0;
SolMean := sumValue/SolSmpCnt;
SolStdDiv := sqrt(sumsqrVal/SolSmpCnt-sqr(SolMean));
end;
end;
//http://wiki.freepascal.org/Generating_Random_Numbers#Normal_.28Gaussian.29_Distribution function rnorm (mean, sd: tTestData): tTestData;
{Calculates Gaussian random numbers according to the Box-Müller approach}
var
u1, u2: extended;
begin
u1 := random;
u2 := random;
rnorm := mean * abs(1 + sqrt(-2 * (ln(u1))) * cos(2 * pi * u2) * sd);
end;
procedure Histo(const sol:TSolution;Colcnt,ColLen :LongInt); var
CntHisto : array of integer; LoLmt,HiLmt,span : tTestData; i, j,cnt,maxCnt: LongInt; sCross : Ansistring;
Begin
setlength(CntHisto,Colcnt);
with Sol do
Begin
span := solHighVal-solLowVal;
LoLmt := solLowVal;
writeln('Count: ',SolSmpCnt:10,' Mean ',SolMean:10:6,' StdDiv ',SolStdDIv:10:6);
writeln('span : ',span:10:5,' Low ',solLowVal:10:6,' high ',solHighVal:10:6);
end;
maxCnt := 0;
For j := 0 to Colcnt-1 do
Begin
HiLmt:= LoLmt+span/Colcnt;
cnt := 0;
with sol do
For i := 0 to High(SolExArr) do
IF (HiLmt > SolExArr[i]) AND (SolExArr[i]>= LoLmt) then
inc(cnt);
CntHisto[j] := cnt;
IF maxCnt < cnt then
maxCnt := cnt;
LoLmt:= HiLmt;
end;
inc(CntHisto[Colcnt]); // for HiLmt itself
writeln;
LoLmt := sol.solLowVal;
For i := 0 to Colcnt-1 do
Begin
Writeln(LoLmt:8:4,': ');
cnt:= Round(CntHisto[i]*ColLen/maxCnt);
setlength(sCross,cnt+3);
fillChar(sCross[1],3,' ');
fillChar(sCross[4],cnt,'#');
writeln(CntHisto[i]:10,sCross);
LoLmt := LoLmt+span/Colcnt;
end;
Writeln(sol.solHighVal:8:4,': ');
end;
const
cHistCnt = 11; cColLen = 65;
cStdDiv = 0.25; cMean = 20*cStdDiv;
var
mySol : tSolution;
begin
Randomize;
// test of randg of unit math
Writeln('function randg');
mySol := getSol(@randg,cMean,cMean*cStdDiv,100000);
Histo(mySol,cHistCnt,cColLen);
writeln;
// test of rnorm from wiki
Writeln('function rnorm');
mySol := getSol(@rnorm,cMean,cStdDiv,1000000);
Histo(mySol,cHistCnt,cColLen);
end.</lang>
- Output:
function randg Count: 100000 Mean 5.000326 StdDiv 1.250027 span : 10.65123 Low -0.333310 high 10.317922
-0.3333: 25 0.6350: 287 # 1.6033: 2291 ##### 2.5716: 9531 ##################### 3.5399: 22608 ################################################# 4.5082: 29953 ################################################################# 5.4765: 22917 ################################################## 6.4447: 9716 ##################### 7.4130: 2352 ##### 8.3813: 295 # 9.3496: 24 10.3179:function rnorm Count: 1000000 Mean 4.998391 StdDiv 1.251103 span : 11.08994 Low 0.001521 high 11.091461
0.0015: 704 1.0097: 7797 ## 2.0179: 49235 ########### 3.0261: 162761 #################################### 4.0342: 293242 ################################################################# 5.0424: 285818 ############################################################### 6.0506: 150781 ################################# 7.0588: 42641 ######### 8.0669: 6467 # 9.0751: 528 10.0833: 25 11.0915:
Perl 6
<lang perl6>constant τ = 2 * pi;
sub normdist ($m, $σ) {
my $r = sqrt -2 * log rand; my $Θ = τ * rand; $r * cos($Θ) * $σ + $m;
}
sub MAIN ($size = 100000, $mean = 50, $stddev = 4) {
my @dataset = normdist($mean,$stddev) xx $size;
my $m = [+](@dataset) / $size; say (:$m);
my $σ = sqrt [+](@dataset X** 2) / $size - $m**2; say (:$σ);
(my %hash){.round}++ for @dataset;
my $scale = 180 * $stddev / $size;
constant @subbar = < ⎸ ▏ ▎ ▍ ▌ ▋ ▊ ▉ █ >;
for %hash.keys».Int.minmax(+*) -> $i {
my $x = (%hash{$i} // 0) * $scale;
my $full = floor $x;
my $part = 8 * ($x - $full);
say $i, "\t", '█' x $full, @subbar[$part];
}
}</lang>
- Output:
"m" => 50.006107405837142e0 "σ" => 4.0814435639885254e0 33 ⎸ 34 ⎸ 35 ⎸ 36 ▏ 37 ▎ 38 ▊ 39 █▋ 40 ███⎸ 41 █████▊ 42 ██████████⎸ 43 ███████████████▋ 44 ███████████████████████▏ 45 ████████████████████████████████▌ 46 ███████████████████████████████████████████▍ 47 ██████████████████████████████████████████████████████▏ 48 ███████████████████████████████████████████████████████████████▏ 49 █████████████████████████████████████████████████████████████████████▋ 50 ███████████████████████████████████████████████████████████████████████▊ 51 █████████████████████████████████████████████████████████████████████▌ 52 ███████████████████████████████████████████████████████████████⎸ 53 ██████████████████████████████████████████████████████▎ 54 ███████████████████████████████████████████⎸ 55 ████████████████████████████████▌ 56 ███████████████████████▍ 57 ███████████████▉ 58 █████████▉ 59 █████▍ 60 ███▍ 61 █▋ 62 ▊ 63 ▍ 64 ▏ 65 ⎸ 66 ⎸ 67 ⎸
Phix
<lang Phix>procedure sample(integer n) -- show mean, standard deviation. Find max, min. sequence dat = repeat(0,n)
for i=1 to n do
dat[i] = sqrt(-2*log(rnd()))*cos(2*PI*rnd())
end for
printf(1,"%d data terms used.\n",{n})
atom mean = sum(dat)/n,
mx = max(dat),
mn = min(dat),
range = mx-mn
printf(1,"Largest term was %g & smallest was %g\n",{mx,mn})
printf(1,"Mean = %g\n",{mean})
printf(1,"Stddev = %g\n",sqrt(sum(sq_mul(dat,dat))/n-mean*mean))
-- show histogram
integer nBins = 50
sequence bins = repeat(0,nBins+1)
for i=1 to n do
bins[floor((dat[i]-mn)/range*nBins)+1] += 1
end for
for b=1 to nBins do
puts(1,repeat('#',floor(nBins*bins[b]/n*30))&"\n")
end for
end procedure
sample(100000)</lang>
- Output:
100000 data terms used. Largest term was 4.30779 & smallest was -4.11902 Mean = -0.00252597 Stddev = 1.00067 # ## #### ###### ########## ############# ################## ######################## ################################# ######################################## #################################################### ############################################################# ###################################################################### ############################################################################### ####################################################################################### ############################################################################################### ################################################################################################# ##################################################################################################### ################################################################################################### ################################################################################################ ######################################################################################## ############################################################################### ####################################################################### ############################################################ ################################################# ####################################### ############################## ######################### ################ ############ ######### ###### #### ## #
<lang Phix>function gaussian(atom mean, atom variance)
return sqrt(-2 * variance * log(rnd())) *
cos(2 * variance * PI * rnd()) + mean
end function
function mean(sequence t)
return sum(t)/length(t)
end function
function std(sequence t)
atom squares = 0,
avg = mean(t)
for i=1 to length(t) do
squares += power(avg-t[i],2)
end for
atom variance = squares/length(t)
return sqrt(variance)
end function
procedure showHistogram(sequence t)
for i=ceil(min(t)) to floor(max(t)) do
integer n = 0
for k=1 to length(t) do
n += ceil(t[k]-0.5)=i
end for
integer l = floor(n/length(t)*200)
printf(1,"%d %s %d\n",{i,repeat('=',l),n})
end for
end procedure
sequence t = repeat(0,100000) integer avg = 50, variance = 10 for i=1 to length(t) do
t[i] = gaussian(avg, variance)
end for printf(1,"Mean: %g, expected %g\n",{mean(t),avg}) printf(1,"StdDev: %g, expected %g\n",{std(t),sqrt(variance)}) showHistogram(t)</lang>
- Output:
Mean: 50.0041, expected 50 StdDev: 3.1673, expected 3.16228 37 2 38 7 39 30 40 92 41 220 42 = 523 43 == 1098 44 ==== 2140 45 ======= 3690 46 =========== 5753 47 =============== 7906 48 ==================== 10299 49 ======================= 11813 50 ========================= 12555 51 ======================= 11934 52 ==================== 10327 53 ================ 8099 54 =========== 5733 55 ======= 3684 56 ==== 2126 57 == 1098 58 487 59 226 60 106 61 36 62 9 63 7
PureBasic
<lang purebasic>Procedure.f randomf(resolution = 2147483647)
ProcedureReturn Random(resolution) / resolution
EndProcedure
Procedure.f normalDist() ;Box Muller method
ProcedureReturn Sqr(-2 * Log(randomf())) * Cos(2 * #PI * randomf())
EndProcedure
Procedure sample(n, nBins = 50)
Protected i, maxBinValue, binNumber
Protected.f d, mean, sum, sumSq, mx, mn, range
Dim dat.f(n)
For i = 1 To n
dat(i) = normalDist()
Next
;show mean, standard deviation, find max & min.
mx = -1000
mn = 1000
sum = 0
sumSq = 0
For i = 1 To n
d = dat(i)
If d > mx: mx = d: EndIf
If d < mn: mn = d: EndIf
sum + d
sumSq + d * d
Next
PrintN(Str(n) + " data terms used.")
PrintN("Largest term was " + StrF(mx) + " & smallest was " + StrF(mn))
mean = sum / n
PrintN("Mean = " + StrF(mean))
PrintN("Stddev = " + StrF((sumSq / n) - Sqr(mean * mean)))
;show histogram
range = mx - mn
Dim bins(nBins)
For i = 1 To n
binNumber = Int(nBins * (dat(i) - mn) / range)
bins(binNumber) + 1
Next
maxBinValue = 1
For i = 0 To nBins
If bins(i) > maxBinValue
maxBinValue = bins(i)
EndIf
Next
#normalizedMaxValue = 70
For binNumber = 0 To nBins
tickMarks = Round(bins(binNumber) * #normalizedMaxValue / maxBinValue, #PB_Round_Nearest)
PrintN(ReplaceString(Space(tickMarks), " ", "#"))
Next
PrintN("")
EndProcedure
If OpenConsole()
sample(100000) Print(#CRLF$ + #CRLF$ + "Press ENTER to exit"): Input() CloseConsole()
EndIf</lang> Sample output:
100000 data terms used. Largest term was 4.5352029800 & smallest was -4.5405135155 Mean = 0.0012346541 Stddev = 0.9959455132 # ### ###### ########## ################## ############################ ######################################### ##################################################### ################################################################ ###################################################################### ###################################################################### ################################################################ ##################################################### ######################################### ############################# ################## ########## ###### ### #
Python
This uses the external matplotlib package as well as the built-in standardlib function random.gauss. <lang python>from __future__ import division import matplotlib.pyplot as plt import random
mean, stddev, size = 50, 4, 100000 data = [random.gauss(mean, stddev) for c in range(size)]
mn = sum(data) / size sd = (sum(x*x for x in data) / size
- (sum(data) / size) ** 2) ** 0.5
print("Sample mean = %g; Stddev = %g; max = %g; min = %g for %i values"
% (mn, sd, max(data), min(data), size))
plt.hist(data,bins=50)</lang>
- Output:
Sample mean = 49.9822; Stddev = 4.00938; max = 66.8091; min = 33.5283 for 100000 values
R
R can generate random normal distributed numbers using the rnorm command: <lang r>n = 100000; X = rnorm(n, mean = 0, sd = 1); mean( X ); sd( X ); hist( X );</lang>
Racket
This shows how one would generate samples from a normal distribution, compute statistics and plot a histogram.
<lang racket>
- lang racket
(require math (planet williams/science/histogram-with-graphics))
(define data (sample (normal-dist 50 4) 100000))
(displayln (~a "Mean:\t" (mean data))) (displayln (~a "Stddev:\t" (stddev data))) (displayln (~a "Max:\t" (apply max data))) (displayln (~a "Min:\t" (apply min data)))
(define h (make-histogram-with-ranges-uniform 40 30 70)) (for ([x data]) (histogram-increment! h x)) (histogram-plot h "Normal distribution μ=50 σ=4") </lang>
The other part of the task was to produce normal distributed numbers from a unit distribution. The following code is an implementation of the polar method. It is a slightly modified version of code originally written by Sebastian Egner. <lang racket>
- lang racket
(require math)
(define random-normal
(let ([unit (uniform-dist)]
[next #f])
(λ (μ σ)
(if next
(begin0
(+ μ (* σ next))
(set! next #f))
(let loop ()
(let* ([v1 (- (* 2.0 (sample unit)) 1.0)]
[v2 (- (* 2.0 (sample unit)) 1.0)]
[s (+ (sqr v1) (sqr v2))])
(cond [(>= s 1) (loop)]
[else (define scale (sqrt (/ (* -2.0 (log s)) s)))
(set! next (* scale v2))
(+ μ (* σ scale v1))])))))))
</lang>
REXX
The REXX language doesn't have any "higher math" BIF functions like SIN/COS/LN/LOG/SQRT/POW/etc,
so we hoi polloi programmers have to roll our own.
<lang rexx>/*REXX program generates 10,000 normally distributed numbers (Gaussian distribution).*/
parse arg n seed . /*obtain optional arguments from the CL*/
if n== | n=="," then n=10000 /*Not specified? Then use the default.*/
if datatype(seed,'W') then call random ,,seed /*seed is for repeatable RANDOM numbers*/
call pi /*call subroutine to define pi constant*/
do g=1 for n /*generate N uniform random numbers. */
#.g=sqrt(-2*ln(rand()))*cos(2*pi*rand()) /*assign a uniform random number to #. */
end /*g*/
mn=#.1; mx=mn; s=0; ss=0; noise=n*.0005 /*calculate the noise: 1/20th % of N.*/
do j=1 for n; _=#.j; s=s+_; ss=ss+_*_ /*the sum, and the sum of squares. */
mn=min(mn,#.j); mx=max(mx,#.j) /*find the minimum and the maximum. */
end /*j*/
!.=0 say 'number of data points = ' aa(n ) say ' minimum = ' aa(mn ) say ' maximum = ' aa(mx ) say ' arithmetic mean = ' aa(s/n) say ' standard deviation = ' aa(sqrt(ss/n - (s/n)**2) ) ?mn=!.1; ?mx=?mn /*define minimum & maximum value so far*/ parse value scrSize() with sd sw . /*obtain the (true) screen size of term*/ /*◄──not all REXXes have this BIF*/ sdE=sd-4 /*the effective (useable) screen depth.*/ swE=sw-1 /* " " " " width.*/ $=1 / max(1,mx-mn) * sdE /*$ is used for scaling depth of histo*/
do i=1 for n; ?=trunc((#.i-mn) * $) /*calculate the relative line.*/
!.?=!.?+1 /*bump the counter. */
?mn=min(?mn, !.?); ?mx=max(?mx, !.?) /*find the minimum and maximum*/
end /*i*/
f=swE/?mx /*limit graph to 1 full screen*/
do h=0 for sdE; _=!.h /*obtain a data point. */
if _>noise then say copies('─', trunc(_*f)) /*display a bar of histogram. */
end /*h*/ /*[↑] use a hyphen for histo.*/
exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ aa: parse arg a; return left(,(a>=0)+2*datatype(a,'W'))a /*prepend a blank if #>=0, add two blanks if its whole*/ e: e =2.7182818284590452353602874713526624977572470936999595749669676277240766303535; return e pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862; return pi rand: return random(, 1e5) / 1e5 /*REXX generates uniform random postive integers.*/ r2r: return arg(1) // (pi()*2) /*normalize the given angle (in radians) to ±2pi.*/ .sincos: parse arg z,_,i; x=x*x; p=z; do k=2 by 2; _=-_*x/(k*(k+i)); z=z+_; if z=p then leave; p=z; end; return z ln: procedure; parse arg x,f; call e; ig=x>1.5; is=1-2*(ig\==1); ii=0; xx=x; return .ln() .ln: do while ig&xx>1.5|\ig&xx<.5;_=e; do k=-1;iz=xx*_**-is;if k>=0&(ig&iz<1|\ig&iz>.5) then leave; _=_*_;izz=iz; end
xx=izz;ii=ii+is*2**k;end; x=x*e**-ii-1;z=0;_=-1;p=z; do k=1;_=-_*x;z=z+_/k;if z=p then leave;p=z;end; return z+ii
cos: procedure; parse arg x; x=r2r(x); a=abs(x); hpi=pi*.5; numeric fuzz min(6,digits()-3); if a=pi() then return -1
if a=hpi|a=hpi*3 then return 0; if a=pi()/3 then return .5; if a=pi()*2/3 then return -.5; return .sinCos(1,1,-1)
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); m.=9; numeric digits; numeric form; h=d+6
parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g*.5'e'_%2; do j=0 while h>9; m.j=h; h=h%2+1; end /*j*/
do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/; numeric digits d; return g/1</lang>
This REXX program makes use of scrsize REXX program (or BIF) which is used to determine the screen size of the terminal (console).
The SCRSIZE.REX REXX program is included here ──► SCRSIZE.REX.
output when using the default input:
(The output shown when the screen size is 50x80.)
number of data points = 10000
minimum = -3.41571894
maximum = 3.96752904
arithmetic mean = -0.0150910306
standard deviation = 0.99056458
─
─
───
─────
───────
───────
────────────
─────────────────
─────────────────────
─────────────────────────────
────────────────────────────────────
───────────────────────────────────────────
───────────────────────────────────────────────
──────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────
───────────────────────────────────────────────────
─────────────────────────────────────────────
───────────────────────────────────────
───────────────────────────────
─────────────────────────────
──────────────────
─────────────
────────
───────
────
───
─
─
output when using the default input:
(The output shown when the screen size is 60x130.)
number of data points = 10000
minimum = -3.83073183
maximum = 3.61051026
arithmetic mean = 0.00421997333
standard deviation = 0.981924955
──
──
────
──
─────
────────
────────
────────
───────────
──────────────
──────────────────
────────────────────
────────────────────────────
─────────────────────────────────
────────────────────────────────
─────────────────────────────────────────────
───────────────────────────────────────────────
───────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────
───────────────────────────────────
──────────────────────────────────────
───────────────────────────────────
───────────────────────────
───────────────────
──────────────────────
─────────────
───────────
───────
──────
───────
───
────
───
──
Run BASIC
<lang runbasic> s = 100000 h$ = "=============================================================" h$ = h$ + h$ dim ndis(s) ' mean and standard deviation. mx = -9999 mn = 9999 sum = 0 sumSqr = 0 for i = 1 to s ' find minimum and maximum ms = rnd(1) ss = rnd(1) nd = (-2 * log(ms))^0.5 * cos(2 *3.14159265 * ss) ' normal distribution ndis(i) = nd mx = max(mx, nd) mn = min(mn, nd) sum = sum + nd sumSqr = sumSqr + nd ^ 2 next i
mean = sum / s range = mx - mn
print "Samples :"; s print "Largest :"; mx print "Smallest :"; mn print "Range :"; range print "Mean :"; mean print "Stand Dev :"; (sumSqr /s -mean^2)^0.5
'Show chart of histogram nBins = 50 dim bins(nBins) for i = 1 to s z = int((ndis(i) -mn) /range *nBins) bins(z) = bins(z) + 1 mb = max(bins(z),mb) next i for b = 0 to nBins -1
print using("##",b);" ";using("#####",bins(b));" ";left$(h$,(bins(b) / mb) * 90)
next b END</lang>
- Output:
Samples :100000 Largest :4.61187177 Smallest :-4.21695424 Range :8.82882601 Mean :-9.25042513e-4 Stand Dev :1.00680067 = == === ===== ======== ============= ================= ======================= ============================== ======================================= =============================================== ========================================================= =================================================================== =========================================================================== =================================================================================== ======================================================================================= ========================================================================================== ======================================================================================== ====================================================================================== ================================================================================= ============================================================================ ================================================================== ======================================================== ============================================== ===================================== ============================ ===================== =============== ========== ======= ===== === = =
SAS
<lang sas>data test; n=100000; twopi=2*constant('pi'); do i=1 to n; u=ranuni(0); v=ranuni(0); r=sqrt(-2*log(u)); x=r*cos(twopi*v); y=r*sin(twopi*v); z=rannor(0); output; end; keep x y z;
proc means mean stddev;
proc univariate; histogram /normal;
run;
/* Variable Mean Std Dev
x -0.0052720 0.9988467 y 0.000023995 1.0019996 z 0.0012857 1.0056536
- /</lang>
Sidef
<lang ruby>define τ = Number.tau
func normdist (m, σ) {
var r = sqrt(-2 * 1.rand.log) var Θ = (τ * 1.rand) r * Θ.cos * σ + m
}
var size = 100_000 var mean = 50 var stddev = 4
var dataset = size.of { normdist(mean, stddev) } var m = (dataset.sum(0) / size) say ("m: #{m}")
var σ = sqrt(dataset »**» 2 -> sum(0) / size - m**2) say ("s: #{σ}")
var hash = Hash() dataset.each { |n| hash{ n.round(0) } := 0 ++ }
var scale = (180 * stddev / size) const subbar = < ⎸ ▏ ▎ ▍ ▌ ▋ ▊ ▉ █ >
for i in (hash.keys.map{.to_i}.sort) {
var x = (hash{i} * scale)
var full = x.int
var part = (8 * (x - full))
say (i, "\t", '█' * full, subbar[part])
}</lang>
- Output:
m: 49.99538275618550306540055142077589 s: 4.00295544816687358837821680496471 33 ⎸ 34 ⎸ 35 ⎸ 36 ▏ 37 ▎ 38 ▊ 39 █▋ 40 ███▏ 41 ██████▏ 42 █████████▍ 43 ███████████████▌ 44 ███████████████████████▋ 45 ████████████████████████████████▍ 46 ████████████████████████████████████████████▎ 47 █████████████████████████████████████████████████████▍ 48 ███████████████████████████████████████████████████████████████▍ 49 █████████████████████████████████████████████████████████████████████▌ 50 ████████████████████████████████████████████████████████████████████████▋ 51 █████████████████████████████████████████████████████████████████████▊ 52 ██████████████████████████████████████████████████████████████▏ 53 ████████████████████████████████████████████████████▉ 54 ███████████████████████████████████████████▉ 55 █████████████████████████████████▎ 56 ███████████████████████⎸ 57 ███████████████▋ 58 █████████▋ 59 █████▍ 60 ███▍ 61 █▊ 62 ▋ 63 ▍ 64 ▏ 65 ⎸ 66 ⎸
Tcl
<lang tcl>package require Tcl 8.5
- Uses the Box-Muller transform to compute a pair of normal random numbers
proc tcl::mathfunc::nrand {mean stddev} {
variable savednormalrandom
if {[info exists savednormalrandom]} {
return [expr {$savednormalrandom*$stddev + $mean}][unset savednormalrandom]
}
set r [expr {sqrt(-2*log(rand()))}]
set theta [expr {2*3.1415927*rand()}]
set savednormalrandom [expr {$r*sin($theta)}]
expr {$r*cos($theta)*$stddev + $mean}
} proc stats {size {slotfactor 10}} {
set sum 0.0
set sum2 0.0
for {set i 0} {$i < $size} {incr i} {
set r [expr { nrand(0.5, 0.2) }]
incr histo([expr {int(floor($r*$slotfactor))}]) set sum [expr {$sum + $r}] set sum2 [expr {$sum2 + $r**2}]
}
set mean [expr {$sum / $size}]
set stddev [expr {sqrt($sum2/$size - $mean**2)}]
puts "$size numbers"
puts "Mean: $mean"
puts "StdDev: $stddev"
foreach i [lsort -integer [array names histo]] {
puts [string repeat "*" [expr {$histo($i)*350/int($size)}]]
}
}
stats 100 puts "" stats 1000 puts "" stats 10000 puts "" stats 100000 20</lang> Sample output:
100 numbers Mean: 0.49355955990390254 StdDev: 0.19651396178121985 *** ******* ************** *********************************** ******************************************************** ****************************************************************** ************************************************************************* ****************************************** ************************************** ************** 1000 numbers Mean: 0.5066940614105869 StdDev: 0.2016794788065389 * ***** ************** **************************** ********************************************************** **************************************************************** ************************************************************* ****************************************************** *********************************** ************ ********* * 10000 numbers Mean: 0.49980964730768285 StdDev: 0.1968441612522318 * ***** *************** ******************************* ***************************************************** ****************************************************************** ******************************************************************* **************************************************** ********************************* *************** ***** * 100000 numbers Mean: 0.49960438950922254 StdDev: 0.20060211160998606 * ** *** ****** ********* ************** ****************** *********************** ***************************** ******************************** ********************************** ********************************** ******************************** **************************** *********************** ****************** ************* ********* ****** *** ** *
The blank lines in the output are where the number of samples is too small to even merit a single unit on the histogram.
zkl
<lang zkl>fcn norm2{ // Box-Muller
const PI2=(0.0).pi*2;; rnd:=(0.0).random.fp(1); // random number in [0,1), using partial application r,a:=(-2.0*rnd().log()).sqrt(), PI2*rnd(); return(r*a.cos(), r*a.sin()); // z0,z1
} const N=100000, BINS=12, SIG=3, SCALE=500; var sum=0.0,sumSq=0.0, h=BINS.pump(List(),0); // (0,0,0,...) fcn accum(v){
sum+=v; sumSq+=v*v; b:=(v + SIG)*BINS/SIG/2; if(0<=b<BINS) h[b]+=1;
};</lang> Partial application: rnd() --> (0.0).random(1). Basically, the fp method fixes the call parameters, which are then used when the partial thing is run. <lang zkl>foreach i in (N/2){ v1,v2:=norm2(); accum(v1); accum(v2); } println("Samples: %,d".fmt(N)); println("Mean: ", m:=sum/N); println("Stddev: ", (sumSq/N - m*m).sqrt()); foreach p in (h){ println("*"*(p/SCALE)) }</lang>
- Output:
Samples: 100,000 Mean: 0.0005999 Stddev: 1.003 * *** ******** ****************** ***************************** ************************************** ************************************** ***************************** ****************** ******** *** *