Random numbers

From Rosetta Code
Revision as of 05:09, 20 June 2009 by rosettacode>Tinku99 (+ AutoHotkey)
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 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
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 ... )

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

C

<lang c>#include <stdlib.h>

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

}</lang>

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

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

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

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

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

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

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

}</lang>

JavaScript

<lang 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;</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. 

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
)

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

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>

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


Perl

Library: Math::CephesCephes

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

map { 

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

} 1..1000</lang>

PHP

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

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

Python

Works with: Python version 2.5

<lang python>import random randList = [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.

R

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

Ruby

<lang ruby>Array.new(1000) { 1 + Math.sqrt(-2 * Math.log(rand)) * Math.cos(2 * Math::PI * rand) }</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
Retrieved from "https://rosettacode.org/wiki/Random_numbers?oldid=28728"