Statistics/Normal distribution: Difference between revisions
(Add R implementation) |
(โโ{{header|Tcl}}: added zkl) |
||
| Line 1,851: | Line 1,851: | ||
</pre> |
</pre> |
||
The blank lines in the output are where the number of samples is too small to even merit a single unit on the histogram. |
The blank lines in the output are where the number of samples is too small to even merit a single unit on the histogram. |
||
=={{header|zkl}}== |
|||
{{trans|Go}} |
|||
<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> |
|||
{{out}} |
|||
<pre> |
|||
Samples: 100,000 |
|||
Mean: 0.0005999 |
|||
Stddev: 1.003 |
|||
* |
|||
*** |
|||
******** |
|||
****************** |
|||
***************************** |
|||
************************************** |
|||
************************************** |
|||
***************************** |
|||
****************** |
|||
******** |
|||
*** |
|||
* |
|||
</pre> |
|||
Revision as of 08:50, 1 September 2016
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++
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:
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 * **** ********* *************** ******************* ****************** ************** ********* **** *
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 โธ
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)) r=mx-mn /*is used for scaling the histogram. */ 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.*/ ?mn=; ?mx=
do i=1 for n; ?=trunc((#.i-mn)/r*sdE)
!.?=!.?+1 /*bump the counter. */
if ?mn== then do; ?mn=!.j; ?mx=!.j; end /*define min, max (1st).*/
?mn=min(?mn, !.?); ?mx=max(?mx, !.?) /*find the min and max. */
end /*i*/
f=swE/?mx /*limit graph to 1 screen*/
do h=0 for sdE; _=!.h /*obtain a data point. */
if _>noise then say copies('โ', trunc(_*f)) /*display the histogram. */
end /*h*/
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(0,1e5) / 1e5 /*REXX generates uniform random postive integers.*/ r2r: return arg(1) // (2*pi()) /*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_comp() .ln_comp: 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 * *** ******** ****************** ***************************** ************************************** ************************************** ***************************** ****************** ******** *** *