Random numbers

Generate a collection filled with 1000 normally distributed random (or pseudo-random) numbers with a mean of 1.0 and a standard deviation of 0.5

Task

Many libraries only generate uniformly distributed random numbers. If so, you may use one of these algorithms.

Related task

Standard deviation

ALGOL 68

Translation of: C

Works with: ALGOL 68 version Revision 1 - no extensions to language used

Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny

<lang algol68>PROC random normal = REAL: # normal distribution, centered on 0, std dev 1 # (

 sqrt(-2*log(random)) * cos(2*pi*random)

);

test:(

 [1000]REAL rands;
 FOR i TO UPB rands DO
   rands[i] := 1 + random normal/2
 OD;
 INT limit=10;
 printf(($"("n(limit-1)(-d.6d",")-d.5d" ... )"$, rands[:limit]))

)</lang>

Output:

( 0.693461, 0.948424, 0.482261, 1.045939, 0.890818, 1.467935, 0.604153, 0.804811, 0.690227, 0.83462 ... )

Arturo

<lang rebol>rnd: function []-> (random 0 10000)//10000

rands: map 1..1000 'x [

   1 + (sqrt neg 2 * ln rnd) * (cos 2 * pi * rnd)

]

print rands</lang>

Output:

0.6219537961087694 1.279396486161406 1.619019280815647 2.119538294228789 -0.1598383851981044 2.67797211803156 0.9304703037226587 1.629254364659528 -0.4171704717398712 0.9082342931486092 -0.5929704390625219 2.117000897984871 -0.1981633787460266 0.01132471973856153 2.102359263212924 0.2408823232884222 2.046195035792376 0.6374831627030295 0.000839808324124558 1.117061838266626 0.7413355299469649 0.4485598815755762 2.999434800016997 2.560580541932842 1.703197984879731 2.889159248353575 -1.800157205708138 1.756810020187321 0.7136708180852145 0.5929151678321705 0.332519993787973 2.660212054362758 0.5835660585480075 0.8527946892567934 1.640573993747053 0.09471843345263908 1.051402997891346 1.116149156137905 -0.7400139019343499 1.782572831979232 2.531779039786426 0.5240268064639871 0.07099232630526586 -0.854892656700071 1.54381929430469 -0.4416899008614745 0.4274356035015117 0.7350027625573482 2.153583935076981 1.461215281535983 -1.041723064151266 2.338060763553139 -0.1492967916030414 0.3799517724040202 0.4577924541353815 0.673317567666373 -2.27731583876462 1.28480355806061 -0.6925811023772748 -0.2642224122781984 0.6590513830891744 2.55537133425143 1.67933335469247 0.8659013395355968 -1.211026941441126 0.9524579534222226 -0.1931750631835656 -0.5119479869237693 0.1814749003063878 3.03139579963414 ...

Avail

<lang Avail>Method "U(_,_)" is [ lower : number, upper : number | divisor ::= ((1<<32)) ÷ (upper - lower)→double; map a pRNG through [i : integer | (i ÷ divisor) + lower] ];

Method "a Marsaglia polar sampler" is [ generator for [ yield : [double]→⊤ | source ::= U(-1, 1); Repeat [ x ::= take 1 from source[1]; y ::= take 1 from source[1]; s ::= x^2 + y^2; If 0 < s < 1 then [ factor ::= ((-2 × ln s) ÷ s) ^ 0.5; yield(x × factor); yield(y × factor); ]; ] ] ];

// the default distribution has mean 0 and std dev 1.0, so we scale the values sampler ::= map a Marsaglia polar sampler through [d : double | d ÷ 2.0 + 1.0]; values ::= take 1000 from sampler;</lang>

AWK

One-liner: <lang awk>$ awk 'func r(){return sqrt(-2*log(rand()))*cos(6.2831853*rand())}BEGIN{for(i=0;i<1000;i++)s=s" "1+0.5*r();print s}'</lang>

Readable version: <lang awk> function r() {

 return sqrt( -2*log( rand() ) ) * cos(6.2831853*rand() )

}

BEGIN {

 n=1000
 for(i=0;i<n;i++) {
   x = 1 + 0.5*r()
   s = s" "x
 }
 print s

} </lang>

Output:

first few values only

0.783753 1.16682 1.17989 1.14975 1.34784 0.29296 0.979227 1.04402 0.567835 1.58812 0.465559 1.27186 0.324533 0.725827 -0.0626549 0.632273 1.0145 1.3387 0.861667 1.04147 1.2576 1.02497 0.58453 0.9619 1.26902 0.851048 -0.126259 0.863256

...

BASIC

Works with: QuickBasic version 4.5

RANDOMIZE TIMER 'seeds random number generator with the system time
pi = 3.141592653589793#
DIM a(1 TO 1000) AS DOUBLE
CLS
FOR i = 1 TO 1000
   a(i) = 1 + SQR(-2 * LOG(RND)) * COS(2 * pi * RND)
NEXT i

Applesoft BASIC

The Commodore BASIC code works in Applesoft BASIC.

BASIC256

Translation of: FreeBASIC

<lang BASIC256># Generates normally distributed random numbers with mean 0 and standard deviation 1 function randomNormal() return cos(2.0 * pi * rand) * sqr(-2.0 * log(rand)) end function

dim r(1000) sum = 0.0

Generate 1000 normally distributed random numbers
with mean 1 and standard deviation 0.5
and calculate their sum

for i = 0 to 999 r[i] = 1.0 + randomNormal() / 2.0 sum += r[i] next i

mean = sum / 1000.0 sum = 0.0

Now calculate their standard deviation

for i = 0 to 999 sum += (r[i] - mean) ^ 2.0 next i sd = sqr(sum/1000.0)

print "Mean is "; mean print "Standard Deviation is "; sd end</lang>

Output:

Mean is               1.002092
Standard Deviation is 0.4838570687

Commodore BASIC

<lang bbcbasic> 10 DIM AR(999): DIM DE(999) 20 FOR N = 0 TO 999 30 AR(N)= 0 + SQR(-1.3*LOG(RND(1))) * COS(1.2*PI*RND(1)) 40 NEXT N 50 : 60 REM SUM 70 LET SU = 0 80 FOR N = 0 TO 999 90 LET SU = SU + AR(N) 100 NEXT N 110 : 120 REM MEAN 130 LET ME= 0 140 LET ME = SU/1000 150 : 160 REM DEVIATION 170 FOR N = 0 TO 999 180 T = AR(N)-ME: REM SUBTRACT MEAN FROM NUMBER 190 T = T * T: REM SQUARE THE RESULT 200 DE(N) = T : REM STORE IN DEVIATION ARRAY 210 NEXT N 220 LET DS=0: REM SUM OF DEVIATION ARRAY 230 FOR N = 0 TO 999 240 LET DS = DS + DE(N) 250 NEXT N 260 LET DM=0: REM MEAN OF DEVIATION ARRAY 270 LET DM = DS / 1000 280 LET DE = 0: 290 LET DE = SQR(DM) 300 : 310 PRINT "MEAN = "ME 320 PRINT "STANDARD DEVIATION ="DE 330 END </lang>

BBC BASIC

<lang bbcbasic> DIM array(999)

     FOR number% = 0 TO 999
       array(number%) = 1.0 + 0.5 * SQR(-2*LN(RND(1))) * COS(2*PI*RND(1))
     NEXT
     
     mean = SUM(array()) / (DIM(array(),1) + 1)
     array() -= mean
     stdev = MOD(array()) / SQR(DIM(array(),1) + 1)
     
     PRINT "Mean = " ; mean
     PRINT "Standard deviation = " ; stdev</lang>

Output:

Mean = 1.01848064
Standard deviation = 0.503551814

C

<lang c>#include <stdlib.h>

include <math.h>
ifndef M_PI
define M_PI 3.14159265358979323846
endif

double drand() /* uniform distribution, (0..1] */ {

 return (rand()+1.0)/(RAND_MAX+1.0);

} double random_normal() /* normal distribution, centered on 0, std dev 1 */ {

 return sqrt(-2*log(drand())) * cos(2*M_PI*drand());

} int main() {

 int i;
 double rands[1000];
 for (i=0; i<1000; i++)
   rands[i] = 1.0 + 0.5*random_normal();
 return 0;

}</lang>

C#

Translation of: JavaScript

<lang csharp> private static double randomNormal() { return Math.Cos(2 * Math.PI * tRand.NextDouble()) * Math.Sqrt(-2 * Math.Log(tRand.NextDouble())); } </lang>

Then the methods in Random numbers#Metafont are used to calculate the average and the Standard Deviation: <lang csharp> static Random tRand = new Random();

static void Main(string[] args) { double[] a = new double[1000];

double tAvg = 0; for (int x = 0; x < a.Length; x++) { a[x] = randomNormal() / 2 + 1; tAvg += a[x]; }

tAvg /= a.Length; Console.WriteLine("Average: " + tAvg.ToString());

double s = 0; for (int x = 0; x < a.Length; x++) { s += Math.Pow((a[x] - tAvg), 2); } s = Math.Sqrt(s / 1000);

Console.WriteLine("Standard Deviation: " + s.ToString());

Console.ReadLine(); } </lang>

An example result:

Average: 1,00510073053613
Standard Deviation: 0,502540443430955

C++

Works with: C++11

The new C++ standard looks very similar to the Boost library example below.

<lang cpp>#include <random>

include <functional>
include <vector>
include <algorithm>

using namespace std;

int main() {

 random_device seed;
 mt19937 engine(seed());
 normal_distribution<double> dist(1.0, 0.5);
 auto rnd = bind(dist, engine);

 vector<double> v(1000);
 generate(v.begin(), v.end(), rnd);
 return 0;

}</lang>

Works with: C++03

<lang cpp>#include <cstdlib> // for rand

include <cmath> // for atan, sqrt, log, cos
include <algorithm> // for generate_n

double const pi = 4*std::atan(1.0);

// simple functor for normal distribution class normal_distribution { public:

 normal_distribution(double m, double s): mu(m), sigma(s) {}
 double operator() const // returns a single normally distributed number
 {
   double r1 = (std::rand() + 1.0)/(RAND_MAX + 1.0); // gives equal distribution in (0, 1]
   double r2 = (std::rand() + 1.0)/(RAND_MAX + 1.0);
   return mu + sigma * std::sqrt(-2*std::log(r1))*std::cos(2*pi*r2);
 }

private:

 const double mu, sigma;

};

int main() {

 double array[1000];
 std::generate_n(array, 1000, normal_distribution(1.0, 0.5));
 return 0;

}</lang>

Library: Boost

This example used Mersenne Twister generator. It can be changed by changing the typedef.

include <vector>
include "boost/random.hpp"
include "boost/generator_iterator.hpp"
include <boost/random/normal_distribution.hpp>
include <algorithm>

