Random numbers

From Rosetta 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 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 ... )

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


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>

Tcl

<lang tcl>package require Tcl 8.5 proc ::tcl::mathfunc::nrand {} {return [expr {sqrt(-2*log(rand()))*cos(4*acos(0)*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=28720"