Random numbers

From Rosetta Code
Revision as of 21:41, 11 July 2013 by Grondilu (talk | contribs) (→‎{{header|Perl 6}}: adding checking code)
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++

Works with: C++11

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

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

Works with: C++03

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

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

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>

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>

МК-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.

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 = randnorm(1, 0.5) xx 1000;

  1. Checking

say my $mean = @nums R/ [+] @nums; say my $stddev = sqrt $mean**2 R- @nums R/ [+] @nums X** 2; </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>

Racket

<lang Racket>

  1. lang racket

(for/list ([i 1000]) 

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

</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 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.=½. */ call pi /*call subroutine to define 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.*/
         #.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>

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

<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=28844"