Statistics/Normal distribution

Statistics/Normal distribution
You are encouraged to solve this task according to the task description, using any language you may know.

The Normal (or Gaussian) distribution is a frequently used distribution in statistics. While most programming languages provide a uniformly distributed random number generator, one can derive normally distributed random numbers from a uniform generator.

The task

1. 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.
2. 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

/*
 * 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;
}

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#

Library: Math.Net

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);
    }
}

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

#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 ;
}

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.

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;
}

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

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)

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: =

Factor

Works with: Factor version 0.99 2020-01-23

USING: assocs formatting kernel math math.functions
math.statistics random sequences sorting ;

2,000,000 [ 0 1 normal-random-float ] replicate   ! make data set
dup [ mean ] [ population-std ] bi                ! calculate and show
"Mean: %f\nStdev: %f\n\n" printf                  ! mean and stddev

[ 10 * floor 10 / ] map                   ! map data to buckets
histogram >alist [ first ] sort-with      ! create histogram sorted by bucket (key)
dup values supremum                       ! find maximum count
[
    [ /f 100 * >integer ] keepd             ! how big should this histogram bar be?
    [ [ CHAR: * ] "" replicate-as ] dip     ! make the bar
    "% 5.2f: %s   %d\n" printf              ! print a line of the histogram
] curry assoc-each

Output:

Mean: 0.000798
Stdev: 1.000549

-4.90:    2
-4.80:    1
-4.70:    1
-4.60:    3
-4.50:    3
-4.40:    6
-4.30:    15
-4.20:    13
-4.10:    16
-4.00:    42
-3.90:    62
-3.80:    68
-3.70:    98
-3.60:    145
-3.50:    205
-3.40:    311
-3.30:    379
-3.20:    580
-3.10:    739
-3.00: *   1002
-2.90: *   1349
-2.80: **   1893
-2.70: ***   2499
-2.60: ****   3211
-2.50: *****   4035
-2.40: ******   5141
-2.30: *******   6392
-2.20: *********   7869
-2.10: ************   9780
-2.00: **************   11787
-1.90: ******************   14483
-1.80: *********************   17183
-1.70: *************************   20387
-1.60: ******************************   24049
-1.50: **********************************   27555
-1.40: ****************************************   32153
-1.30: *********************************************   36707
-1.20: ***************************************************   40921
-1.10: *********************************************************   45928
-1.00: ***************************************************************   50707
-0.90: *********************************************************************   55697
-0.80: ***************************************************************************   60377
-0.70: ********************************************************************************   64358
-0.60: ************************************************************************************   67928
-0.50: *****************************************************************************************   71911
-0.40: *********************************************************************************************   75054
-0.30: ************************************************************************************************   77073
-0.20: **************************************************************************************************   78768
-0.10: ***************************************************************************************************   79732
 0.00: ****************************************************************************************************   79952
 0.10: ***************************************************************************************************   79412
 0.20: ************************************************************************************************   77511
 0.30: *********************************************************************************************   74487
 0.40: ******************************************************************************************   72250
 0.50: **************************************************************************************   68789
 0.60: ********************************************************************************   64408
 0.70: ***************************************************************************   60122
 0.80: *********************************************************************   55619
 0.90: ***************************************************************   50869
 1.00: *********************************************************   45883
 1.10: ****************************************************   41586
 1.20: **********************************************   37145
 1.30: ***************************************   31715
 1.40: **********************************   27779
 1.50: ******************************   24270
 1.60: *************************   20516
 1.70: *********************   17221
 1.80: *****************   14344
 1.90: **************   11789
 2.00: ************   9796
 2.10: *********   7922
 2.20: *******   6331
 2.30: ******   5138
 2.40: *****   4044
 2.50: ***   3065
 2.60: **   2397
 2.70: **   1846
 2.80: *   1462
 2.90: *   1001
 3.00:    765
 3.10:    587
 3.20:    393
 3.30:    299
 3.40:    197
 3.50:    132
 3.60:    100
 3.70:    74
 3.80:    59
 3.90:    32
 4.00:    29
 4.10:    12
 4.20:    15
 4.30:    6
 4.40:    3
 4.50:    4
 4.60:    3
 4.70:    2
 4.80:    1

Fortran

Works with: Fortran version 95 and later

Using the Marsaglia polar method

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

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

' 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

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.

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))
    }
}

Output:

mean: -0.0034970888831523488
stddev: 1.0040682925006286

*
****
*********
***************
*******************
******************
**************
*********
****
*

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

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

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.)

Example Usage

   DataSet=: rnorm01 1e5
   (mean , stddev) DataSet
0.000781667 1.00154
   require 'plot'
   plot (5 %~ i: 25) ([;histogram) DataSet

Java

Translation of: D

Works with: Java version 8

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());
    }
}

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: **

jq

Adapted from Wren

Works with: jq

Works with gojq, the Go implementation of jq (*)

Since jq does not have a built-in PRNG, this entry uses an external source for entropy. For the sake of illustration, we will use /dev/urandom as follows:

  cat /dev/urandom | tr -cd '0-9' | fold -w 10 |
     jq -nRr -f normal-distribution.jq

To save space, the function that generates the sample does not retain the observations, and for simplicity, computes the sum of squared observations on the assumption that overflow will not be an an issue, which is reasonable as jq arithmetic uses IEEE 754 64-bit numbers.

(*) gojq requires an enormous amount of memory to complete the task for N=100,000, and takes about 20 times longer.

Preliminaries

# Pretty print a number to facilitate alignment of the decimal point.
# Input: a number without an exponent
# Output: a string holding the reformatted number so that there are at least `left` characters
# to the left of the decimal point, and exactly `right` characters to its right.
# Spaces are used for padding on the left, and zeros for padding on the right.
# No left-truncation occurs, so `left` can be specified as 0 to prevent left-padding.
def pp(left; right):
  def lpad: tostring | (left - length) as $l | (" " * $l)[:$l] + .;
  def rpad:
    if (right > length) then . + ((right - length) * "0")
    else .[:right]
    end;
  tostring as $s
  | $s
  | index(".") as $ix
    | ((if $ix then $s[0:$ix] else $s end) | lpad) + "." +
       (if $ix then $s[$ix+1:] | .[:right] else "" end | rpad) ;

def sigma( stream ): reduce stream as $x (0; . + $x);

# Input: {n, sum, ss}
# Output: augmented object with {mean, variance}
def sample_mean_and_variance:
  .mean = (.sum/.n)
  | .variance = ((.ss / .n) - .mean*.mean);

The Task

# Task parameters
def parameters: {
    N:         100000,
    NUM_BINS:  12,
    HIST_CHAR: "■",
    HIST_CHAR_ALT: "-",
    HIST_CHAR_SIZE: null,  # null means compute dynamically
    binSize:   0.1,
    mu:        0.5,
    sigma:     0.25 }
    | .bins = [range(0; .NUM_BINS) | 0] ;

# input: an array of two iid rvs on [0,1]
# output: [z0, z1] as per the Box-Muller method -- see
# https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
def normal(mu; sigma):
  def pi: (1|atan) * 4;
  . as [$u1, $u2]
  | pi as $pi
  | (sigma * ((-2 * ($u1|log))|sqrt)) as $mag
  | [ $mag * ((2 * $pi * $u2)|cos) + mu,
      $mag * ((2 * $pi * $u2)|sin) + mu ] ;

# Generate a random sample as specified by ., the task object (see `parameters`).
# Output: updated task object with sample statistics and .bins for creating a histogram.
# Each call to `input` should yield a string of random decimal digits 
# such that the ensemble of ("0." + input | tonumber) can be considered to be iid rv on [0,1].
def generate:
  # uniformly distributed random variable on [0,1]:
  def udrv: "0." + input | tonumber;
  # Maybe compute the bucket size:
  (.HIST_CHAR_SIZE = (.HIST_CHAR_SIZE // (.N / (.NUM_BINS * 20) | ceil))) as $p
  | reduce range(0; $p.N/2) as $i ($p;
      ([udrv, udrv] | normal($p.mu; $p.sigma)) as $rns
      | reduce (0,1) as $j (.;
          $rns[$j] as $rn
	  | .n += 1
	  | .sum += $rn
	  | .ss  += ($rn*$rn)
          | (if $rn < 0 then 0
 	     elif $rn >= 1 then ($p.NUM_BINS - 1)
             else  ($rn/.binSize)|floor + 1
	     end ) as $bn
          | .bins[$bn] += 1
          # to retain the observations: .samples[$i*2 + $j] = $rn
	  )) ;

# Input: an object with
# {NUM_BINS, HIST_CHAR_SIZE, HIST_CHAR, HIST_CHAR_ALT, binSize, bins}
# Output: a stream of strings
def histogram:
  def tidy: pp(2;1);
  range(0; .NUM_BINS) as $i
  | ((.bins[$i] / .HIST_CHAR_SIZE)|floor) as $bs
  | (if $i == 0 or $i == .NUM_BINS -1
     then .HIST_CHAR_ALT else .HIST_CHAR end) as $char
  | (if $bs == 0 then "" else $char * $bs end) as $hist
  | if $i == 0
    then " -∞  ..< 0.0 \($hist)"    #   .bins[0]
    elif ($i < .NUM_BINS - 1)
    then "\(.binSize * ($i-1) | tidy) ..<\(.binSize * $i|tidy) \($hist)"  # .bins[$i]]
    else " 1.0 ..  +∞  \($hist)"    #   .bins[.NUM_BINS - 1]
    end;

def task:
  parameters
  | generate
  | sample_mean_and_variance
  | (if .HIST_CHAR_SIZE == 1 then "" else "s" end) as $plural
  | "Summary statistics for \(.N) observations from N(\(.mu), \(.sigma)):",
     "    mean:              \(.mean | pp(2;4))",
     "    variance:          \(.variance | pp(2;4))",
     "    unadjusted stddev: \(.variance | sqrt | pp(2;4))",
     "    Range               Number of observations (each \(.HIST_CHAR) represents \(.HIST_CHAR_SIZE) observation\($plural))",
     histogram ;

task

Output:

Summary statistics for 100000 observations from N(0.5, 0.25):
    mean:               0.5001
    variance:           0.0622
    unadjusted stddev:  0.2495
    Range               Number of observations (each ■ represents 417 observations)
 -∞  ..< 0.0 -----
 0.0 ..< 0.1 ■■■■■■■■
 0.1 ..< 0.2 ■■■■■■■■■■■■■■
 0.2 ..< 0.3 ■■■■■■■■■■■■■■■■■■■■■■■
 0.3 ..< 0.4 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
 0.4 ..< 0.5 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
 0.5 ..< 0.6 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
 0.6 ..< 0.7 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
 0.7 ..< 0.8 ■■■■■■■■■■■■■■■■■■■■■■■
 0.8 ..< 0.9 ■■■■■■■■■■■■■■
 0.9 ..< 1.0 ■■■■■■■■
 1.0 ..  +∞  ------

Julia

Julia has the builtin package "Distributions" to generate random numbers from a standard distribution (Normal, Chisq etc.).

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)

Output:

N = 1000
μ = 0.02	σ = 0.97
range = (-2.76, 3.42)

Kotlin

Translation of: FreeBASIC

// 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)
}

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

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'
^}

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

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

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.

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)

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:

with(Statistics):
n := 100000:
X := Sample( Normal(0,1), n );
Mean( X );
StandardDeviation( X );
Histogram( X );

Mathematica/Wolfram Language

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

MATLAB / Octave

  N = 100000;	
  x = randn(N,1);
  mean(x)
  std(x) 
  [nn,xx] = hist(x,100);
  bar(xx,nn);

Nim

In module “random”, Nim provides two procedures named gauss to generate random values following normal distribution and following Gauss distribution with given mean and standard deviation.

Here is a way to generate random values following normal distribution from random values following uniform distribution. It uses the Basic form of the Box-Muller transform.

import math, random, sequtils, stats, strformat, strutils

proc drawHistogram(ns: seq[float]) =

  # Distribute values in bins.
  const NBins = 50
  var minval = min(ns)
  var maxval = max(ns)
  var h = newSeq[int](NBins + 1)
  for n in ns:
    let pos = ((n - minval) * NBins / (maxval - minval)).toInt
    inc h[pos]

  # Eliminate extremes values.
  const MaxWidth = 50
  let mx = max(h)
  var first = 0
  while (h[first] / mx * MaxWidth).toInt == 0: inc first
  var last = h.high
  while (h[last] / mx * MaxWidth).toInt == 0: dec last

  # Draw the histogram.
  echo ""
  for n in first..last:
    echo repeat('+', (h[n] / mx * MaxWidth).toInt)
  echo ""


const N = 100_000

randomize()

let u1, u2 = newSeqWith(N, rand(1.0))

var z = newSeq[float](N)
for i in 0..<N:
  z[i] = sqrt(-2 * ln(u1[i])) * cos(2 * PI * u2[i])

echo &"μ = {z.mean:.12f}   σ = {z.standardDeviation:.12f}"
z.drawHistogram()

Output:

μ = -0.001105836229   σ = 0.999906544722

+
+
++
+++
+++++
+++++++
+++++++++
++++++++++++
++++++++++++++++
+++++++++++++++++++++
++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++
+++++++++++++++++++++++++
++++++++++++++++++++
++++++++++++++++
++++++++++++
+++++++++
++++++
+++++
+++
++
+
+

PARI/GP

Works with: PARI/GP version 2.4.3 and above

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)

For versions before 2.4.3, define

rreal()={
  my(pr=32*ceil(default(realprecision)*log(10)/log(4294967296))); \\ Current precision
  random(2^pr)*1.>>pr
};

and use rreal() in place of random(1.).

A PARI implementation:

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;
}

Use mpsincos and caching to generate two values at nearly the same cost.

Pascal

Works with: free 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

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.

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

Translation of: Raku

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";
    }

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

Translation of: Liberty_BASIC

with javascript_semantics
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
        integer bdx = floor((dat[i]-mn)/range*nBins)+1
        bins[bdx] += 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)

Output:

100000 data terms used.
Largest term was 4.30779 & smallest was -4.11902
Mean = -0.00252597
Stddev = 1.00067

#
##
####
######
##########
#############
##################
########################
#################################
########################################
####################################################
#############################################################
######################################################################
###############################################################################
#######################################################################################
###############################################################################################
#################################################################################################
#####################################################################################################
###################################################################################################
################################################################################################
########################################################################################
###############################################################################
#######################################################################
############################################################
#################################################
#######################################
##############################
#########################
################
############
#########
######
####
##
#

Translation of: Lua

with javascript_semantics
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)

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

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

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.

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)

Output:

Sample mean = 49.9822; Stddev = 4.00938; max = 66.8091; min = 33.5283 for 100000 values

R

Generate normal random numbers from uniform random numbers, using the Box-Muller transform. Both x and y are normally distributed, and independent.

n <- 100000
u <- sqrt(-2*log(runif(n)))
v <- 2*pi*runif(n)
x <- u*cos(v)
y <- v*sin(v)
hist(x)

R can generate random normal distributed numbers using the rnorm function:

n <- 100000
x <- rnorm(n, mean=0, sd=1)
mean(x)
sd(x)
hist(x)

Racket

This shows how one would generate samples from a normal distribution, compute statistics and plot a histogram.

#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")

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
(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))])))))))

Raku

(formerly Perl 6)

Works with: Rakudo version 2018.03

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];
    }
}

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.

/*REXX program generates  10,000  normally distributed numbers  (Gaussian distribution).*/
numeric digits 20                                /*use twenty decimal digs 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            /*the sum,  and  the sum of squares.   */
        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 = '   fmt(n  )