typedef boost::mt19937 RNGType; ///< mersenne twister generator

int main() {

   RNGType rng;
   boost::normal_distribution<> rdist(1.0,0.5); /**< normal distribution 
                          with mean of 1.0 and standard deviation of 0.5 */

   boost::variate_generator< RNGType, boost::normal_distribution<> >
                   get_rand(rng, rdist);

   std::vector<double> v(1000);
   generate(v.begin(),v.end(),get_rand);
   return 0;

} </lang>

Clojure

<lang lisp>(import '(java.util Random)) (def normals

 (let [r (Random.)]
   (take 1000 (repeatedly #(-> r .nextGaussian (* 0.5) (+ 1.0))))))</lang>

COBOL

      IDENTIFICATION DIVISION.
      PROGRAM-ID. RANDOM.
      AUTHOR.  Bill Gunshannon
      INSTALLATION.  Home.
      DATE-WRITTEN.  14 January 2022.
     ************************************************************
     ** Program Abstract:
     **   Able to get the Mean to be really close to 1.0 but
     **     couldn't get the Standard Deviation any closer than
     **     .3 to .4.
     ************************************************************
      
      DATA DIVISION.
      
      WORKING-STORAGE SECTION.
      
      01  Sample-Size          PIC 9(5)         VALUE 1000.
      01  Total                PIC 9(10)V9(5)  VALUE 0.0.
      01  Arith-Mean           PIC 999V999  VALUE 0.0.
      01  Std-Dev              PIC 999V999  VALUE 0.0.
      01  Seed                 PIC 999V999.
      01  TI                   PIC 9(8).

      01  Idx                  PIC 99999     VALUE 0.
      01  Intermediate         PIC 9(10)V9(5)  VALUE 0.0.
      01  Rnd-Work.
          05  Rnd-Tbl 
                  OCCURS 1 TO 99999 TIMES DEPENDING ON Sample-Size.
              10  Rnd              PIC 9V9999999  VALUE 0.0.
      
      PROCEDURE DIVISION.
      
      Main-Program.
          ACCEPT TI FROM TIME.
          MOVE FUNCTION RANDOM(TI) TO Seed.
          PERFORM WITH TEST AFTER VARYING Idx 
                  FROM 1 BY 1 
                  UNTIL Idx = Sample-Size
             COMPUTE Intermediate = 
                          (FUNCTION RANDOM() * 2.01) 
             MOVE Intermediate TO Rnd(Idx)
          END-PERFORM.
          PERFORM WITH TEST AFTER VARYING Idx 
                  FROM 1 BY 1 
                  UNTIL Idx = Sample-Size 
             COMPUTE Total = Total + Rnd(Idx)
          END-PERFORM.

          COMPUTE Arith-Mean = Total / Sample-Size.
          DISPLAY "Mean: " Arith-Mean.

          PERFORM WITH TEST AFTER VARYING Idx
                  FROM 1 BY 1
                  UNTIL Idx = Sample-Size
             COMPUTE Intermediate = 
                     Intermediate + (Rnd(Idx) - Arith-Mean) ** 2
          END-PERFORM.
             COMPUTE Std-Dev = Intermediate / Sample-Size.

          DISPLAY "Std-Dev: " Std-Dev.

          STOP RUN.
      
      END PROGRAM RANDOM.

</lang>

Common Lisp

<lang lisp>(loop for i from 1 to 1000

     collect (1+ (* (sqrt (* -2 (log (random 1.0)))) (cos (* 2 pi (random 1.0))) 0.5)))</lang>

Crystal

Translation of: Phix

<lang ruby>n, mean, sd, tau = 1000, 1, 0.5, (2 * Math::PI) array = Array.new(n) { mean + sd * Math.sqrt(-2 * Math.log(rand)) * Math.cos(tau * rand) }

mean = array.sum / array.size standev = Math.sqrt( array.sum{ |x| (x - mean) ** 2 } / array.size ) puts "mean = #{mean}, standard deviation = #{standev}"</lang>

Output:

mean = 1.0093442539237896, standard deviation = 0.504694489463623

D

<lang d>import std.stdio, std.random, std.math;

struct NormalRandom {

   double mean, stdDev;

   // Necessary because it also defines an opCall.
   this(in double mean_, in double stdDev_) pure nothrow {
       this.mean = mean_;
       this.stdDev = stdDev_;
   }

   double opCall() const /*nothrow*/ {
       immutable r1 = uniform01, r2 = uniform01; // Not nothrow.
       return mean + stdDev * sqrt(-2 * r1.log) * cos(2 * PI * r2);
   }

}

void main() {

   double[1000] array;
   auto nRnd = NormalRandom(1.0, 0.5);
   foreach (ref x; array)
       //x = nRnd;
       x = nRnd();

}</lang>

Alternative Version

(Untested)

Library: tango

<lang d>import tango.math.random.Random;

void main() {

   double[1000] list;
   auto r = new Random();
   foreach (ref l; list) {
       r.normalSource!(double)()(l);
       l = 1.0 + 0.5 * l;
   }

}</lang>

Delphi

Delphi has RandG function which generates random numbers with normal distribution using Marsaglia-Bray algorithm:

<lang Delphi>program Randoms;

{$APPTYPE CONSOLE}

uses

 Math;

var

 Values: array[0..999] of Double;
 I: Integer;

begin // Randomize; Commented to obtain reproducible results

 for I:= Low(Values) to High(Values) do
   Values[I]:= RandG(1.0, 0.5);  // Mean = 1.0, StdDev = 0.5
 Writeln('Mean          = ', Mean(Values):6:4);
 Writeln('Std Deviation = ', StdDev(Values):6:4);
 Readln;

end.</lang>

Output:

Mean          = 1.0098
Std deviation = 0.5016

DWScript

<lang delphi>var values : array [0..999] of Float; var i : Integer;

for i := values.Low to values.High do

  values[i] := RandG(1, 0.5);</lang>

E

<lang e>accum [] for _ in 1..1000 { _.with(entropy.nextGaussian()) }</lang>

EasyLang

<lang>for i range 1000

 a[] &= 1 + 0.5 * sqrt (-2 * logn randomf) * cos (360 * randomf)

. print a[]</lang>

Eiffel

<lang eiffel> class APPLICATION

inherit ARGUMENTS

create make

feature {NONE} -- Initialization

l_time: TIME l_seed: INTEGER math:DOUBLE_MATH rnd:RANDOM Size:INTEGER once Result:= 1000 end

make -- Run application. local ergebnis:ARRAY[DOUBLE] tavg: DOUBLE x: INTEGER tmp: DOUBLE text : STRING

do -- initialize random generator create l_time.make_now

    		        l_seed := l_time.hour
     		        l_seed := l_seed * 60 + l_time.minute
     		        l_seed := l_seed * 60 + l_time.second
     		        l_seed := l_seed * 1000 + l_time.milli_second
     		        create rnd.set_seed (l_seed)

-- initialize random number container and math create ergebnis.make_filled (0.0, 1, size) tavg := 0; create math

from x := 1 until x > ergebnis.count loop tmp := randomNormal / 2 + 1 tavg := tavg + tmp ergebnis.enter (tmp , x) x := x + 1 end

tavg := tavg / ergebnis.count text := "Average: " text.append_double (tavg) text.append ("%N") print(text)

tmp := 0 from x:= 1 until x > ergebnis.count loop tmp := tmp + (ergebnis.item (x) - tavg)^2 x := x + 1 end

tmp := math.sqrt (tmp / ergebnis.count) text := "Standard Deviation: " text.append_double (tmp) text.append ("%N") print(text)

end

randomNormal:DOUBLE

local

     		        first: DOUBLE
     		        second: DOUBLE

do

                       rnd.forth
                       first := rnd.double_item
                       rnd.forth
                       second := rnd.double_item

                       Result := math.cosine (2 * math.pi * first) * math.sqrt (-2 * math.log (second))

end end </lang>

Example Result

Average: 1.0079398405028137
Standard Deviation: 0.49042787564453988

Elena

Translation of: C#

ELENA 4.1 : <lang elena>import extensions; import extensions'math;

randomNormal() {

   ^ cos(2 * Pi_value * randomGenerator.nextReal()) 
                     * sqrt(-2 * ln(randomGenerator.nextReal()))

}

public program() {

   real[] a := new real[](1000);

   real tAvg := 0;
   for (int x := 0, x < a.Length, x += 1)
   {
       a[x] := (randomNormal()) / 2 + 1;
       tAvg += a[x]
   };

   tAvg /= a.Length;
   console.printLine("Average: ", tAvg);

   real s := 0;
   for (int x := 0, x < a.Length, x += 1)
   {
       s += power(a[x] - tAvg, 2)
   };

   s := sqrt(s / 1000);

   console.printLine("Standard Deviation: ", s);

   console.readChar()

}</lang>

Output:

Average: 0.9842420481571
Standard Deviation: 0.5109070975558

Elixir

<lang elixir>defmodule Random do

 def normal(mean, sd) do
   {a, b} = {:rand.uniform, :rand.uniform}
   mean + sd * (:math.sqrt(-2 * :math.log(a)) * :math.cos(2 * :math.pi * b)) 
 end

end

std_dev = fn (list) ->

           mean = Enum.sum(list) / length(list)
           sd = Enum.reduce(list, 0, fn x,acc -> acc + (x-mean)*(x-mean) end) / length(list)
                |> :math.sqrt
           IO.puts "Mean: #{mean},\tStdDev: #{sd}"
         end

xs = for _ <- 1..1000, do: Random.normal(1.0, 0.5) std_dev.(xs)</lang>

Output:

Mean: 1.009079383094275,        StdDev: 0.4991894476975088

used Erlang function :rand.normal <lang elixir>xs = for _ <- 1..1000, do: 1.0 + :rand.normal * 0.5 std_dev.(xs)</lang>

Output:

Mean: 0.9955701150615597,       StdDev: 0.5036412260426065

Erlang

Works with: Erlang

<lang erlang> mean(Values) ->

   mean(tl(Values), hd(Values), 1).

mean([], Acc, Length) ->

   Acc / Length;

mean(Values, Acc, Length) ->

   mean(tl(Values), hd(Values)+Acc, Length+1).

variance(Values) ->

   Mean = mean(Values),
   variance(Values, Mean, 0) / length(Values).

variance([], _, Acc) ->

   Acc;

variance(Values, Mean, Acc) ->

   Diff = hd(Values) - Mean,
   DiffSqr = Diff * Diff,
   variance(tl(Values), Mean, Acc + DiffSqr).

stddev(Values) ->

   math:sqrt(variance(Values)).

normal(Mean, StdDev) ->

   U = random:uniform(),
   V = random:uniform(),
   Mean + StdDev * ( math:sqrt(-2 * math:log(U)) * math:cos(2 * math:pi() * V) ).  % Erlang's math:log is the natural logarithm.

main(_) ->

   X = [ normal(1.0, 0.5) || _ <- lists:seq(1, 1000) ],
   io:format("mean = ~w\n", [mean(X)]),
   io:format("stddev = ~w\n", [stddev(X)]).

</lang>

Output:

mean = 1.0118289913718608
stddev = 0.5021636849524854

ERRE

<lang> PROGRAM DISTRIBUTION

! ! for rosettacode.org !

! formulas taken from TI-59 Master Library manual

CONST NUM_ITEM=1000

!VAR SUMX#,SUMX2#,R1#,R2#,Z#,I%

DIM A#[1000]

BEGIN ! seeds random number generator with system time

  RANDOMIZE(TIMER)

  PRINT(CHR$(12);)  !CLS
  SUMX#=0  SUMX2#=0

  FOR I%=1 TO NUM_ITEM DO
     R1#=RND(1)  R2#=RND(1)
     Z#=SQR(-2*LOG(R1#))*COS(2*π*R2#)
     A#[I%]=Z#/2+1   ! I want a normal distribution with
                     !      mean=1 and std.dev=0.5
     SUMX#+=A#[I%]  SUMX2#+=A#[I%]*A#[I%]
  END FOR

  Z#=SUMX#/NUM_ITEM

  PRINT("Average is";Z#)
  PRINT("Standard dev. is";SQR(SUMX2#/NUM_ITEM-Z#*Z#))

END PROGRAM </lang>

Euler Math Toolbox

<lang Euler Math Toolbox> >v=normal(1,1000)*0.5+1; >mean(v), dev(v)

1.00291801071
0.498226876528

</lang>

Euphoria

Translation of: PureBasic

<lang euphoria>include misc.e

function RandomNormal()

   atom x1, x2
   x1 = rand(999999) / 1000000
   x2 = rand(999999) / 1000000
   return sqrt(-2*log(x1)) * cos(2*PI*x2)

end function

constant n = 1000 sequence s s = repeat(0,n) for i = 1 to n do

   s[i] = 1 + 0.5 * RandomNormal()

end for</lang>

F#

<lang fsharp> let n = MathNet.Numerics.Distributions.Normal(1.0,0.5) List.init 1000 (fun _->n.Sample()) </lang>

Output:

  [0.734433576; 1.54225304; 0.4407528678; 1.177675412; 0.4318617021;
   0.6026656337; 0.769764924; 1.104693934; 0.6297500925; 0.9594598077;
   1.684736389; 1.160376323; 0.883354356; 0.9513968363; 0.9727698268;
   0.5315570949; 0.9599239266; 1.564976755; 0.7232002879; 1.084139442;
   1.220914517; 0.3553085946; 1.112549824; 1.989443553; 0.5752307543;
   1.156682549; 0.7886670467; 0.02050745923; 1.532060208; 1.18789591;
   1.408946777; 1.038812004; 1.724679503; 1.671565045; 1.266831442;
   1.363611654; 1.705819067; 0.5772366328; 0.4503488498; 1.496891481;
   0.9831877282; 0.3845460366; 0.8253240671; 1.298969969; 0.4265904553;
   0.9303696876; 0.445003361; 0.753175816; 0.6143534043; 1.059982235;
   0.7143206784; 0.2233328038; 1.005178481; 0.7697392436; 0.5904948577;
   0.5127953044; 0.6467346747; 0.7929387604; -0.1501790761; 0.8750780903;
   0.941704369; 1.37941579; 0.4739006145; 1.998886344; 1.219428519;
   0.06270791476; 1.097739804; 0.7584232803; 1.042177217; 1.166561247;
   1.502357164; 1.171525776; 0.1528807432; 0.2289389756; 1.36208422;
   0.3714421124; 1.299571092; 1.171553369; 1.317807265; 1.616662281;
   1.724223246; 1.059580642; 1.270520918; -0.1827677907; 1.938593232;
   1.420362143; 1.888357595; 0.7851629936; 0.7080554899; 0.7747215818;
   1.403719877; 0.5765950249; 1.275206565; 0.6292054813; 1.525562798;
   0.6224640457; 0.8524078517; 0.7646595627; 0.6799834691; 0.773111053; ...]

==F#==

<lang fsharp>let gaussianRand count =

   let o = new System.Random()
   let pi = System.Math.PI
   let gaussrnd = 
       (fun _ -> 1. + 0.5 * sqrt(-2. * log(o.NextDouble())) * cos(2. * pi * o.NextDouble()))
   [ for i in {0 .. (int count)} -> gaussrnd() ]</lang>

Factor

<lang factor>USING: random ; 1000 [ 1.0 0.5 normal-random-float ] replicate</lang>

Falcon

<lang falcon>a = [] for i in [0:1000] : a+= norm_rand_num()

function norm_rand_num()

  pi = 2*acos(0)
  return 1 + (cos(2 * pi * random()) * pow(-2 * log(random()) ,1/2)) /2

end</lang>

Fantom

Two solutions. The first uses Fantom's random-number generator, which produces a uniform distribution. So, convert to a normal distribution using a formula:

<lang fantom> class Main {

 static const Float PI := 0.0f.acos * 2  // we need to precompute PI

 static Float randomNormal ()
 {
   return (Float.random * PI * 2).cos * (Float.random.log * -2).sqrt
 }

 public static Void main ()
 {
   mean := 1.0f
   sd := 0.5f
   Float[] values := [,] // this is the collection to fill with random numbers
   1000.times { values.add (randomNormal * sd + mean) }
 }

} </lang>

The second calls out to Java's Gaussian random-number generator:

<lang fantom> using [java] java.util::Random

class Main {

 Random generator := Random()

 Float randomNormal ()
 {
   return generator.nextGaussian
 }

 public static Void main ()
 {
   rnd := Main()  // create an instance of Main class, which holds the generator
   mean := 1.0f
   sd := 0.5f
   Float[] values := [,] // this is the collection to fill with random numbers
   1000.times { values.add (rnd.randomNormal * sd + mean) }
 }

} </lang>

Forth

Works with: gforth version 0.6.2

<lang forth>require random.fs here to seed

-1. 1 rshift 2constant MAX-D \ or s" MAX-D" ENVIRONMENT? drop

frnd ( -- f ) \ uniform distribution 0..1

 rnd rnd dabs d>f MAX-D d>f f/ ;

frnd-normal ( -- f ) \ centered on 0, std dev 1

 frnd pi f* 2e f* fcos
 frnd fln -2e f* fsqrt f* ;

,normals ( n -- ) \ store many, centered on 1, std dev 0.5

 0 do frnd-normal 0.5e f* 1e f+ f, loop ;

create rnd-array 1000 ,normals</lang>

For newer versions of gforth (tested on 0.7.3), it seems you need to use HERE SEED ! instead of HERE TO SEED, because SEED has been made a variable instead of a value.

<lang>rnd rnd dabs d>f</lang> is necessary, but surprising and definitely not well documented / perhaps not compliant.

Fortran

Works with: Fortran version 90 and later

<lang fortran>PROGRAM Random

 INTEGER, PARAMETER :: n = 1000
 INTEGER :: i
 REAL :: array(n), pi, temp, mean = 1.0, sd = 0.5

 pi = 4.0*ATAN(1.0)
 CALL RANDOM_NUMBER(array) ! Uniform distribution

! Now convert to normal distribution

 DO i = 1, n-1, 2
   temp = sd * SQRT(-2.0*LOG(array(i))) * COS(2*pi*array(i+1)) + mean
   array(i+1) = sd * SQRT(-2.0*LOG(array(i))) * SIN(2*pi*array(i+1)) + mean
   array(i) = temp
 END DO

! Check mean and standard deviation

 mean = SUM(array)/n
 sd = SQRT(SUM((array - mean)**2)/n)

 WRITE(*, "(A,F8.6)") "Mean = ", mean
 WRITE(*, "(A,F8.6)") "Standard Deviation = ", sd

END PROGRAM Random</lang>

Output:

 Mean = 0.995112
 Standard Deviation = 0.503373

Free Pascal

Free Pascal provides the randg function in the RTL math unit that produces Gaussian-distributed random numbers with the Box-Müller algorithm.

<lang pascal> function randg(mean,stddev: float): float; </lang>

FreeBASIC

<lang freebasic>' FB 1.05.0 Win64

Const pi As Double = 3.141592653589793 Randomize

' Generates normally distributed random numbers with mean 0 and standard deviation 1 Function randomNormal() As Double

 Return Cos(2.0 * pi * Rnd) * Sqr(-2.0 * Log(Rnd))

End Function

Dim r(0 To 999) As Double Dim sum As Double = 0.0

' Generate 1000 normally distributed random numbers ' with mean 1 and standard deviation 0.5 ' and calculate their sum For i As Integer = 0 To 999

  r(i) = 1.0 + randomNormal/2.0
  sum += r(i)

Dim sd As Double sum = 0.0 ' Now calculate their standard deviation For i As Integer = 0 To 999

 sum += (r(i) - mean) ^ 2.0

Next sd = Sqr(sum/1000.0)

Print "Mean is "; mean Print "Standard Deviation is"; sd Print Print "Press any key to quit" Sleep</lang> Sample result:

Output:

Mean is               1.000763573902885
Standard Deviation is 0.500653063426955

Frink

<lang frink>a = new array[[1000], {|x| randomGaussian[1, 0.5]}]</lang>

FutureBasic

Note: To generate the random number, rather than using FB's native "rnd" function, this code wraps C code into the RandomZeroToOne function. <lang futurebasic> window 1

local fn RandomZeroToOne as double dim as double result BeginCCode

 result = (double)( (rand() % 100000 ) * 0.00001 );

EndC end fn = result

local fn RandomGaussian as double dim as double r

r = fn RandomZeroToOne end fn = 1 + .5 * ( sqr( -2 * log(r) ) * cos( 2 * pi * r ) )

long i double mean, std, a(1000)

for i = 1 to 1000 a(i) = fn RandomGaussian mean += a(i) next mean = mean / 1000

for i = 1 to 1000 std += ( a(i) - mean )^2 next std = std / 1000

print " Average: "; mean print "Standard Deviation: "; std

HandleEvents

Go

This solution uses math/rand package in the standard library. See also though the subrepository rand package at https://godoc.org/golang.org/x/exp/rand, which also has a NormFloat64 and has a rand source with a number of advantages over the one in standard library. <lang go>package main

import (

   "fmt"
   "math"
   "math/rand"
   "strings"
   "time"

)

const mean = 1.0 const stdv = .5 const n = 1000

func main() {

   var list [n]float64
   rand.Seed(time.Now().UnixNano())
   for i := range list {
       list[i] = mean + stdv*rand.NormFloat64()
   }
   // show computed mean and stdv of list
   var s, sq float64
   for _, v := range list {
       s += v
   }
   cm := s / n
   for _, v := range list {
       d := v - cm
       sq += d * d
   }
   fmt.Printf("mean %.3f, stdv %.3f\n", cm, math.Sqrt(sq/(n-1)))
   // show histogram by hdiv divisions per stdv over +/-hrange stdv
   const hdiv = 3
   const hrange = 2
   var h [1 + 2*hrange*hdiv]int
   for _, v := range list {
       bin := hrange*hdiv + int(math.Floor((v-mean)/stdv*hdiv+.5))
       if bin >= 0 && bin < len(h) {
           h[bin]++
       }
   }
   const hscale = 10
   for _, c := range h {
       fmt.Println(strings.Repeat("*", (c+hscale/2)/hscale))
   }

}</lang>

Output:

mean 0.995, stdv 0.503
**
****
******
********
************
************
*************
************
**********
********
*****
***
**

Groovy

<lang groovy>rnd = new Random() result = (1..1000).inject([]) { r, i -> r << rnd.nextGaussian() }</lang>

Haskell

<lang haskell>import System.Random

pairs :: [a] -> [(a,a)] pairs (x:y:zs) = (x,y):pairs zs pairs _ = []

gauss mu sigma (r1,r2) =

 mu + sigma * sqrt (-2 * log r1) * cos (2 * pi * r2)

gaussians :: (RandomGen g, Random a, Floating a) => Int -> g -> [a] gaussians n g = take n $ map (gauss 1.0 0.5) $ pairs $ randoms g

result :: IO [Double] result = getStdGen >>= \g -> return $ gaussians 1000 g</lang>

Or using Data.Random from random-fu package: <lang haskell>replicateM 1000 $ normal 1 0.5</lang> To print them: <lang haskell>import Data.Random import Control.Monad

thousandRandomNumbers :: RVar [Double] thousandRandomNumbers = replicateM 1000 $ normal 1 0.5

main = do

  x <- sample thousandRandomNumbers
  print x</lang>

HicEst

<lang hicest>REAL :: n=1000, m=1, s=0.5, array(n)

pi = 4 * ATAN(1) array = s * (-2*LOG(RAN(1)))^0.5 * COS(2*pi*RAN(1)) + m </lang>

Icon and Unicon

The seed &random may be assigned in either language; either to randomly seed or to pick a fixed starting point. ?i is the random number generator, returning an integer from 0 to i - 1 for non-zero integer i. As a special case, ?0 yields a random floating point number from 0.0 <= r < 1.0

Note that Unicon randomly seeds it's generator. <lang icon> procedure main()

   local L
   L := list(1000)
   every L[1 to 1000] := 1.0 + 0.5 * sqrt(-2.0 * log(?0)) * cos(2.0 * &pi * ?0)

   every write(!L)

end </lang>

IDL

<lang idl>result = 1.0 + 0.5*randomn(seed,1000)</lang>

J

Solution: <lang j>urand=: ?@$ 0: zrand=: (2 o. 2p1 * urand) * [: %: _2 * [: ^. urand

1 + 0.5 * zrand 100</lang>

Alternative Solution:
Using the normal script from the stats/distribs addon. <lang j> require 'stats/distribs/normal'

  1 0.5 rnorm 1000

1.44868803 1.21548637 0.812460657 1.54295452 1.2470606 ...</lang>

Java

<lang java>double[] list = new double[1000]; double mean = 1.0, std = 0.5; Random rng = new Random(); for(int i = 0;i<list.length;i++) {

 list[i] = mean + std * rng.nextGaussian();

}</lang>

JavaScript

<lang javascript>function randomNormal() {

 return Math.cos(2 * Math.PI * Math.random()) * Math.sqrt(-2 * Math.log(Math.random()))

}

var a = [] for (var i=0; i < 1000; i++){

 a[i] = randomNormal() / 2 + 1

}</lang>

jq

Works with: jq version 1.4

Since jq is a purely functional language, it is convenient to define the pseudo-random number generator functions as filters whose inputs and outputs are arrays containing a "seed".

The following uses the same pseudo-random number generator as the Microsoft C Runtime (see Linear congruential generator).

'A Pseudo-Random Number Generator' <lang jq># 15-bit integers generated using the same formula as rand() from the Microsoft C Runtime.

The random numbers are in [0 -- 32767] inclusive.
Input: an array of length at least 2 interpreted as [count, state, ...]
Output: [count+1, newstate, r] where r is the next pseudo-random number.

def next_rand_Microsoft:

 .[0] as $count | .[1] as $state
 | ( (214013 * $state) + 2531011) % 2147483648 # mod 2^31
 | [$count+1 , ., (. / 65536 | floor) ] ;</lang>

'Box-Muller Method' <lang jq># Generate a single number following the normal distribution with mean 0, variance 1,

using the Box-Muller method: X = sqrt(-2 ln U) * cos(2 pi V) where U and V are uniform on [0,1].
Input: [n, state]
Output [n+1, nextstate, r]

def next_rand_normal:

 def u: next_rand_Microsoft | .[2] /= 32767; 
 u as $u1
 | ($u1 | u) as $u2
 | ((( (8*(1|atan)) * $u1[2]) | cos)
    * ((-2 * (($u2[2]) | log)) | sqrt)) as $r
 | [ (.[0]+1), $u2[1], $r] ;

Generate "count" arrays, each containing a random normal variate with the given mean and standard deviation.
Input: [count, state]
Output: [updatedcount, updatedstate, rnv]
where "state" is a seed and "updatedstate" can be used as a seed.

def random_normal_variate(mean; sd; count):

 next_rand_normal
 | recurse( if .[0] < count then next_rand_normal else empty end)
 | .[2] = (.[2] * sd) + mean;</lang>

Example The task can be completed using: [0,1] | random_normal_variate(1; 0.5; 1000) | .[2]

We show just the sample average and standard deviation: <lang jq>def summary:

 length as $l | add as $sum | ($sum/$l) as $a
 | reduce .[] as $x (0; . + ( ($x - $a) | .*. ))
 | [ $a, (./$l | sqrt)] ;

[ [0,1] | random_normal_variate(1; 0.5; 1000) | .[2] ] | summary</lang>

Output:

$ jq -n -c -f Random_numbers.jq
[0.9932830741018853,0.4977760644490579]

Julia

Julia's standard library provides a randn function to generate normally distributed random numbers (with mean 0 and standard deviation 0.5, which can be easily rescaled to any desired values): <lang julia>randn(1000) * 0.5 + 1</lang>

Kotlin

<lang scala>// version 1.0.6

import java.util.Random

fun main(args: Array<String>) {

   val r = Random()
   val da = DoubleArray(1000)
   for (i in 0 until 1000)  da[i] = 1.0 + 0.5 * r.nextGaussian()
   // now check actual mean and SD
   val mean = da.average()
   val sd = Math.sqrt(da.map { (it - mean) * (it - mean) }.average())
   println("Mean is $mean")
   println("S.D. is $sd")

}</lang> Sample output:

Output:

Mean is 1.0071784073168768
S.D. is 0.48567118114896807

LabVIEW

Works with: LabVIEW version 8.6

Liberty BASIC

<lang lb>dim a(1000) mean =1 sd =0.5 for i = 1 to 1000 ' throw 1000 normal variates

  a( i)  =mean +sd *( sqr( -2 * log( rnd( 0))) * cos( 2 * pi * rnd( 0)))

next i</lang>

Lingo

<lang Lingo>-- Returns a random float value in range 0..1 on randf ()

   n = random(the maxinteger)-1
   return n / float(the maxinteger-1)

end</lang>

<lang Lingo>normal = [] repeat with i = 1 to 1000

   normal.add(1 + sqrt(-2 * log(randf())) * cos(2 * PI * randf()) / 2)

end repeat</lang>

Lobster

Uses built-in rnd_gaussian <lang Lobster> let mean = 1.0 let stdv = 0.5 let count = 1000

// stats computes a running mean and variance // See Knuth TAOCP vol 2, 3rd edition, page 232

def stats(xs: [float]) -> float, float: // variance, mean

   var M = xs[0]
   var S = 0.0
   var n = 1.0
   for(xs.length - 1) i:
       let x = xs[i + 1]
       n = n + 1.0
       let mm = (x - M)
       M += mm / n
       S += mm * (x - M)
   return (if n > 0.0: S / n else: 0.0), M

def test_random_normal() -> [float]:

   rnd_seed(floor(seconds_elapsed() * 1000000))
   let r = vector_reserve(typeof return, count)
   for (count):
       r.push(rnd_gaussian() * stdv + mean)
   let cvar, cmean = stats(r)
   let cstdv = sqrt(cvar)
   print concat_string(["Mean: ", string(cmean), ", Std.Deviation: ", string(cstdv)], "")

test_random_normal() </lang>

Logo

Works with: UCB Logo

The earliest Logos only have a RANDOM function for picking a random non-negative integer. Many modern Logos have floating point random generators built-in. <lang logo>to random.float ; 0..1

 localmake "max.int lshift -1 -1
 output quotient random :max.int :max.int

end

to random.gaussian

 output product cos random 360  sqrt -2 / ln random.float

end

make "randoms cascade 1000 [fput random.gaussian / 2 + 1 ?] []</lang>

Lua

<lang lua>local list = {} for i = 1, 1000 do

 list[i] = 1 + math.sqrt(-2 * math.log(math.random())) * math.cos(2 * math.pi * math.random()) / 2

end</lang>

M2000 Interpreter

M2000 use a Wichmann - Hill Pseudo Random Number Generator. <lang M2000 Interpreter> Module CheckIt {

     Function StdDev (A()) {
         \\ A()  has a copy of values
           N=Len(A())
           if N<1 then Error "Empty Array"
           M=Each(A()) 
           k=0
           \\ make sum, dev same type as A(k)
           sum=A(k)-A(k)
           dev=sum
           \\ find mean
           While M {
                 sum+=Array(M)
           }
           Mean=sum/N
           \\ make a pointet to A()
           P=A()
           \\ subtruct from each item
           P-=Mean
           
           M=Each(P) 
           While M {
                 dev+=Array(M)*Array(M)
           }
           \\ as pointer to arrray
            =(if(dev>0->Sqrt(dev/N), 0), Mean)
     }
     Function randomNormal {
           \\ by default all numbers are double
           \\ cos() get degrees
         =1+Cos(360 * rnd) * Sqrt(-2 * Ln(rnd)) /2
     }  
     \\ fill array calling  randomNormal() for each item
     Dim A(1000)<<randomNormal()
     \\ we can pass a pointer to array and place it to stack of values
     DisplayMeanAndStdDeviation(A())  ' mean ~ 1 std deviation ~0.5
     \\ check M2000 rnd only
     Dim B(1000)<<rnd
     DisplayMeanAndStdDeviation(B())  ' mean ~ 0.5 std deviation ~0.28

     DisplayMeanAndStdDeviation((0,0,14,14))  ' mean = 7 std deviation = 7
     DisplayMeanAndStdDeviation((0,6,8,14))  ' mean = 7 std deviation = 5
     DisplayMeanAndStdDeviation((6,6,8,8))  ' mean = 7 std deviation = 1
     
     Sub DisplayMeanAndStdDeviation(A)
           \\ push to stack all items of an array (need an array pointer)
           Push  ! StdDev(A) 
           \\ read from strack two numbers
           Print "Mean is               "; Number
           Print "Standard Deviation is "; Number
     End Sub

} Checkit </lang>

Maple

<lang maple>with(Statistics): Sample(Normal(1, 0.5), 1000);</lang>

or

<lang maple>1+0.5*ArrayTools[RandomArray](1000,1,distribution=normal);</lang>

Mathematica/Wolfram Language

Built-in function RandomReal with built-in distribution NormalDistribution as an argument: <lang Mathematica>RandomReal[NormalDistribution[1, 1/2], 1000]</lang>

MATLAB

Native support : <lang MATLAB> mu = 1; sd = 0.5;

   x = randn(1000,1) * sd + mu;

</lang>

The statistics toolbox provides this function <lang MATLAB> x = normrnd(mu, sd, [1000,1]); </lang>

This script uses the Box-Mueller Transform to transform a number from the uniform distribution to a normal distribution of mean = mu0 and standard deviation = chi2.

<lang MATLAB>function randNum = randNorm(mu0,chi2, sz)

   radiusSquared = +Inf;

   while (radiusSquared >= 1)
       u = ( 2 * rand(sz) ) - 1;
       v = ( 2 * rand(sz) ) - 1;

       radiusSquared = u.^2 + v.^2;
   end

   scaleFactor = sqrt( ( -2*log(radiusSquared) )./ radiusSquared );
   randNum = (v .* scaleFactor .* chi2) + mu0;

end</lang>

Output: <lang MATLAB>>> randNorm(1,.5, [1000,1])

ans =

  0.693984121077029</lang>

Maxima

<lang maxima>load(distrib)$

random_normal(1.0, 0.5, 1000);</lang>

MAXScript

<lang maxscript>arr = #() for i in 1 to 1000 do (

   a = random 0.0 1.0
   b = random 0.0 1.0
   c = 1.0 + 0.5 * sqrt (-2*log a) * cos (360*b) -- Maxscript cos takes degrees
   append arr c

)</lang>

Metafont

Metafont has normaldeviate which produces pseudorandom normal distributed numbers with mean 0 and variance one. So the following complete the task:

<lang metafont>numeric col[];

m := 0; % m holds the mean, for testing purposes for i = 1 upto 1000:

 col[i] := 1 + .5normaldeviate;
 m := m + col[i];

endfor

% testing m := m / 1000; % finalize the computation of the mean

s := 0; % in s we compute the standard deviation for i = 1 upto 1000:

 s := s + (col[i] - m)**2;

endfor s := sqrt(s / 1000);

show m, s; % and let's show that really they get what we wanted end</lang>

A run gave

>> 0.99947
>> 0.50533

Assigning a value to the special variable randomseed will allow to have always the same sequence of pseudorandom numbers

MiniScript

<lang MiniScript>randNormal = function(mean=0, stddev=1)

   return mean + sqrt(-2 * log(rnd,2.7182818284)) * cos(2*pi*rnd) * stddev

end function

x = [] for i in range(1,1000)

   x.push randNormal(1, 0.5)

end for</lang>

Mirah

<lang mirah>import java.util.Random

list = double[999] mean = 1.0 std = 0.5 rng = Random.new 0.upto(998) do | i |

   list[i] = mean + std * rng.nextGaussian

end </lang>

МК-61/52

<lang>П7 <-> П8 1/x П6 ИП6 П9 СЧ П6 1/x ln ИП8 * 2 * КвКор ИП9 2 * пи

sin * ИП7 + С/П БП 05</lang>

Input: РY - variance, РX - expectation.

Or:

<lang>3 10^x П0 ПП 13 2 / 1 + С/П L0 03 С/П СЧ lg 2 /-/ * КвКор 2 пи ^ СЧ * * cos * В/О</lang>

to generate 1000 numbers with a mean of 1.0 and a standard deviation of 0.5.

Modula-3

Translation of: C

<lang modula3>MODULE Rand EXPORTS Main;

IMPORT Random; FROM Math IMPORT log, cos, sqrt, Pi;

VAR rands: ARRAY [1..1000] OF LONGREAL;

(* Normal distribution. *) PROCEDURE RandNorm(): LONGREAL =

 BEGIN
   WITH rand = NEW(Random.Default).init() DO
     RETURN 
       sqrt(-2.0D0 * log(rand.longreal())) * cos(2.0D0 * Pi * rand.longreal());
   END;
 END RandNorm;

BEGIN

 FOR i := FIRST(rands) TO LAST(rands) DO
   rands[i] := 1.0D0 + 0.5D0 * RandNorm();
 END;

END Rand.</lang>

Nanoquery

Translation of: Java

<lang nanoquery>list = {0} * 1000 mean = 1.0; std = 0.5 rng = new(Nanoquery.Util.Random)

for i in range(0, len(list) - 1)

       list[i] = mean + std * rng.getGaussian()

end</lang>

NetRexx

<lang NetRexx>/* NetRexx */ options replace format comments java crossref symbols nobinary

import java.math.BigDecimal import java.math.MathContext

-- prologue numeric digits 20

-- get input, set defaults parse arg dp mu sigma ec . if mu = | mu = '.' then mean = 1.0; else mean = mu if sigma = | sigma = '.' then stdDeviation = 0.5; else stdDeviation = sigma if dp = | dp = '.' then displayPrecision = 1; else displayPrecision = dp if ec = | ec = '.' then elements = 1000; else elements = ec

-- set up RNG = Random() numberList = java.util.List numberList = ArrayList()

-- generate list of random numbers loop for elements

 rn = mean + stdDeviation * RNG.nextGaussian()
 numberList.add(BigDecimal(rn, MathContext.DECIMAL128))
 end

-- report say "Mean: " mean say "Standard Deviation:" stdDeviation say "Precision: " displayPrecision say drawBellCurve(numberList, displayPrecision)

return

-- ----------------------------------------------------------------------------- method drawBellCurve(numberList = java.util.List, precision) static

 Collections.sort(numberList)
 val = BigDecimal
 lastN = 
 nextN = 
 loop val over numberList
   nextN = Rexx(val.toPlainString()).format(5, precision)
   select
     when lastN =  then nop
     when lastN \= nextN then say lastN
     otherwise nop
     end
   say '*\-'
   lastN = nextN
   end val
 say lastN

 return

</lang>

Output:

Mean:               1.0
Standard Deviation: 0.5
Precision:          1

*    2.7
**    2.5
*    2.4
***    2.3
*****    2.2
*******    2.1
*************    2.0
*************    1.9
*****************************    1.8
*************************    1.7
*************************************    1.6
******************************************************    1.5
********************************************    1.4
********************************************************************    1.3
*****************************************************************    1.2
**************************************************************************    1.1
*********************************************************************************************    1.0
*************************************************************    0.9
**********************************************************************    0.8
**************************************************************    0.7
***********************************************************************    0.6
**************************************************************    0.5
******************************************    0.4
*******************************    0.3
***************************    0.2
***************    0.1
*********    0.0
******   -0.1
***   -0.2
***   -0.3
*   -0.4
*   -0.6
**   -0.7

NewLISP

<lang NewLISP>(normal 1 .5 1000)</lang>

Nim

<lang nim>import random, stats, strformat

var rs: RunningStat

randomize()

for _ in 1..5:

 for _ in 1..1000: rs.push gauss(1.0, 0.5)
 echo &"mean: {rs.mean:.5f}    stdDev: {rs.standardDeviation:.5f}"

</lang>

Output:

mean: 1.01294    stdDev: 0.49692
mean: 1.00262    stdDev: 0.50028
mean: 0.99878    stdDev: 0.49662
mean: 0.99830    stdDev: 0.49820
mean: 1.00658    stdDev: 0.49703

Objeck

<lang objeck>bundle Default {

 class RandomNumbers {
   function : Main(args : String[]) ~ Nil {
     rands := Float->New[1000];
     for(i := 0; i < rands->Size(); i += 1;) {
       rands[i] := 1.0 + 0.5 * RandomNormal();
     };

     each(i : rands) {
       rands[i]->PrintLine();
     };
   }

   function : native : RandomNormal() ~ Float {
     return (2 * Float->Pi() * Float->Random())->Cos() * (-2 * (Float->Random()->Log()))->SquareRoot();
   }
 }

}</lang>

OCaml

<lang ocaml>let pi = 4. *. atan 1.;; let random_gaussian () =

 1. +. sqrt (-2. *. log (Random.float 1.)) *. cos (2. *. pi *. Random.float 1.);;

let a = Array.init 1000 (fun _ -> random_gaussian ());;</lang>

Octave

<lang octave>p = normrnd(1.0, 0.5, 1000, 1); disp(mean(p)); disp(sqrt(sum((p - mean(p)).^2)/numel(p)));</lang>

Output:

1.0209
0.51048

ooRexx

Translation of: REXX

version 1

<lang oorexx>/*REXX pgm gens 1,000 normally distributed #s: mean=1, standard dev.=0.5*/

 pi=RxCalcPi()                     /* get value of pi                */
 Parse Arg n seed .                /* allow specification of N & seed*/
 If n==|n==',' Then
   n=1000                          /* N  is the size of the array.   */
 If seed\== Then
   Call random,,seed               /* use seed for repeatable RANDOM#*/
 mean=1                            /* desired new mean (arith. avg.) */
 sd=1/2                            /* desired new standard deviation.*/
 Do g=1 For n                      /* generate N uniform random nums.*/
   n.g=random(0,1e5)/1e5           /* REXX gens uniform rand integers*/
   End

 Say '              old mean=' mean()
 Say 'old standard deviation=' stddev()
 Say
 Do j=1 To n-1 By 2
   m=j+1
                                   /*use Box-Muller method           */
   _=sd*RxCalcPower(-2*RxCalcLog(n.j),.5)*RxCalcCos(2*pi*n.m,,'R')+mean
   n.m=sd*RxCalcpower(-2*RxCalcLog(n.j),.5)*RxCalcSin(2*pi*n.m,,'R')+,
 mean                              /* rand # must be 0???1.          */
   n.j=_
   End                             /* j                              */
 Say '              new mean=' mean()
 Say 'new standard deviation=' stddev()
 Exit

mean:

 _=0
 Do k=1 For n
   _=_+n.k
   End
 Return _/n

stddev:

 _avg=mean()
 _=0
 Do k=1 For n
   _=_+(n.k-_avg)**2
   End
 Return RxCalcPower(_/n,.5)

requires rxmath library</lang>

Output:

              old mean= 0.49830002
old standard deviation= 0.283199568

              new mean= 1.00377404
new standard deviation= 0.501444536

version 2

Using the nice function names in the algorithm. <lang oorexx>/*REXX pgm gens 1,000 normally distributed #s: mean=1, standard dev.=0.5*/

 pi=RxCalcPi()                     /* get value of pi                */
 Parse Arg n seed .                /* allow specification of N & seed*/
 If n==|n==',' Then
   n=1000                          /* N  is the size of the array.   */
 If seed\== Then
   Call random,,seed               /* use seed for repeatable RANDOM#*/
 mean=1                            /* desired new mean (arith. avg.) */
 sd=1/2                            /* desired new standard deviation.*/
 Do g=1 For n                      /* generate N uniform random nums.*/
   n.g=random(0,1e5)/1e5           /* REXX gens uniform rand integers*/
   End

 Say '              old mean=' mean()
 Say 'old standard deviation=' stddev()
 Say
 Do j=1 To n-1 By 2
   m=j+1
                                   /*use Box-Muller method           */
   _=sd*sqrt(-2*ln(n.j))*cos(2*pi*n.m)+mean
   n.m=sd*sqrt(-2*ln(n.j))*sin(2*pi*n.m)+mean
   n.j=_
   End
 Say '              new mean=' mean()
 Say 'new standard deviation=' stddev()
 Exit

mean:

 _=0
 Do k=1 For n
   _=_+n.k
   End
 Return _/n

stddev:

 _avg=mean()
 _=0
 Do k=1 For n
   _=_+(n.k-_avg)**2
   End
 Return sqrt(_/n)

sqrt: Return RxCalcSqrt(arg(1)) ln: Return RxCalcLog(arg(1)) cos: Return RxCalcCos(arg(1),,'R') sin: Return RxCalcSin(arg(1),,'R')

requires rxmath library</lang>

PARI/GP

<lang parigp>rnormal()={ my(pr=32*ceil(default(realprecision)*log(10)/log(4294967296)),u1=random(2^pr)*1.>>pr,u2=random(2^pr)*1.>>pr); sqrt(-2*log(u1))*cos(2*Pi*u2) \\ in previous version "u1" instead of "u2" was used --> has given crap distribution \\ Could easily be extended with a second normal at very little cost. }; vector(1000,unused,rnormal()/2+1)</lang>

Pascal

The following function calculates Gaussian-distributed random numbers with the Box-Müller algorithm: <lang pascal> function rnorm (mean, sd: real): real;

{Calculates Gaussian random numbers according to the Box-Müller approach}

var

 u1, u2: real;

begin

 u1 := random;
 u2 := random;
 rnorm := mean * abs(1 + sqrt(-2 * (ln(u1))) * cos(2 * pi * u2) * sd); 
      /* error !?! Shouldn't it be "mean +" instead of "mean *" ? */

end; </lang>

Delphi and Free Pascal support implement a randg function that delivers Gaussian-distributed random numbers.

Perl

<lang perl>my $PI = 2 * atan2 1, 0;

my @nums = map {

   1 + 0.5 * sqrt(-2 * log rand) * cos(2 * $PI * rand)

} 1..1000;</lang>

Phix

Translation of: Euphoria

function RandomNormal()
    return sqrt(-2*log(rnd())) * cos(2*PI*rnd())
end function
 
sequence s = repeat(0,1000)
for i=1 to length(s) do
    s[i] = 1 + 0.5 * RandomNormal()
end for

Phixmonti

<lang Phixmonti>include ..\Utilitys.pmt

def RandomNormal

   drop rand log -2 * sqrt 2 pi * rand * cos * 0.5 * 1 +

enddef

1000 var n 0 n repeat

getid RandomNormal map

dup sum n / var mean "Mean: " print mean print nl

0 swap n for

   get mean - 2 power rot + swap

endfor swap n / sqrt "Standard deviation: " print print</lang>

PHP

<lang php>function random() {

   return mt_rand() / mt_getrandmax();

}

$pi = pi(); // Set PI

$a = array(); for ($i = 0; $i < 1000; $i++) {

   $a[$i] = 1.0 + ((sqrt(-2 * log(random())) * cos(2 * $pi * random())) * 0.5);

}</lang>

Picat

  _ = random2(), % random seed
  G = [gaussian_dist(1,0.5) : _ in 1..1000],
  println(first_10=G[1..10]),
  println([mean=avg(G),stdev=stdev(G)]),
  nl.

% Gaussian (Normal) distribution, Box-Muller algorithm gaussian01() = Y =>

   U = frand(0,1),
   V = frand(0,1),
   Y = sqrt(-2*log(U))*sin(2*math.pi*V).

gaussian_dist(Mean,Stdev) = Mean + (gaussian01() * Stdev).

% Variance of Xs variance(Xs) = Variance =>

   Mu = avg(Xs),
   N  = Xs.len,
   Variance = sum([ (X-Mu)**2 : X in Xs ]) / N.

% Standard deviation stdev(Xs) = sqrt(variance(Xs)).</lang>

Output:

first_10 = [1.639965415776091,0.705425965005482,0.981532402477848,0.309148743347499,1.252800181962738,0.098829881195179,0.74888084504147,0.181494956495445,1.304931340021904,0.595939453660087]
[mean = 0.99223677282248,stdev = 0.510336641737154]

PicoLisp

Translation of: C

<lang PicoLisp>(load "@lib/math.l")

(de randomNormal () # Normal distribution, centered on 0, std dev 1

  (*/
     (sqrt (* -2.0 (log (rand 0 1.0))))
     (cos (*/ 2.0 pi (rand 0 1.0) `(* 1.0 1.0)))
     1.0 ) )

(seed (time)) # Randomize

(let Result

  (make                                           # Build list
     (do 1000                                     # of 1000 elements
        (link (+ 1.0 (/ (randomNormal) 2))) ) )
  (for N (head 7 Result)                          # Print first 7 results
     (prin (format N *Scl) " ") ) )</lang>

Output:

1.500334 1.212931 1.095283 0.433122 0.459116 1.302446 0.402477

PL/I

<lang PL/I> /* CONVERTED FROM WIKI FORTRAN */ Normal_Random: procedure options (main);

  declare (array(1000), pi, temp,
           mean initial (1.0), sd initial (0.5)) float (18);
  declare (i, n) fixed binary;
  
  n = hbound(array, 1);
  pi = 4.0*ATAN(1.0);
  array = random(); /* Uniform distribution */
  /* Now convert to normal distribution */
  DO i = 1 to n-1 by 2;
     temp = sd * SQRT(-2.0*LOG(array(i))) * COS(2*pi*array(i+1)) + mean;
     array(i+1) = sd * SQRT(-2.0*LOG(array(i))) * SIN(2*pi*array(i+1)) + mean;
     array(i) = temp;
  END;
  /* Check mean and standard deviation */
  mean = SUM(array)/n;
  sd = SQRT(SUM((array - mean)**2)/n);
  put skip edit ( "Mean = ", mean ) (a, F(18,16) );
  put skip edit ( "Standard Deviation = ", sd) (a, F(18,16));

END Normal_Random; </lang>

Output:

Mean = 1.0125630677913652  Standard Deviation = 0.5067289784535284
3 runs with different seeds to random():
Mean = 1.0008390411168471  Standard Deviation = 0.5095810511317908
Mean = 0.9754351286894838  Standard Deviation = 0.4804376530558166
Mean = 1.0177411222687990  Standard Deviation = 0.5165899662493400

PL/SQL

<lang PL/SQL> DECLARE

 --The desired collection
 type t_coll is table of number index by binary_integer; 
 l_coll t_coll;

 c_max pls_integer := 1000;

BEGIN

  FOR l_counter IN 1 .. c_max
  LOOP
     -- dbms_random.normal delivers normal distributed random numbers with a mean of 0 and a variance of 1
     -- We just adjust the values and get the desired result:
     l_coll(l_counter) := DBMS_RANDOM.normal * 0.5 + 1;
     DBMS_OUTPUT.put_line (l_coll(l_counter));
  END LOOP;

END; </lang>

Pop11

<lang pop11>;;; Choose radians as arguments to trigonometic functions true -> popradians;

define random_normal() -> result; lvars r1 = random0(1.0), r2 = random0(1.0);

    cos(2*pi*r1)*sqrt(-2*log(r2)) -> result

enddefine;

lvars array, i;

for i from 1 to 1000 do 1.0+0.5*random_normal() endfor;

consvector(1000) -> array;</lang>

PowerShell

Equation adapted from Liberty BASIC <lang powershell>function Get-RandomNormal

   {
   [CmdletBinding()]
   Param ( [double]$Mean, [double]$StandardDeviation )

   $RandomNormal = $Mean + $StandardDeviation * [math]::Sqrt( -2 * [math]::Log( ( Get-Random -Minimum 0.0 -Maximum 1.0 ) ) ) * [math]::Cos( 2 * [math]::PI * ( Get-Random -Minimum 0.0 -Maximum 1.0 ) )

   return $RandomNormal
   }

Standard deviation function for testing

function Get-StandardDeviation

   {
   [CmdletBinding()]
   param ( [double[]]$Numbers )

   $Measure = $Numbers | Measure-Object -Average
   $PopulationDeviation = 0
   ForEach ($Number in $Numbers) { $PopulationDeviation += [math]::Pow( ( $Number - $Measure.Average ), 2 ) }
   $StandardDeviation = [math]::Sqrt( $PopulationDeviation / ( $Measure.Count - 1 ) )
   return $StandardDeviation
   }

Test

$RandomNormalNumbers = 1..1000 | ForEach { Get-RandomNormal -Mean 1 -StandardDeviation 0.5 }

$Measure = $RandomNormalNumbers | Measure-Object -Average

$Stats = [PSCustomObject]@{

   Count             = $Measure.Count
   Average           = $Measure.Average
   StandardDeviation = Get-StandardDeviation -Numbers $RandomNormalNumbers

}

$Stats | Format-List </lang>

Output:

Count             : 1000
Average           : 1.01206560135809
StandardDeviation : 0.489099623426272

PureBasic

<lang PureBasic>Procedure.f RandomNormal()

  ; This procedure can return any real number.
  Protected.f x1, x2

  ; random numbers from the open interval ]0, 1[
  x1 = (Random(999998)+1) / 1000000       ; must be > 0 because of Log(x1)
  x2 = (Random(999998)+1) / 1000000

  ProcedureReturn Sqr(-2*Log(x1)) * Cos(2*#PI*x2)

EndProcedure

Define i, n=1000

Dim a.q(n-1) For i = 0 To n-1

  a(i) = 1 + 0.5 * RandomNormal()

Next</lang>

Python

Using random.gauss

<lang python>>>> import random >>> values = [random.gauss(1, .5) for i in range(1000)] >>> </lang>

Quick check of distribution

<lang python>>>> def quick_check(numbers):

   count = len(numbers)
   mean = sum(numbers) / count
   sdeviation = (sum((i - mean)**2 for i in numbers) / count)**0.5
   return mean, sdeviation

>>> quick_check(values) (1.0140373306786599, 0.49943411329234066) >>> </lang>

Note that the random module in the Python standard library supports a number of statistical distribution methods.

Alternatively using random.normalvariate

<lang python>>>> values = [ random.normalvariate(1, 0.5) for i in range(1000)] >>> quick_check(values) (0.990099111944864, 0.5029847005836282) >>> </lang>

R

<lang r># For reproducibility, set the seed: set.seed(12345L)

result <- rnorm(1000, mean = 1, sd = 0.5)</lang>

Racket

lang racket

(for/list ([i 1000])

 (add1 (* (sqrt (* -2 (log (random)))) (cos (* 2 pi (random))) 0.5)))

</lang>

Alternative: <lang Racket>

lang racket

(require math/distributions) (sample (normal-dist 1.0 0.5) 1000) </lang>

Raku

(formerly Perl 6)

Works with: Rakudo version #22 "Thousand Oaks"

<lang perl6>sub randnorm ($mean, $stddev) {

   $mean + $stddev * sqrt(-2 * log rand) * cos(2 * pi * rand)

}

my @nums = randnorm(1, 0.5) xx 1000;

Checking

say my $mean = @nums R/ [+] @nums; say my $stddev = sqrt $mean**2 R- @nums R/ [+] @nums X** 2; </lang>

Raven

<lang raven>define PI

  -1 acos

define rand1

  9999999 choose 1 + 10000000.0 /

define randNormal

  rand1 PI * 2 * cos
  rand1 log -2 * sqrt
  *
  2 / 1 +

1000 each drop randNormal "%f\n" print</lang> Quick Check (on linux with code in file rand.rv) <lang raven>raven rand.rv | awk '{sum+=$1; sumsq+=$1*$1;} END {print "stdev = " sqrt(sumsq/NR - (sum/NR)**2); print "mean = " sum/NR}' stdev = 0.497773 mean = 1.01497</lang>

REXX

The REXX language doesn't have any "higher math" functions like SQRT/SIN/COS/LN/LOG/EXP/POW/etc.,
so we hoi polloi REXX programmers have to roll our own.

Programming note: note the range of the random numbers: (0,1]
(that is, random numbers from zero──►unity, excluding zero, including unity). <lang rexx>/*REXX pgm generates 1,000 normally distributed numbers: mean=1, standard deviation=½.*/ numeric digits 20 /*the default decimal digit precision=9*/ parse arg n seed . /*allow specification of N and the seed*/ if n== | n=="," then n=1000 /*N: is the size of the array. */ if datatype(seed,'W') then call random ,,seed /*SEED: for repeatable random numbers. */ newMean=1 /*the desired new mean (arithmetic avg)*/ sd=1/2 /*the desired new standard deviation. */

      do g=1  for n                             /*generate  N uniform random #'s (0,1].*/
      #.g = random(1, 1e5)  /  1e5              /*REXX's RANDOM BIF generates integers.*/
      end   /*g*/                               /* [↑]  random integers ──► fractions. */

say ' old mean=' mean() say 'old standard deviation=' stdDev() call pi; pi2=pi * 2 /*define pi and also 2 * pi. */ say

      do j=1  to n-1  by 2;    m=j+1            /*step through the iterations by two.  */
          _=sd *  sqrt(ln(#.j) * -2)            /*calculate the  used-twice expression.*/
      #.j=_ * cos(pi2 * #.m)  +  newMean        /*utilize the  Box─Muller method.      */
      #.m=_ * sin(pi2 * #.m)  +  newMean        /*random number must be:      (0,1]    */
      end   /*j*/

say ' new mean=' mean() say 'new standard deviation=' stdDev() exit /*stick a fork in it, we're all done. */ /*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ mean: _=0; do k=1 for n; _=_ + #.k; end; return _/n stdDev: _avg=mean(); _=0; do k=1 for n; _=_ + (#.k - _avg)**2; end; return sqrt(_/n) e: e =2.7182818284590452353602874713526624977572470936999595749669676277240766303535; return e /*digs overkill*/ pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862; return pi /* " " */ r2r: return arg(1) // (pi() * 2) /*normalize ang*/ sin: procedure; parse arg x;x=r2r(x);numeric fuzz min(5,digits()-3);if abs(x)=pi then return 0;return .sincos(x,x,1) .sincos:parse arg z,_,i; x=x*x; p=z; do k=2 by 2; _=-_*x/(k*(k+i)); z=z+_; if z=p then leave; p=z; end; return z /*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ ln: procedure; parse arg x,f; call e; ig= x>1.5; is=1 - 2 * (ig\==1); ii=0; xx=x

         do while ig&xx>1.5|\ig&xx<.5;_=e;do k=-1;iz=xx*_**-is;if k>=0&(ig&iz<1|\ig&iz>.5) then leave;_=_*_;izz=iz;end
       xx=izz;ii=ii+is*2**k;end;x=x*e**-ii-1;z=0;_=-1;p=z;do k=1;_=-_*x;z=z+_/k;if z=p then leave;p=z;end; return z+ii

/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ cos: procedure; parse arg x; x=r2r(x); a=abs(x); hpi=pi * .5

           numeric fuzz min(6, digits() - 3);      if a=pi    then return -1
           if a=hpi | a=hpi*3  then return 0;      if a=pi/3  then return .5
           if a=pi * 2/3       then return -.5;                    return .sinCos(1,1,-1)

/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits; h=d+6

       numeric form; parse value format(x,2,1,,0) 'E0'  with  g 'E' _ .;  g=g * .5'e'_ %2
       m.=9;     do j=0  while h>9;       m.j=h;                h=h%2 + 1;      end /*j*/
                 do k=j+5  to 0  by -1;   numeric digits m.k;   g=(g+x/g)*.5;   end /*k*/
       numeric digits d;     return g/1</lang>

output when using the default inputs:

              old mean= 0.5015724
old standard deviation= 0.28652466389342471402

              new mean= 0.98807025356443262689
new standard deviation= 0.50002924192766720838

Ring

<lang ring> for i = 1 to 10

   see random(i) + nl

next i </lang>

Ruby

<lang ruby>Array.new(1000) { 1 + Math.sqrt(-2 * Math.log(rand)) * Math.cos(2 * Math::PI * rand) }</lang>

Run BASIC

<lang runbasic>dim a(1000) pi = 22/7 for i = 1 to 1000

  a( i)  = 1 + .5 * (sqr(-2 * log(rnd(0))) * cos(2 * pi * rnd(0)))

next i</lang>

Rust

Library: rand

Using a for-loop: <lang rust>extern crate rand; use rand::distributions::{Normal, IndependentSample};

fn main() {

   let mut rands = [0.0; 1000];
   let normal = Normal::new(1.0, 0.5);
   let mut rng = rand::thread_rng();
   for num in rands.iter_mut() {
       *num = normal.ind_sample(&mut rng);
   }

}</lang>

Using iterators: <lang rust>extern crate rand; use rand::distributions::{Normal, IndependentSample};

fn main() {

   let rands: Vec<_> = {
       let normal = Normal::new(1.0, 0.5);
       let mut rng = rand::thread_rng();
       (0..1000).map(|_| normal.ind_sample(&mut rng)).collect()
   };

}</lang>

SAS

<lang SAS> /* Generate 1000 random numbers with mean 1 and standard deviation 0.5.

 SAS version 9.2 was used to create this code.*/

data norm1000;

 call streaminit(123456);

/* Set the starting point, so we can replicate results.

  If you want different results each time, comment the above line. */
 do i=1 to 1000;
   r=rand('normal',1,0.5);
   output;
 end;

run; </lang> Results:

 The MEANS Procedure

                     Analysis Variable : r

                          Mean         Std Dev
                  ----------------------------
                     0.9907408       0.4844051
                  ----------------------------

Sather

<lang sather>class MAIN is

 main is
   a:ARRAY{FLTD} := #(1000);
   i:INT;

   RND::seed(2010);
   loop i := 1.upto!(1000) - 1;
     a[i] := 1.0d + 0.5d * RND::standard_normal;
   end;

   -- testing the distribution
   mean ::= a.reduce(bind(_.plus(_))) / a.size.fltd;
   #OUT + "mean " + mean + "\n";
   a.map(bind(_.minus(mean)));
   a.map(bind(_.pow(2.0d)));
   dev ::= (a.reduce(bind(_.plus(_))) / a.size.fltd).sqrt;
   #OUT + "dev  " + dev + "\n";
 end;

end;</lang>

Scala

One liner

<lang scala>List.fill(1000)(1.0 + 0.5 * scala.util.Random.nextGaussian)</lang>

Academic

<lang scala> object RandomNumbers extends App {

 val distribution: LazyList[Double] = {
   def randomNormal: Double = 1.0 + 0.5 * scala.util.Random.nextGaussian

   def normalDistribution(a: Double): LazyList[Double] = a #:: normalDistribution(randomNormal)

   normalDistribution(randomNormal)
 }

 /*
  * Let's test it
  */
 def calcAvgAndStddev[T](ts: Iterable[T])(implicit num: Fractional[T]): (T, Double) = {
   val mean: T =
     num.div(ts.sum, num.fromInt(ts.size)) // Leaving with type of function T

   // Root of mean diffs
   val stdDev = Math.sqrt(ts.map { x =>
     val diff = num.toDouble(num.minus(x, mean))
     diff * diff
   }.sum / ts.size)

   (mean, stdDev)
 }

 println(calcAvgAndStddev(distribution.take(1000))) // e.g. (1.0061433267806525,0.5291834867560893)

} </lang>

Scheme

<lang scheme>; linear congruential generator given in C99 section 7.20.2.1 (define ((c-rand seed)) (set! seed (remainder (+ (* 1103515245 seed) 12345) 2147483648)) (quotient seed 65536))

uniform real numbers in open interval (0, 1)

(define (unif-rand seed) (let ((r (c-rand seed))) (lambda () (/ (+ (r) 1) 32769.0))))

Box-Muller method to generate normal distribution

(define (normal-rand unif m s) (let ((? #t) (! 0.0) (twopi (* 2.0 (acos -1.0)))) (lambda ()

  (set! ? (not ?))
  (if ? !
        (let ((a (sqrt (* -2.0 (log (unif))))) (b (* twopi (unif))))
             (set! ! (+ m (* s a (sin b))))
             (+ m (* s a (cos b))))))))

(define rnorm (normal-rand (unif-rand 0) 1.0 0.5))

auxiliary function to get a list of 'n random numbers from generator 'r

(define (rand-list r n) = (if (zero? n) '() (cons (r) (rand-list r (- n 1)))))

(define v (rand-list rnorm 1000))

v

|

(-0.27965824722565835

-0.8870860825789542
0.6499618744638194
0.31336141955110863
...
0.5648743998193049
0.8282656735558756
0.6399951934564637
0.7699535302478072)

|#

check mean and standard deviation

(define (mean-sdev v) (let loop ((v v) (a 0) (b 0) (n 0)) (if (null? v)

   (let ((mean (/ a n)))
        (list mean (sqrt (/ (- b (* n mean mean)) (- n 1)))))
   (let ((x (car v)))
        (loop (cdr v) (+ a x) (+ b (* x x)) (+ n 1))))))

(mean-sdev v)

(0.9562156817697293 0.5097087109575911)</lang>

Seed7

<lang seed7>$ include "seed7_05.s7i";

 include "float.s7i";
 include "math.s7i";

const func float: frand is func # Uniform distribution, (0..1]

 result
   var float: frand is 0.0;
 begin
   repeat
     frand := rand(0.0, 1.0);
   until frand <> 0.0;
 end func;

const func float: randomNormal is # Normal distribution, centered on 0, std dev 1

 return sqrt(-2.0 * log(frand)) * cos(2.0 * PI * frand);

const proc: main is func

 local
   var integer: i is 0;
   var array float: rands is 1000 times 0.0;
 begin
   for i range 1 to length(rands) do
     rands[i] := 1.0 + 0.5 * randomNormal;
   end for;
 end func;</lang>

Sidef

Standard ML

Works with: SML/NJ

SML/NJ has two structures for random numbers:

1) Rand (a linear congruential generator). You create the generator by calling Rand.mkRandom with a seed (of word type). You can call the generator with () repeatedly to get a word in the range [Rand.randMin, Rand.randMax]. You can use the Rand.norm function to transform the output into a real from 0 to 1, or use the Rand.range (i,j) function to transform the output into an int of the given range. <lang sml>val seed = 0w42; val gen = Rand.mkRandom seed; fun random_gaussian () =

 1.0 + Math.sqrt (~2.0 * Math.ln (Rand.norm (gen ()))) * Math.cos (2.0 * Math.pi * Rand.norm (gen ()));

val a = List.tabulate (1000, fn _ => random_gaussian ());</lang>

2) Random (a subtract-with-borrow generator). You create the generator by calling Random.rand with a seed (of a pair of ints). You can use the Random.randInt function to generate a random int over its whole range; Random.randNat to generate a non-negative random int; Random.randReal to generate a real between 0 and 1; or Random.randRange (i,j) to generate an int in the given range. <lang sml>val seed = (47,42); val gen = Random.rand seed; fun random_gaussian () =

 1.0 + Math.sqrt (~2.0 * Math.ln (Random.randReal gen)) * Math.cos (2.0 * Math.pi * Random.randReal gen);

val a = List.tabulate (1000, fn _ => random_gaussian ());</lang>

Other implementations of Standard ML have their own random number generators. For example, Moscow ML has a Random structure that is different from the one from SML/NJ.

Works with: PolyML

The SML Basis Library does not provide a routine for uniform deviate generation, and PolyML does not have one. Using a routine from "Monte Carlo" by Fishman (Springer), in the function uniformdeviate, and avoiding the slow IntInf's: <lang smlh> val urandomlist = fn seed => fn n => let val uniformdeviate = fn seed => let val in31m = (Real.fromInt o Int32.toInt ) (getOpt (Int32.maxInt,0) ); val in31 = in31m +1.0; val s1 = 41160.0; val s2 = 950665216.0; val v = Real.realFloor seed; val val1 = v*s1; val val2 = v*s2; val next1 = Real.fromLargeInt (Real.toLargeInt IEEEReal.TO_NEGINF (val1/in31)) ; val next2 = Real.rem(Real.realFloor(val2/in31) , in31m ); val valt = val1+val2 - (next1+next2)*in31m; val nextt = Real.realFloor(valt/in31m); val valt = valt - nextt*in31m; in (valt/in31m,valt) end; val store = ref (0.0,0.0); val rec u = fn S => fn 0 => [] | n=> (store:=uniformdeviate S; (#1 (!store)):: (u (#2 (!store)) (n-1))) ; in u seed n end;

local open Math in val bmconv = fn urand => fn vrand => 1.0+0.5*(sqrt(~2.0*ln urand)*cos (2.0*pi*vrand) ) end;

val rec makeNormals = fn once => fn u::v::[] => [once u v] | u::v::rm => (once u v )::(makeNormals once rm );

val anyrealseed=1009.0 ; makeNormals bmconv (urandomlist anyrealseed 2000); </lang>

Stata

<lang stata>clear all set obs 1000 gen x=rnormal(1,0.5)</lang>

Mata

<lang stata>a = rnormal(1000,1,1,0.5)</lang>

Tcl

<lang tcl>package require Tcl 8.5 variable ::pi [expr acos(0)] proc ::tcl::mathfunc::nrand {} {

   expr {sqrt(-2*log(rand())) * cos(2*$::pi*rand())}

}

set mean 1.0 set stddev 0.5 for {set i 0} {$i < 1000} {incr i} {

   lappend result [expr {$mean + $stddev*nrand()}]

}</lang>

TI-83 BASIC

Builtin function: randNorm()

 randNorm(1,.5)

Or by a program:

Calculator symbol translations:

"STO" arrow: →

Square root sign: √

ClrList L₁
Radian
For(A,1,1000)
√(-2*ln(rand))*cos(2*π*A)→L₁(A)
End

TorqueScript

Ursala

There are two ways of interpreting the task, either to simulate sampling a population described by the given statistics, or to construct a sample exhibiting the given statistics. Both are illustrated below. The functions parameterized by the mean and standard deviation take a sample size and return a sample of that size, represented as a list of floating point numbers. The Z library function simulates a draw from a standard normal distribution. Mean and standard deviation library functions are also used in this example. <lang Ursala>#import nat

import flo

pop_stats("mu","sigma") = plus/*"mu"+ times/*"sigma"+ Z*+ iota

sample_stats("mu","sigma") = plus^*D(minus/"mu"+ mean,~&)+ vid^*D(div\"sigma"+ stdev,~&)+ Z*+ iota

cast %eWL

test =

^(mean,stdev)* <

  pop_stats(1.,0.5) 1000,
  sample_stats(1.,0.5) 1000></lang>

The output shows the mean and standard deviation for both sample vectors, the latter being exact by construction.

<
   (1.004504e+00,4.915525e-01),
   (1.000000e+00,5.000000e-01)>

Visual FoxPro

<lang vfp> LOCAL i As Integer, m As Double, n As Integer, sd As Double py = PI() SET TALK OFF SET DECIMALS TO 6 CREATE CURSOR gdev (deviate B(6)) RAND(-1) n = 1000 m = 1 sd = 0.5 CLEAR FOR i = 1 TO n INSERT INTO gdev VALUES (GaussDev(m, 1/sd)) ENDFOR CALCULATE AVG(deviate), STD(deviate) TO m, sd ? "Mean", m, "Std Dev", sd SET TALK ON SET DECIMALS TO

FUNCTION GaussDev(mean As Double, sdev As Double) As Double LOCAL z As Double z = SQRT(-2*LOG(RAND()))*COS(py*RAND()) IF sdev # 0 z = mean + z/sdev ENDIF RETURN z ENDFUNC </lang>

Wren

<lang ecmascript>import "random" for Random

var rand = Random.new()

var randNormal = Fn.new { (-2 * rand.float().log).sqrt * (2 * Num.pi * rand.float()).cos }

var stdDev = Fn.new { |a, m|

   var c = a.count
   return ((a.reduce(0) { |acc, x| acc + x*x } - m*m*c) / c).sqrt

}

var n = 1000 var numbers = List.filled(n, 0) var mu = 1 var sigma = 0.5 var sum = 0 for (i in 0...n) {

   numbers[i] = mu + sigma*randNormal.call()
   sum = sum + numbers[i]

} var mean = sum / n System.print("Actual mean : %(mean)") System.print("Actual std dev: %(stdDev.call(numbers, mean))")</lang>

Output:

Sample run:

Actual mean   : 1.0053988699746
Actual std dev: 0.4961645117026

Yorick

Returns array of count random numbers with mean 0 and standard deviation 1. <lang yorick>func random_normal(count) {

  return sqrt(-2*log(random(count))) * cos(2*pi*random(count));

}</lang>

Example of basic use:

> nums = random_normal(1000); // create an array 1000 random numbers
> nums(avg); // show the mean
0.00901216
> nums(rms); // show the standard deviation
0.990265

Example with a mean of 1.0 and a standard deviation of 0.5:

> nums = random_normal(1000) * 0.5 + 1;
> nums(avg);
1.00612
> nums(rms);
0.496853

zkl

<lang zkl>fcn mkRand(mean,sd){ //normally distributed random w/mean & standard deviation

  pi:=(0.0).pi;    // using the Box–Muller transform
  rz1:=fcn{1.0-(0.0).random(1)}  // from [0,1) to (0,1]
  return('wrap(){((-2.0*rz1().log()).sqrt() * (2.0*pi*rz1()).cos())*sd + mean })

}</lang> This creates a new random number generator, now to use it: <lang zkl>var g=mkRand(1,0.5); ns:=(0).pump(1000,List,g); // 1000 rands with mean==1 & sd==1/2 mean:=(ns.sum(0.0)/1000); //-->1.00379

  // calc sd of list of numbers:

(ns.reduce('wrap(p,n){p+(n-mean).pow(2)},0.0)/1000).sqrt() //-->0.494844</lang>

ZX Spectrum Basic

Here we have converted the QBasic code to suit the ZX Spectrum:

<lang zxbasic>10 RANDOMIZE 0 : REM seeds random number generator based on uptime 20 DIM a(1000) 30 CLS 40 FOR i = 1 TO 1000 50 LET a(i) = 1 + SQR(-2 * LN(RND)) * COS(2 * PI * RND) 60 NEXT i</lang>