Random numbers

The goal of this task is to generate a collection filled with 1000 normally distributed random 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

<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; </Ada>

ALGOL 68

Translation of: C

PROC random normal = REAL:  # normal distribution, centered on 0, std dev 1 #
(
  sqrt(-2*log(random)) * cos(2*pi*random)
);
main:
(
  [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]))
)

Output sample:

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

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

C

<c>#include <stdlib.h>

include <math.h>

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;

}</c>

C++

<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() // 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:

 double mu, sigma;

};

int main() {

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

}</cpp>

Common Lisp

<lisp>(loop for i from 1 to 1000

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

E

accum [] for _ in 1..1000 { _.with(entropy.nextGaussian()) }

Forth

Works with: gforth version 0.6.2

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

Fortran

Works with: Fortran version 90 and later

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

Output

Mean = 0.995112
Standard Deviation = 0.503373

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

Groovy

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

IDL

result = 1.0 + 0.5*randomn(seed,1000)

J

urand=: ?@$ 0: 
zrand=: (2 o. 2p1 * urand) * [: %: _2 * [: ^. urand

1 + 0.5 * zrand 100

Java

<java>double[] list = new double[1000]; Random rng = new Random(); for(int i = 0;i<list.length;i++) {

 list[i] = 1.0 + 0.5 * rng.nextGaussian()

}</java>

JavaScript

<javascript>function randomNormal() {

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

}

var a = new Array(1000); for (var i=0; i<a.length; i++)

 a[i] = randomNormal() / 2 + 1;</javascript>

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.

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 ?] []

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
)

OCaml

<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 ());;</ocaml>

Perl

Library: Math::Cephes Cephes

<perl>use Math::Cephes qw($PI);

map {

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

} 1..1000</perl>

PHP

<php>$pi = pi(); // Set PI $a = range(1,1000); // Create array // Cycle array values foreach(range(1,1000) as $i){

     $a[$i] = 1 + sqrt(-2 * log(mt_rand())) * cos(2 * $pi * mt_rand());

}</php>

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;

Python

Works with: Python version 2.5

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

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

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

R

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

Ruby

(1..1000).map { 1 + Math.sqrt(-2 * Math.log(rand)) * Math.cos(2 * Math::PI * rand) }

Tcl

proc nrand {} {return [expr sqrt(-2*log(rand()))*cos(4*acos(0)*rand())]}
for {set i 0} {$i < 1000} {incr i} {lappend result [expr 1+.5*nrand()]}

TI-83 BASIC

Calculator symbol translations:

"STO" arrow: →

Square root sign: √

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