say '              minimum = '   fmt(mn )
say '              maximum = '   fmt(mx )
say '      arithmetic mean = '   fmt(s/n)
say '   standard deviation = '   fmt(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 (usable) 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, !.?)                   /*find the minimum.                    */
            ?mx= max(?mx, !.?)                   /*  "   "  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. */
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/
fmt:  parse arg @; return left('', (@>=0) + 2 * datatype(@, 'W'))@  /*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.*/
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/
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;   z= 1;   _= 1
      x= x*x;  p= z;      do k=2  by 2; _= -_ * x / (k*(k-1));   z= z + _;  if z=p  then leave;  p= z; end;    return z
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/
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

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   ³/₄   size.

number of data points =     10000
              minimum =  -3.8181072371544448250
              maximum =   3.5445917138265268562
      arithmetic mean =  -0.01406470979976873427
   standard deviation =   0.99486092191249231518
─
─
───
────
─────
─────
────────
───────────
──────────────
─────────────────────
──────────────────────
──────────────────────────────────
────────────────────────────────────────
───────────────────────────────────────────────
─────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────────────────────────────────────────────
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────
─────────────────────────────────────────────────
─────────────────────────────────
──────────────────────────────────
───────────────────────
──────────────────────
──────────────────
───────────
──────────
──────
───
────
──
─

