Statistics/Normal distribution: Difference between revisions
Thundergnat (talk | contribs) (Rename Perl 6 -> Raku, alphabetize, minor clean-up) |
SqrtNegInf (talk | contribs) m (→{{header|Perl}}: Fix link: Perl 6 --> Raku) |
||
Line 1,918:
=={{header|Perl}}==
{{trans|
<lang perl>use constant pi => 3.14159265;
use List::Util qw(sum reduce min max);
|
Revision as of 15:29, 15 March 2020
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 standard deviation, and show a histogram of the data.
- 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 C>/*
* RosettaCode example: Statistics/Normal distribution in C * * The random number generator rand() of the standard C library is obsolete * and should not be used in more demanding applications. There are plenty * libraries with advanced features (eg. GSL) with functions to calculate * the mean, the standard deviation, generating random numbers etc. * However, these features are not the core of the standard C library. */
- include <stdio.h>
- include <stdlib.h>
- include <math.h>
- include <string.h>
- include <time.h>
- define NMAX 10000000
double mean(double* values, int n)
{
int i; double s = 0;
for ( i = 0; i < n; i++ ) s += values[i]; return s / n;
}
double stddev(double* values, int n)
{
int i; double average = mean(values,n); double s = 0;
for ( i = 0; i < n; i++ ) s += (values[i] - average) * (values[i] - average); return sqrt(s / (n - 1));
}
/*
* Normal random numbers generator - Marsaglia algorithm. */
double* generate(int n) {
int i; int m = n + n % 2; double* values = (double*)calloc(m,sizeof(double)); double average, deviation;
if ( values ) { for ( i = 0; i < m; i += 2 ) { double x,y,rsq,f; do { x = 2.0 * rand() / (double)RAND_MAX - 1.0; y = 2.0 * rand() / (double)RAND_MAX - 1.0; rsq = x * x + y * y; }while( rsq >= 1. || rsq == 0. ); f = sqrt( -2.0 * log(rsq) / rsq ); values[i] = x * f; values[i+1] = y * f; } } return values;
}
void printHistogram(double* values, int n)
{
const int width = 50; int max = 0;
const double low = -3.05; const double high = 3.05; const double delta = 0.1;
int i,j,k; int nbins = (int)((high - low) / delta); int* bins = (int*)calloc(nbins,sizeof(int)); if ( bins != NULL ) { for ( i = 0; i < n; i++ ) { int j = (int)( (values[i] - low) / delta ); if ( 0 <= j && j < nbins ) bins[j]++; }
for ( j = 0; j < nbins; j++ ) if ( max < bins[j] ) max = bins[j];
for ( j = 0; j < nbins; j++ ) { printf("(%5.2f, %5.2f) |", low + j * delta, low + (j + 1) * delta ); k = (int)( (double)width * (double)bins[j] / (double)max ); while(k-- > 0) putchar('*'); printf(" %-.1f%%", bins[j] * 100.0 / (double)n); putchar('\n'); }
free(bins); }
}
int main(void)
{
double* seq;
srand((unsigned int)time(NULL));
if ( (seq = generate(NMAX)) != NULL ) { printf("mean = %g, stddev = %g\n\n", mean(seq,NMAX), stddev(seq,NMAX)); printHistogram(seq,NMAX); free(seq);
printf("\n%s\n", "press enter"); getchar(); return EXIT_SUCCESS; } return EXIT_FAILURE;
}</lang>
- Output:
mean = 0.000477941, stddev = 0.999945 (-3.05, -2.95) | 0.1% (-2.95, -2.85) | 0.1% (-2.85, -2.75) |* 0.1% (-2.75, -2.65) |* 0.1% (-2.65, -2.55) |* 0.1% (-2.55, -2.45) |** 0.2% (-2.45, -2.35) |** 0.2% (-2.35, -2.25) |*** 0.3% (-2.25, -2.15) |**** 0.4% (-2.15, -2.05) |***** 0.4% (-2.05, -1.95) |****** 0.5% (-1.95, -1.85) |******** 0.7% (-1.85, -1.75) |********* 0.8% (-1.75, -1.65) |*********** 0.9% (-1.65, -1.55) |************* 1.1% (-1.55, -1.45) |**************** 1.3% (-1.45, -1.35) |****************** 1.5% (-1.35, -1.25) |********************* 1.7% (-1.25, -1.15) |************************ 1.9% (-1.15, -1.05) |*************************** 2.2% (-1.05, -0.95) |****************************** 2.4% (-0.95, -0.85) |********************************* 2.7% (-0.85, -0.75) |************************************ 2.9% (-0.75, -0.65) |*************************************** 3.1% (-0.65, -0.55) |***************************************** 3.3% (-0.55, -0.45) |******************************************** 3.5% (-0.45, -0.35) |********************************************** 3.7% (-0.35, -0.25) |*********************************************** 3.8% (-0.25, -0.15) |************************************************* 3.9% (-0.15, -0.05) |************************************************* 4.0% (-0.05, 0.05) |************************************************** 4.0% ( 0.05, 0.15) |************************************************* 4.0% ( 0.15, 0.25) |************************************************* 3.9% ( 0.25, 0.35) |*********************************************** 3.8% ( 0.35, 0.45) |********************************************** 3.7% ( 0.45, 0.55) |******************************************** 3.5% ( 0.55, 0.65) |***************************************** 3.3% ( 0.65, 0.75) |*************************************** 3.1% ( 0.75, 0.85) |************************************ 2.9% ( 0.85, 0.95) |********************************* 2.7% ( 0.95, 1.05) |****************************** 2.4% ( 1.05, 1.15) |*************************** 2.2% ( 1.15, 1.25) |************************ 1.9% ( 1.25, 1.35) |********************* 1.7% ( 1.35, 1.45) |****************** 1.5% ( 1.45, 1.55) |**************** 1.3% ( 1.55, 1.65) |************* 1.1% ( 1.65, 1.75) |*********** 0.9% ( 1.75, 1.85) |********* 0.8% ( 1.85, 1.95) |******** 0.7% ( 1.95, 2.05) |****** 0.5% ( 2.05, 2.15) |***** 0.4% ( 2.15, 2.25) |**** 0.4% ( 2.25, 2.35) |*** 0.3% ( 2.35, 2.45) |** 0.2% ( 2.45, 2.55) |** 0.2% ( 2.55, 2.65) |* 0.1% ( 2.65, 2.75) |* 0.1% ( 2.75, 2.85) |* 0.1% ( 2.85, 2.95) | 0.1% press enter
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: **
Julia
Julia has the builtin package "Distributions" to generate random numbers from a standard distribution (Normal, Chisq etc.). <lang julia>using Printf, Distributions, Gadfly
data = rand(Normal(0, 1), 1000) @printf("N = %i\n", length(data)) @printf("μ = %2.2f\tσ = %2.2f\n", mean(data), std(data)) @printf("range = (%2.2f, %2.2f\n)", minimum(data), maximum(data)) h = plot(x=data, Geom.histogram) draw(PNG("norm_hist.png", 10cm, 10cm), h)</lang>
- Output:
N = 1000 μ = 0.02 σ = 0.97 range = (-2.76, 3.42)
Kotlin
<lang scala>// version 1.1.2
val rand = java.util.Random()
fun normalStats(sampleSize: Int) {
if (sampleSize < 1) return val r = DoubleArray(sampleSize) val h = IntArray(12) // all zero by default /* Generate 'sampleSize' normally distributed random numbers with mean 0.5 and SD 0.25 and calculate in which box they will fall when drawing the histogram */ for (i in 0 until sampleSize) { r[i] = 0.5 + rand.nextGaussian() / 4.0 when { r[i] < 0.0 -> h[0]++ r[i] >= 1.0 -> h[11]++ else -> h[1 + (r[i] * 10).toInt()]++ } }
// adjust one of the h[] values if necessary to ensure they sum to sampleSize val adj = sampleSize - h.sum() if (adj != 0) { for (i in 0..11) { h[i] += adj if (h[i] >= 0) break h[i] -= adj } }
val mean = r.average() val sd = Math.sqrt(r.map { (it - mean) * (it - mean) }.average()) // Draw a histogram of the data with interval 0.1 var numStars: Int // If sample size > 300 then normalize histogram to 300 val scale = if (sampleSize <= 300) 1.0 else 300.0 / sampleSize println("Sample size $sampleSize\n") println(" Mean ${"%1.6f".format(mean)} SD ${"%1.6f".format(sd)}\n") for (i in 0..11) { when (i) { 0 -> print("< 0.00 : ") 11 -> print(">=1.00 : ") else -> print(" %1.2f : ".format(i / 10.0)) } print("%5d ".format(h[i])) numStars = (h[i] * scale + 0.5).toInt() println("*".repeat(numStars)) } println()
}
fun main(args: Array<String>) {
val sampleSizes = intArrayOf(100, 1_000, 10_000, 100_000) for (sampleSize in sampleSizes) normalStats(sampleSize)
}</lang>
- Output:
Sample size 100 Mean 0.525211 SD 0.266316 < 0.00 : 3 *** 0.10 : 1 * 0.20 : 3 *** 0.30 : 11 *********** 0.40 : 14 ************** 0.50 : 13 ************* 0.60 : 15 *************** 0.70 : 13 ************* 0.80 : 10 ********** 0.90 : 11 *********** 1.00 : 4 **** >=1.00 : 2 ** Sample size 1000 Mean 0.500948 SD 0.255757 < 0.00 : 29 ********* 0.10 : 35 *********** 0.20 : 70 ********************* 0.30 : 71 ********************* 0.40 : 138 ***************************************** 0.50 : 139 ****************************************** 0.60 : 168 ************************************************** 0.70 : 123 ************************************* 0.80 : 110 ********************************* 0.90 : 62 ******************* 1.00 : 32 ********** >=1.00 : 23 ******* Sample size 10000 Mean 0.501376 SD 0.248317 < 0.00 : 240 ******* 0.10 : 305 ********* 0.20 : 617 ******************* 0.30 : 927 **************************** 0.40 : 1291 *************************************** 0.50 : 1554 *********************************************** 0.60 : 1609 ************************************************ 0.70 : 1319 **************************************** 0.80 : 983 ***************************** 0.90 : 609 ****************** 1.00 : 324 ********** >=1.00 : 222 ******* Sample size 100000 Mean 0.499427 SD 0.250533 < 0.00 : 2341 ******* 0.10 : 3246 ********** 0.20 : 6005 ****************** 0.30 : 9718 ***************************** 0.40 : 13247 **************************************** 0.50 : 15595 *********************************************** 0.60 : 15271 ********************************************** 0.70 : 13405 **************************************** 0.80 : 9653 ***************************** 0.90 : 5990 ****************** 1.00 : 3257 ********** >=1.00 : 2272 *******
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 * 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
<lang perl>use constant pi => 3.14159265; use List::Util qw(sum reduce min max);
sub normdist {
my($m, $sigma) = @_; my $r = sqrt -2 * log rand; my $theta = 2 * pi * rand; $r * cos($theta) * $sigma + $m;
}
$size = 100000; $mean = 50; $stddev = 4;
push @dataset, normdist($mean,$stddev) for 1..$size;
my $m = sum(@dataset) / $size; print "m = $m\n";
my $sigma = sqrt( (reduce { $a + $b **2 } 0,@dataset) / $size - $m**2 ); print "sigma = $sigma\n";
$hash{int $_}++ for @dataset; my $scale = 180 * $stddev / $size; my @subbar = < ⎸ ▏ ▎ ▍ ▌ ▋ ▊ ▉ █ >; for $i (min(@dataset)..max(@dataset)) { my $x = ($hash{$i} // 0) * $scale; my $full = int $x; my $part = 8 * ($x - $full); my $t1 = '█' x $full; my $t2 = $subbar[$part]; print "$i\t$t1$t2\n"; }
</lang>
- Output:
32 ⎸ 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 ⎸
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>
Raku
(formerly Perl 6)
<lang perl6>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 ⎸
REXX
The REXX language doesn't have any "higher math" BIF functions like SIN, COS, LN, LOG, SQRT, EXP, POW, etc,
so we hoi polloi programmers have to roll our own.
<lang rexx>/*REXX program generates 10,000 normally distributed numbers (Gaussian distribution).*/
numeric digits 20 /*use 20 decimal digits for accuracy.*/
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; #.g= sqrt( -2 * ln( rand() ) ) * cos( 2 * pi * rand() ) end /*g*/ /* [↑] uniform random number ───► #.g */
s= 0 mn= #.1; mx= mn; noise= n * .0005 /*calculate the noise: 1/20th % of N.*/ ss= 0
do j=1 for n; _=#.j; s=s+_; ss=ss+_*_ /*the sum, and the sum of squares. */ mn= min(mn, _); mx= max(mx, _) /*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 2 blanks if whole.*/ e: e = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535; return e pi: pi= 3.1415926535897932384626433832795028841971693993751058209749445923078164062862; return pi r2r: return arg(1) // (pi() * 2) /*normalize the given angle (in radians) to ±2pi.*/ rand: return random(1, 1e5) / 1e5 /*REXX generates uniform random postive integers.*/ .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; 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); this is to aid in maximizing the width of the horizontal histogram.
The SCRSIZE.REX REXX program is included here ──► SCRSIZE.REX.
- output when using the default input:
(The output shown when the screen size is 62x140.)
The graph is shown at 3/4 size.
number of data points = 10000 minimum = -3.8181072371544448250 maximum = 3.5445917138265268562 arithmetic mean = -0.01406470979976873427 standard deviation = 0.99486092191249231518 ─ ─ ─── ──── ───── ───── ──────── ─────────── ────────────── ───────────────────── ────────────────────── ────────────────────────────────── ──────────────────────────────────────── ─────────────────────────────────────────────── ───────────────────────────────────────────────────── ───────────────────────────────────────────────────────────────────────── ───────────────────────────────────────────────────────────────── ───────────────────────────────────────────────────────────────────────────────────── ────────────────────────────────────────────────────────────────────────────────────────────────── ────────────────────────────────────────────────────────────────────────────────────────────── ──────────────────────────────────────────────────────────────────────────────────────────────────────────────── ──────────────────────────────────────────────────────────────────────────────────────────────────────────────── ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ────────────────────────────────────────────────────────────────────────────────────────────────────────── ────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ──────────────────────────────────────────────────────────────────────────────────────────────────────────────── ─────────────────────────────────────────────────────────────────────────────────────────────── ──────────────────────────────────────────────────────────────────────────── ─────────────────────────────────────────────────────────────────────── ─────────────────────────────────────────────────────────────────── ───────────────────────────────────────────────────── ───────────────────────────────────────────────── ───────────────────────────────── ────────────────────────────────── ─────────────────────── ────────────────────── ────────────────── ─────────── ────────── ────── ─── ──── ── ─
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 τ = Num.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 / size) say ("m: #{m}")
var σ = sqrt(dataset »**» 2 -> sum / size - m**2) say ("s: #{σ}")
var hash = Hash() dataset.each { |n| hash{ n.round } := 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 ⎸
Stata
Pairs of normal numbers are generated from pairs of uniform numbers using the Box-Muller method. A normal density is added to the histogram for comparison. See histogram in Stata help. A Q-Q plot is also drawn.
<lang stata>clear all set obs 100000 gen u=runiform() gen v=runiform() gen r=sqrt(-2*log(u)) gen x=r*cos(2*_pi*v) gen y=r*sin(2*_pi*v) gen z=rnormal() sum x y z
Variable | Obs Mean Std. Dev. Min Max
+---------------------------------------------------------
x | 100,000 .0025861 1.002346 -4.508192 4.164336 y | 100,000 .0017389 1.001586 -4.631144 4.460274 z | 100,000 .005054 .9998861 -5.134265 4.449522
hist x, normal hist y, normal hist z, normal qqplot x z, msize(tiny)</lang>
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.
VBA
<lang vb>Public Sub standard_normal()
Dim s() As Variant, bins(71) As Single ReDim s(20000) For i = 1 To 20000 s(i) = WorksheetFunction.Norm_S_Inv(Rnd()) Next i For i = -35 To 35 bins(i + 36) = i / 10 Next i Debug.Print "sample size"; UBound(s), "mean"; mean(s), "standard deviation"; standard_deviation(s) t = WorksheetFunction.Frequency(s, bins) For i = -35 To 35 Debug.Print Format((i - 1) / 10, "0.00"); Debug.Print "-"; Format(i / 10, "0.00"), Debug.Print String$(t(i + 36, 1) / 10, "X"); Debug.Print Next i
End Sub</lang>
- Output:
sample size 20000 mean-5,26306310478751E-03 standard deviation 1,00355037427319 -3,60--3,50 -3,50--3,40 -3,40--3,30 -3,30--3,20 -3,20--3,10 -3,10--3,00 -3,00--2,90 XX -2,90--2,80 X -2,80--2,70 XX -2,70--2,60 XX -2,60--2,50 XXX -2,50--2,40 XXXX -2,40--2,30 XXXXX -2,30--2,20 XXXXXXXX -2,20--2,10 XXXXXXXX -2,10--2,00 XXXXXXXXXXX -2,00--1,90 XXXXXXXXXXXXX -1,90--1,80 XXXXXXXXXXXXXXX -1,80--1,70 XXXXXXXXXXXXXXXX -1,70--1,60 XXXXXXXXXXXXXXXXXXXX -1,60--1,50 XXXXXXXXXXXXXXXXXXXXXXXX -1,50--1,40 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -1,40--1,30 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -1,30--1,20 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -1,20--1,10 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -1,10--1,00 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -1,00--0,90 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,90--0,80 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,80--0,70 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,70--0,60 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,60--0,50 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,50--0,40 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,40--0,30 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,30--0,20 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,20--0,10 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,10-0,00 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,00-0,10 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,10-0,20 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,20-0,30 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,30-0,40 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,40-0,50 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,50-0,60 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,60-0,70 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,70-0,80 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,80-0,90 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,90-1,00 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 1,00-1,10 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 1,10-1,20 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 1,20-1,30 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 1,30-1,40 XXXXXXXXXXXXXXXXXXXXXXXXXXXXX 1,40-1,50 XXXXXXXXXXXXXXXXXXXXXXXXXX 1,50-1,60 XXXXXXXXXXXXXXXXXXXXXXXXX 1,60-1,70 XXXXXXXXXXXXXXXXXXXXXX 1,70-1,80 XXXXXXXXXXXXXXXXXX 1,80-1,90 XXXXXXXXXXXXXXX 1,90-2,00 XXXXXXXXXXX 2,00-2,10 XXXXXXXXXXXX 2,10-2,20 XXXXXXX 2,20-2,30 XXXXXX 2,30-2,40 XXXXX 2,40-2,50 XXX 2,50-2,60 XXXX 2,60-2,70 XX 2,70-2,80 XX 2,80-2,90 X 2,90-3,00 X 3,00-3,10 X 3,10-3,20 X 3,20-3,30 3,30-3,40 3,40-3,50
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 * *** ******** ****************** ***************************** ************************************** ************************************** ***************************** ****************** ******** *** *