Random numbers

From Rosetta Code
Revision as of 20:08, 20 November 2012 by rosettacode>Gerard Schildberger (→‎{{header|REXX}}: changed pronoun (in section header comments). -- ~~~~)
Task
Random numbers
You are encouraged to solve this task according to the task description, using any language you may know.

The goal of this task is to generate a collection filled with 1000 normally distributed random (or pseudorandom) numbers with a mean of 1.0 and a standard deviation of 0.5

Many libraries only generate uniformly distributed random numbers. If so, use this formula to convert them to a normal distribution.

Ada

<lang ada>with Ada.Numerics; use Ada.Numerics; with Ada.Numerics.Float_Random; use Ada.Numerics.Float_Random; with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;

procedure Normal_Random is 

  function Normal_Distribution
           (  Seed  : Generator;
              Mu    : Float := 1.0;
              Sigma : Float := 0.5
           )  return Float is 
  begin
     return
        Mu + (Sigma * Sqrt (-2.0 * Log (Random (Seed), 10.0)) * Cos (2.0 * Pi * Random (Seed)));
  end Normal_Distribution;
     
  Seed         : Generator;
  Distribution : array (1..1_000) of Float;

begin 

  Reset (Seed);
  for I in Distribution'Range loop
     Distribution (I) := Normal_Distribution (Seed);
  end loop;

end Normal_Random;</lang>

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

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

AutoHotkey

contributed by Laszlo on the ahk forum <lang AutoHotkey>Loop 40 

  R .= RandN(1,0.5) "`n"  ; mean = 1.0, standard deviation = 0.5

MsgBox %R%

RandN(m,s) { ; Normally distributed random numbers of mean = m, std.dev = s by Box-Muller method 

  Static i, Y
  If (i := !i) { ; every other call
     Random U, 0, 1.0
     Random V, 0, 6.2831853071795862
     U := sqrt(-2*ln(U))*s
     Y := m + U*sin(V)
     Return m + U*cos(V)
  }
  Return Y

}</lang>

AWK

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

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

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>

  1. include <math.h>
  2. ifndef M_PI
  3. define M_PI 3.14159265358979323846
  4. 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++

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

  1. include <cmath> // for atan, sqrt, log, cos
  2. 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.

<lang cpp>

  1. include <vector>
  2. include "boost/random.hpp"
  3. include "boost/generator_iterator.hpp"
  4. include <boost/random/normal_distribution.hpp>
  5. 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>

Works with: C++11

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

<lang cpp>#include <random>

  1. include <functional>
  2. include <vector>
  3. include <algorithm>

using namespace std;

int main() { 

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

 vector<double> v(1000);
 generate(v.begin(), v.end(), rnd);
 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>

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>

D

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

struct NormalRandom { 

   double mean, stdDev;

   // needed 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 double r1 = uniform(0.0, 1.0);
       immutable double r2 = uniform(0.0, 1.0);
       return mean + stdDev * sqrt(-2 * log(r1)) * cos(2 * PI * r2);
   }

}

void main() { 

   double[1000] array;
   auto nrnd = NormalRandom(1.0, 0.5);
   foreach (ref x; array)
       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 := RandG(1, 0.5);</lang>

E

<lang e>accum [] for _ in 1..1000 { _.with(entropy.nextGaussian()) }</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

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

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>

Factor

<lang factor>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>

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

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>

Go

<lang go>package main

import ( 

   "math/rand"
   "time"

)

const mean = 1.0 const stdv = .5

func main() { 

   var list [1000]float64
   rand.Seed(time.Now().UnixNano())
   for i := range list {
       list[i] = mean + stdv*rand.NormFloat64()
   }

}</lang>

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>

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>

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>

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>

Mathematica

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

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>

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>

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>

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


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*u1) \\ Could easily be extended with a second normal at very little cost. }; vector(1000,unused,rnormal()/2+1)</lang>

Pascal