Ruby

Works with: Ruby version 2.7

The Implementation

# Class to implement a Normal distribution, generated from a Uniform distribution.
# Uses the Marsaglia polar method.

class NormalFromUniform
  # Initialize an instance.
  def initialize()
    @next = nil
  end
  # Generate and return the next Normal distribution value.
  def rand()
    if @next
      retval, @next = @next, nil
      return retval
    else
      u = v = s = nil
      loop do
        u = Random.rand(-1.0..1.0)
        v = Random.rand(-1.0..1.0)
        s = u**2 + v**2
        break if (s > 0.0) && (s <= 1.0)
      end
      f = Math.sqrt(-2.0 * Math.log(s) / s)
      @next = v * f
      return u * f
    end
  end
end

The Task

Library: enumerable-statistics

require('enumerable/statistics')

def show_stats_and_histogram(data, bins)
  puts("size = #{data.length}  mean = #{data.mean()}  stddev = #{data.stdev()}")
  hist = data.histogram(bins)
  scale = 100.0 / hist.weights.max
  inx_beg = nil
  inx_end = nil
  hist.weights.length.times do |inx|
    histstars = (0.5 + (scale * hist.weights[inx])).to_i
    inx_beg = inx if !inx_beg && (histstars > 0)
    inx_end = inx if (histstars > 0)
  end
  (inx_beg..inx_end).each do |inx|
    bincenter = 0.5 * (hist.edges[inx] + hist.edges[inx + 1])
    histstars = (0.5 + (scale * hist.weights[inx])).to_i
    puts('%6.2f: %s' % [bincenter, '*' * histstars])
  end
end

puts
puts('Uniform random number generator:')
show_stats_and_histogram(1000000.times.map { Random.rand(-1.0..1.0) }, 20)

puts
puts('Normal random numbers using the Marsaglia polar method:')
gen_normal = NormalFromUniform.new
show_stats_and_histogram(100.times.map { gen_normal.rand }, 40)
show_stats_and_histogram(10000.times.map { gen_normal.rand }, 60)
show_stats_and_histogram(1000000.times.map { gen_normal.rand }, 120)

Output:

Uniform random number generator:
size = 1000000  mean = 0.0001483724103528628  stddev = 0.5773085740222473
 -0.95: ****************************************************************************************************
 -0.85: ****************************************************************************************************
 -0.75: ***************************************************************************************************
 -0.65: ***************************************************************************************************
 -0.55: ***************************************************************************************************
 -0.45: ****************************************************************************************************
 -0.35: ****************************************************************************************************
 -0.25: ****************************************************************************************************
 -0.15: ****************************************************************************************************
 -0.05: ***************************************************************************************************
  0.05: ****************************************************************************************************
  0.15: ****************************************************************************************************
  0.25: ****************************************************************************************************
  0.35: ***************************************************************************************************
  0.45: ***************************************************************************************************
  0.55: ****************************************************************************************************
  0.65: ****************************************************************************************************
  0.75: ****************************************************************************************************
  0.85: ****************************************************************************************************
  0.95: ***************************************************************************************************