See Delphi

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>

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 = map { randnorm 1, 0.5 }, ^1000;</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>

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>create or replace PROCEDURE PROCEDURE1 AS 

   TYPE numsColl is TABLE OF NUMBER;
   nums numsColl;

   FUNCTION GenNums(n IN NUMBER) RETURN numsColl AS 
       PI NUMBER := ACOS (-1);
   BEGIN
       nums := numsColl();
       nums.extend(n);

       FOR i in 1 .. n LOOP
           nums(i) := 1 + .5 * (sqrt(-2 * log(dbms_random.value, 10)) * cos(2 * PI * dbms_random.value));
       END LOOP;
       
       RETURN nums;
   END GenNums;

BEGIN 

   nums := GenNums(10);
   FOR i in 1 .. 10 LOOP
       DBMS_OUTPUT.PUT_LINE(nums(i));
   END LOOP;

END PROCEDURE1;</lang>

Pop11

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

procedure generating standard normal distribution

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;

Put numbers on the stack

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

collect them into array

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

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

Works with: Python version 2.5

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

  1. or [ random.normalvariate(1, 0.5) for i in range(1000)]</lang>

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

Quick check <lang python>>>> mean = sum(values)/1000 >>> sdeviation = (sum((i - mean)**2 for i in values)/1000)**0.5 >>> mean, sdeviation (1.0127861555468178, 0.5006682783828207) </lang>

R

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

REXX

The REXX language doesn't have any "higher math" functions like SIN/COS/LN/LOG/SQRT/POW/etc,
so we hoi polloi programmers have roll our own. <lang rexx>/*REXX pgm gens 1,000 normally distributed #s: mean=1, standard dev.=½. */

parse arg n .; if n== then n=1000 /* N is the size of the array. */ call pi /*call subroutine to define pi. */ mean=1 /*desired new mean (arith. avg.) */ sd=1/2 /*desired new standard deviation.*/ 

 do g=1 for n                         /*generate  N  uniform random #s.*/
 #.g=random(0,1e5)/1e5                /*REXX gens uniform rand integers*/
 end   /*g*/

say ' old mean=' mean() say 'old standard deviation=' stddev() say 

     do j=1 to n-1 by 2
     m=j+1
       _=sd*sqrt(-2*ln(#.j))*cos(2*pi*#.m)+mean /*use Box-Muller method*/
     #.m=sd*sqrt(-2*ln(#.j))*sin(2*pi*#.m)+mean /*rand # must be 0──>1.*/
     #.j=_
     end    /*j*/

say ' new mean=' mean() say 'new standard deviation=' stddev() exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────subroutines─────────────────────────*/ 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 pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862; return pi r2r: return arg(1)//(2*pi()) sqrt: procedure;parse arg x; if x=0 then return 0; d=digits(); numeric digits 11; g=.sqrtGuess() 

       do j=0 while p>9; m.j=p; p=p%2+1; end; do k=j+5 to 0 by -1; if m.k>11 then numeric digits m.k
       g=.5*(g+x/g); end; numeric digits d; return g/1

.sqrtGuess: numeric form; m.=11; p=d+d%4+2 

         parse value format(x,2,1,,0) 'E0' with g 'E' _ .; return g*.5'E'_%2

cos: procedure; arg x; x=r2r(x); a=abs(x); numeric fuzz min(9,digits()-9); if a=pi() then return -1 

       if a=pi()/2|a=2*pi() then return 0;if a=pi()/3 then return .5;if a=2*pi()/3 then return -.5;return .sincos(1,1,-1)

sin: procedure; 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; return .ln_comp() .ln_comp: do while ig&xx>1.5|\ig&xx<.5;_=e;do k=-1;iz=xx*_**-is;if k>=0&(ig&iz<1|\ig&iz>.5) then leave;_=_*_;izz=iz;end 

         xx=izz;ii=ii+is*2**k;end;x=x*e**-ii-1;z=0;_=-1;p=z;do k=1;_=-_*x;z=z+_/k;if z=p then leave;p=z;end;return z+ii</lang>

output when using the default input 

              old mean= 0.50398978
old standard deviation= 0.292088074

              new mean= 0.999477923
new standard deviation= 0.505691743

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>

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

<lang scala>List.fill(1000)(1.0 + 0.5 * scala.util.Random.nextGaussian)</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

  1. |

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

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.

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

Calculator symbol translations:

"STO" arrow: →

Square root sign: √ 

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

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

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

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

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

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>

Retrieved from "https://rosettacode.org/wiki/Random_numbers?oldid=28832"