Normal random numbers using the Marsaglia polar method:
size = 100  mean = 0.1611961166227389  stddev = 0.9827581078952005
 -2.30: **********
 -2.10: 
 -1.90: **********
 -1.70: ********************
 -1.50: 
 -1.30: **********
 -1.10: **********************************************************************
 -0.90: ************************************************************
 -0.70: **********************************************************************
 -0.50: ************************************************************
 -0.30: ********************************************************************************
 -0.10: ********************
  0.10: ********************************************************************************
  0.30: ****************************************************************************************************
  0.50: ******************************************************************************************
  0.70: ********************************************************************************
  0.90: ************************************************************
  1.10: ******************************
  1.30: **************************************************
  1.50: 
  1.70: ********************
  1.90: **************************************************
  2.10: **********
  2.30: ********************
size = 10000  mean = -0.004863294071004369  stddev = 0.9984395022628252
 -3.30: *
 -3.10: *
 -2.90: **
 -2.70: **
 -2.50: ****
 -2.30: *********
 -2.10: ***********
 -1.90: ****************
 -1.70: ***********************
 -1.50: *****************************
 -1.30: *****************************************
 -1.10: *************************************************
 -0.90: ******************************************************************
 -0.70: ******************************************************************************
 -0.50: ***************************************************************************************
 -0.30: *********************************************************************************************
 -0.10: ***********************************************************************************************
  0.10: ****************************************************************************************************
  0.30: **************************************************************************************
  0.50: ************************************************************************************
  0.70: *******************************************************************************
  0.90: ****************************************************************
  1.10: ***************************************************
  1.30: ********************************************
  1.50: *******************************
  1.70: **********************
  1.90: ****************
  2.10: **********
  2.30: ******
  2.50: *****
  2.70: **
  2.90: *
  3.10: *
size = 1000000  mean = 0.0007049206911295587  stddev = 1.0000356747411552
 -3.15: *
 -3.05: *
 -2.95: *
 -2.85: **
 -2.75: **
 -2.65: ***
 -2.55: ****
 -2.45: *****
 -2.35: ******
 -2.25: ********
 -2.15: **********
 -2.05: ************
 -1.95: ***************
 -1.85: ******************
 -1.75: *********************
 -1.65: *************************
 -1.55: ******************************
 -1.45: ***********************************
 -1.35: ****************************************
 -1.25: *********************************************
 -1.15: ***************************************************
 -1.05: *********************************************************
 -0.95: ***************************************************************
 -0.85: *********************************************************************
 -0.75: **************************************************************************
 -0.65: *********************************************************************************
 -0.55: *************************************************************************************
 -0.45: *****************************************************************************************
 -0.35: ********************************************************************************************
 -0.25: ************************************************************************************************
 -0.15: **************************************************************************************************
 -0.05: ****************************************************************************************************
  0.05: ***************************************************************************************************
  0.15: **************************************************************************************************
  0.25: ************************************************************************************************
  0.35: *********************************************************************************************
  0.45: ******************************************************************************************
  0.55: *************************************************************************************
  0.65: ********************************************************************************
  0.75: **************************************************************************
  0.85: *********************************************************************
  0.95: ***************************************************************
  1.05: *********************************************************
  1.15: ****************************************************
  1.25: **********************************************
  1.35: ****************************************
  1.45: **********************************
  1.55: ******************************
  1.65: *************************
  1.75: **********************
  1.85: ******************
  1.95: ***************
  2.05: ************
  2.15: **********
  2.25: ********
  2.35: ******
  2.45: *****
  2.55: ****
  2.65: ***
  2.75: **
  2.85: **
  2.95: *
  3.05: *
  3.15: *

Run BASIC

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

Output:

Samples   :100000
Largest   :4.61187177
Smallest  :-4.21695424
Range     :8.82882601
Mean      :-9.25042513e-4
Stand Dev :1.00680067

=
==
===
=====
========
=============
=================
=======================
==============================
=======================================
===============================================
=========================================================
===================================================================
===========================================================================
===================================================================================
=======================================================================================
==========================================================================================
========================================================================================
======================================================================================
=================================================================================
============================================================================
==================================================================
========================================================
==============================================
=====================================
============================
=====================
===============
==========
=======
=====
===
=
=

Rust

Library: math

Library: rand

Library: rand_distr

//! Rust rosetta example for normal distribution
use math::{histogram::Histogram, traits::ToIterator};
use rand;
use rand_distr::{Distribution, Normal};

/// Returns the mean of the provided samples
///
/// ## Arguments
/// * data -- reference to float32 array
fn mean(data: &[f32]) -> Option<f32> {
    let sum: f32 = data.iter().sum();
    Some(sum / data.len() as f32)
}

/// Returns standard deviation of the provided samples
///
/// ## Arguments
/// * data -- reference to float32 array
fn standard_deviation(data: &[f32]) -> Option<f32> {
    let mean = mean(data).expect("invalid mean");
    let sum = data.iter().fold(0.0, |acc, &x| acc + (x - mean).powi(2));
    Some((sum / data.len() as f32).sqrt())
}

/// Prints a histogram in the shell
///
/// ## Arguments
/// * data -- reference to float32 array
/// * maxwidth -- the maxwidth of the histogram in # of characters
/// * bincount -- number of bins in the histogram
/// * ch -- character used to plot the graph
fn print_histogram(data: &[f32], maxwidth: usize, bincount: usize, ch: char) {
    let min_val = data.iter().cloned().fold(f32::NAN, f32::min);
    let max_val = data.iter().cloned().fold(f32::NAN, f32::max);
    let histogram = Histogram::new(Some(&data.to_vec()), bincount, min_val, max_val).unwrap();
    let max_bin_value = histogram.get_counters().iter().max().unwrap();
    println!();
    for x in histogram.to_iter() {
        let (bin_min, bin_max, freq) = x;
        let bar_width = (((freq as f64) / (*max_bin_value as f64)) * (maxwidth as f64)) as u32;
        let bar_as_string = (1..bar_width).fold(String::new(), |b, _| b + &ch.to_string());
        println!(
            "({:>6},{:>6}) |{} {:.2}%",
            format!("{:.2}", bin_min),
            format!("{:.2}", bin_max),
            bar_as_string,
            (freq as f64) * 100.0 / (data.len() as f64)
        );
    }
    println!();
}

/// Runs the demo to generate normal distribution of three different sample sizes
fn main() {
    let expected_mean: f32 = 0.0;
    let expected_std_deviation: f32 = 4.0;
    let normal = Normal::new(expected_mean, expected_std_deviation).unwrap();

    let mut rng = rand::thread_rng();
    for &number_of_samples in &[1000, 10_000, 1_000_000] {
        let data: Vec<f32> = normal
            .sample_iter(&mut rng)
            .take(number_of_samples)
            .collect();
        println!("Statistics for sample size {}:", number_of_samples);
        println!("\tMean: {:?}", mean(&data).expect("invalid mean"));
        println!(
            "\tStandard deviation: {:?}",
            standard_deviation(&data).expect("invalid standard deviation")
        );
        print_histogram(&data, 80, 40, '-');
    }
}

Output:

Statistics for sample size 1000:
	Mean: -0.11356559
	Standard deviation: 4.0244555

(-13.81,-13.12) | 0.10%
(-13.12,-12.44) | 0.00%
(-12.44,-11.75) | 0.10%
(-11.75,-11.06) | 0.20%
(-11.06,-10.38) |- 0.30%
(-10.38, -9.69) | 0.10%
( -9.69, -9.01) |--- 0.50%
( -9.01, -8.32) |---- 0.60%
( -8.32, -7.64) |------ 0.80%
( -7.64, -6.95) |-------------- 1.60%
( -6.95, -6.27) |----------------- 1.90%
( -6.27, -5.58) |------------------------ 2.60%
( -5.58, -4.90) |----------------------- 2.50%
( -4.90, -4.21) |---------------------------------------- 4.20%
( -4.21, -3.53) |------------------------------------- 3.90%
( -3.53, -2.84) |------------------------------------------------- 5.10%
( -2.84, -2.15) |---------------------------------------------- 4.80%
( -2.15, -1.47) |------------------------------------------------------------------------ 7.40%
( -1.47, -0.78) |---------------------------------------------------------- 6.00%
( -0.78, -0.10) |----------------------------------------------------------------------- 7.30%
( -0.10,  0.59) |------------------------------------------------------------------------------- 8.10%
(  0.59,  1.27) |----------------------------------------------------------------------- 7.30%
(  1.27,  1.96) |------------------------------------------------- 5.10%
(  1.96,  2.64) |------------------------------------------------------------ 6.20%
(  2.64,  3.33) |----------------------------------------- 4.30%
(  3.33,  4.01) |----------------------------- 3.10%
(  4.01,  4.70) |------------------------------------- 3.90%
(  4.70,  5.39) |-------------------------- 2.80%
(  5.39,  6.07) |---------------------- 2.40%
(  6.07,  6.76) |---------------- 1.80%
(  6.76,  7.44) |---------------- 1.80%
(  7.44,  8.13) |--------- 1.10%
(  8.13,  8.81) |---------- 1.20%
(  8.81,  9.50) | 0.20%
(  9.50, 10.18) | 0.00%
( 10.18, 10.87) | 0.10%
( 10.87, 11.55) |- 0.30%
( 11.55, 12.24) | 0.10%
( 12.24, 12.92) | 0.10%
( 12.92, 13.61) | 0.10%

Statistics for sample size 10000:
	Mean: 0.02012564
	Standard deviation: 3.9697735

(-15.80,-14.99) | 0.02%
(-14.99,-14.18) | 0.04%
(-14.18,-13.37) | 0.04%
(-13.37,-12.56) | 0.04%
(-12.56,-11.75) | 0.09%
(-11.75,-10.94) | 0.08%
(-10.94,-10.13) |- 0.25%
(-10.13, -9.32) |--- 0.42%
( -9.32, -8.51) |----- 0.67%
( -8.51, -7.70) |--------- 1.10%
( -7.70, -6.89) |------------- 1.48%
( -6.89, -6.08) |------------------ 1.98%
( -6.08, -5.27) |-------------------------- 2.82%
( -5.27, -4.45) |------------------------------------ 3.80%
( -4.45, -3.64) |--------------------------------------------- 4.66%
( -3.64, -2.83) |------------------------------------------------------- 5.72%
( -2.83, -2.02) |------------------------------------------------------------------ 6.85%
( -2.02, -1.21) |---------------------------------------------------------------------------- 7.80%
( -1.21, -0.40) |---------------------------------------------------------------------------- 7.79%
( -0.40,  0.41) |------------------------------------------------------------------------------- 8.09%
(  0.41,  1.22) |----------------------------------------------------------------------------- 7.89%
(  1.22,  2.03) |--------------------------------------------------------------------------- 7.76%
(  2.03,  2.84) |-------------------------------------------------------------------- 7.06%
(  2.84,  3.65) |------------------------------------------------------- 5.74%
(  3.65,  4.46) |-------------------------------------------- 4.64%
(  4.46,  5.27) |-------------------------------------- 4.00%
(  5.27,  6.08) |---------------------------- 3.03%
(  6.08,  6.89) |------------------- 2.07%
(  6.89,  7.71) |-------------- 1.54%
(  7.71,  8.52) |-------- 0.97%
(  8.52,  9.33) |----- 0.61%
(  9.33, 10.14) |-- 0.36%
( 10.14, 10.95) |- 0.25%
( 10.95, 11.76) | 0.19%
( 11.76, 12.57) | 0.08%
( 12.57, 13.38) | 0.02%
( 13.38, 14.19) | 0.01%
( 14.19, 15.00) | 0.03%
( 15.00, 15.81) | 0.00%
( 15.81, 16.62) | 0.01%

Statistics for sample size 1000000:
	Mean: -0.004743685
	Standard deviation: 4.0006065

(-20.00,-19.02) | 0.00%
(-19.02,-18.04) | 0.00%
(-18.04,-17.06) | 0.00%
(-17.06,-16.07) | 0.00%
(-16.07,-15.09) | 0.00%
(-15.09,-14.11) | 0.01%
(-14.11,-13.13) | 0.03%
(-13.13,-12.15) | 0.06%
(-12.15,-11.16) | 0.14%
(-11.16,-10.18) |- 0.28%
(-10.18, -9.20) |--- 0.53%
( -9.20, -8.22) |------ 0.95%
( -8.22, -7.24) |----------- 1.51%
( -7.24, -6.25) |------------------ 2.40%
( -6.25, -5.27) |--------------------------- 3.48%
( -5.27, -4.29) |-------------------------------------- 4.82%
( -4.29, -3.31) |-------------------------------------------------- 6.27%
( -3.31, -2.32) |------------------------------------------------------------- 7.62%
( -2.32, -1.34) |----------------------------------------------------------------------- 8.77%
( -1.34, -0.36) |----------------------------------------------------------------------------- 9.58%
( -0.36,  0.62) |------------------------------------------------------------------------------- 9.74%
(  0.62,  1.60) |---------------------------------------------------------------------------- 9.39%
(  1.60,  2.59) |-------------------------------------------------------------------- 8.49%
(  2.59,  3.57) |---------------------------------------------------------- 7.30%
(  3.57,  4.55) |----------------------------------------------- 5.86%
(  4.55,  5.53) |----------------------------------- 4.45%
(  5.53,  6.51) |------------------------ 3.16%
(  6.51,  7.50) |---------------- 2.09%
(  7.50,  8.48) |--------- 1.34%
(  8.48,  9.46) |----- 0.81%
(  9.46, 10.44) |-- 0.46%
( 10.44, 11.42) | 0.23%
( 11.42, 12.41) | 0.11%
( 12.41, 13.39) | 0.06%
( 13.39, 14.37) | 0.02%
( 14.37, 15.35) | 0.01%
( 15.35, 16.34) | 0.00%
( 16.34, 17.32) | 0.00%
( 17.32, 18.30) | 0.00%
( 18.30, 19.28) | 0.00%

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
*/

Sidef

Translation of: Raku

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])
}

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.

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)

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

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

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

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  

Wren

Library: Wren-fmt

Library: Wren-math

import "random" for Random
import "/fmt" for Fmt
import "/math" for Nums

var rgen = Random.new()

// Box-Muller method from Wikipedia
var normal  = Fn.new { |mu, sigma|
    var u1  = rgen.float()
    var u2  = rgen.float()
    var mag = sigma * (-2 * u1.log).sqrt
    var z0  = mag * (2 * Num.pi * u2).cos + mu
    var z1  = mag * (2 * Num.pi * u2).sin + mu
    return [z0, z1]
}

var N = 100000
var NUM_BINS = 12
var HIST_CHAR = "■"
var HIST_CHAR_SIZE = 250
var bins = List.filled(NUM_BINS, 0)
var binSize = 0.1
var samples = List.filled(N, 0)
var mu = 0.5
var sigma = 0.25
for (i in 0...N/2) {
    var rns = normal.call(mu, sigma)
    for (j in 0..1) {
        var rn = rns[j]
        var bn
        if (rn < 0) {
            bn = 0
        } else if (rn >= 1) {
            bn = 11
        } else {
            bn = (rn/binSize).floor + 1
        }
        bins[bn] = bins[bn] + 1
        samples[i*2 + j] = rn
    }
}

Fmt.print("Normal distribution with mean $0.2f and S/D $0.2f for $,d samples:\n", mu, sigma, N)
System.print("    Range           Number of samples within that range")
for (i in 0...NUM_BINS) {
    var hist = HIST_CHAR * (bins[i] / HIST_CHAR_SIZE).round
    if (i == 0) {
        Fmt.print("  -∞ ..< 0.00  $s $,d", hist, bins[0])
    } else if (i < NUM_BINS - 1) {
        Fmt.print("$4.2f ..< $4.2f  $s $,d", binSize * (i-1), binSize * i, hist, bins[i])
    } else {
        Fmt.print("1.00 ... +∞    $s $,d", hist, bins[NUM_BINS - 1])
    }
}
Fmt.print("\nActual mean for these samples : $0.5f", Nums.mean(samples))
Fmt.print("Actual S/D  for these samples : $0.5f", Nums.stdDev(samples))

Output:

Specimen run:

Normal distribution with mean 0.50 and S/D 0.25 for 100,000 samples:

    Range           Number of samples within that range
  -∞ ..< 0.00  ■■■■■■■■■ 2,243
0.00 ..< 0.10  ■■■■■■■■■■■■■ 3,250
0.10 ..< 0.20  ■■■■■■■■■■■■■■■■■■■■■■■■ 5,977
0.20 ..< 0.30  ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 9,723
0.30 ..< 0.40  ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 13,104
0.40 ..< 0.50  ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 15,601
0.50 ..< 0.60  ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 15,469
0.60 ..< 0.70  ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 13,334
0.70 ..< 0.80  ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 9,659
0.80 ..< 0.90  ■■■■■■■■■■■■■■■■■■■■■■■■ 6,119
0.90 ..< 1.00  ■■■■■■■■■■■■■ 3,173
1.00 ... +∞    ■■■■■■■■■ 2,348

Actual mean for these samples : 0.50099
Actual S/D  for these samples : 0.25051

zkl

Translation of: Go

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;
};

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.

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)) }

Output:

Samples: 100,000
Mean:    0.0005999
Stddev:  1.003
*
***
********
******************
*****************************
**************************************
**************************************
*****************************
******************
********
***
*