Random numbers

From Rosetta Code
Jump to: navigation, search
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.

Contents

[edit] 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;

[edit] 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
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]))
)

Output sample:

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

[edit] AutoHotkey

contributed by Laszlo on the ahk forum

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
}

[edit] 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}'

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

[edit] BBC BASIC

      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

Output:

Mean = 1.01848064
Standard deviation = 0.503551814

[edit] C

#include <stdlib.h>
#include <math.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#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;
}

[edit] C#

Translation of: JavaScript
 
private static double randomNormal()
{
return Math.Cos(2 * Math.PI * tRand.NextDouble()) * Math.Sqrt(-2 * Math.Log(tRand.NextDouble()));
}
 

Then the methods in Random numbers#Metafont are used to calculate the average and the Standard Deviation:

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

An example result:

Average: 1,00510073053613
Standard Deviation: 0,502540443430955

[edit] C++

Works with: C++11

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

#include <random>
#include <functional>
#include <vector>
#include <algorithm>
using namespace std;
 
int main()
{
random_device seed;
mt19937 engine(seed());
normal_distribution<double> dist(1.0, 0.5);
auto rnd = bind(dist, engine);
 
vector<double> v(1000);
generate(v.begin(), v.end(), rnd);
return 0;
}
Works with: C++03
#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() 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;
}
Library: Boost

This example used Mersenne Twister generator. It can be changed by changing the typedef.

 
#include <vector>
#include "boost/random.hpp"
#include "boost/generator_iterator.hpp"
#include <boost/random/normal_distribution.hpp>
#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;
}
 

[edit] Clojure

(import '(java.util Random))
(def normals
(let [r (Random.)]
(take 1000 (repeatedly #(-> r .nextGaussian (* 0.5) (+ 1.0))))))

[edit] Common Lisp

(loop for i from 1 to 1000
collect (1+ (* (sqrt (* -2 (log (random 1.0)))) (cos (* 2 pi (random 1.0))) 0.5)))

[edit] D

import std.stdio, std.random, std.math;
 
struct NormalRandom {
double mean, stdDev;
 
// Necessary 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 r1 = uniform01, r2 = uniform01; // Not nothrow.
return mean + stdDev * sqrt(-2 * r1.log) * cos(2 * PI * r2);
}
}
 
void main() {
double[1000] array;
auto nRnd = NormalRandom(1.0, 0.5);
foreach (ref x; array)
//x = nRnd;
x = nRnd();
}

[edit] Alternative Version

(Untested)

Library: tango
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;
}
}

[edit] Delphi

Delphi has RandG function which generates random numbers with normal distribution using Marsaglia-Bray algorithm:

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.

Output:

Mean          = 1.0098
Std deviation = 0.5016

[edit] DWScript

var values : array [0..999] of Float;
var i : Integer;
 
for i := values.Low to values.High do
values := RandG(1, 0.5);

[edit] E

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

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

Example Result

Average: 1.0079398405028137
Standard Deviation: 0.49042787564453988

[edit] Elixir

 
defmodule Random do
def init() do
 :random.seed(:erlang.now())
end
def normal(mean, sd) do
{a, b} = {:random.uniform(), :random.uniform()}
mean + sd * (:math.sqrt(-2 * :math.log(a)) * :math.cos(2 * :math.pi * b))
end
end
 
Random.init()
xs = for _ <- 1..1000, do: Random.normal(1.0, 0.5)
 

[edit] Erlang

Works with: 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)]).
 

Output:

mean = 1.0118289913718608
stddev = 0.5021636849524854

[edit] Euler Math Toolbox

 
>v=normal(1,1000)*0.5+1;
>mean(v), dev(v)
1.00291801071
0.498226876528
 

[edit] Euphoria

Translation of: PureBasic
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

[edit] Factor

1000 [ 1.0 0.5 normal-random-float ] replicate

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

[edit] 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:

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

The second calls out to Java's Gaussian random-number generator:

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

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

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

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

[edit] 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()
}
}

[edit] Groovy

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

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

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

[edit] 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.

 
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
 

[edit] IDL

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

[edit] J

Solution:

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

Alternative Solution:
Using the normal script from the stats/distribs addon.

   require 'stats/distribs/normal'
1 0.5 rnorm 1000
1.44868803 1.21548637 0.812460657 1.54295452 1.2470606 ...

[edit] 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();
}

[edit] 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
}

[edit] Julia

Julia's standard library provides a randn function to generate normally distributed random numbers (with mean 0 and standard deviation 0.5, which can be easily rescaled to any desired values):

randn(1000) * 0.5 + 1

[edit] LabVIEW

Works with: LabVIEW version 8.6

LV array of randoms with given mean and stdev.png

[edit] Liberty BASIC

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

[edit]

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


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

[edit] Mathematica

Built-in function RandomReal with built-in distribution NormalDistribution as an argument:

RandomReal[NormalDistribution[1, 1/2], 1000]

[edit] MATLAB

Native support :

    mu = 1; sd = 0.5; 
x = randn(1000,1) * sd + mu;
 

The statistics toolbox provides this function

   x = normrnd(mu, sd, [1000,1]); 

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.

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

Output:

>> randNorm(1,.5, [1000,1])
 
ans =
 
0.693984121077029

[edit] Maxima

load(distrib)$
 
random_normal(1.0, 0.5, 1000);

[edit] 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
)

[edit] Metafont

Metafont has normaldeviate which produces pseudorandom normal distributed numbers with mean 0 and variance one. So the following complete the task:

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

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

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

[edit] МК-61/52

П7	<->	П8	1/x	П6	ИП6	П9	СЧ	П6	1/x
ln ИП8 * 2 * КвКор ИП9 2 * пи
* sin * ИП7 + С/П БП 05

Input: РY - variance, РX - expectation.

Or:

3	10^x	П0	ПП	13	2	/	1	+	С/П	L0	03	С/П
СЧ lg 2 /-/ * КвКор 2 пи ^ СЧ * * cos * В/О

to generate 1000 numbers with a mean of 1.0 and a standard deviation of 0.5.

[edit] Modula-3

Translation of: C
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.

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

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

[edit] NewLISP

(normal 1 .5 1000)

[edit] Nimrod

import math, strutils
 
const precisn = 5
var rs: TRunningStat
 
proc normGauss: float {.inline.} = 1 + 0.76 * cos(2*PI*random(1.0)) * sqrt(-2*log10(random(1.0)))
 
randomize()
 
for j in 0..5:
for i in 0..1000:
rs.push(normGauss())
echo("mean: ", $formatFloat(rs.mean,ffDecimal,precisn),
" stdDev: ", $formatFloat(rs.standardDeviation(),ffDecimal,precisn))
Output:
mean: 1.01703 stdDev: 0.50324
mean: 1.01187 stdDev: 0.50060
mean: 1.00216 stdDev: 0.49969
mean: 1.00335 stdDev: 0.50184
mean: 1.00120 stdDev: 0.49830
mean: 1.00217 stdDev: 0.49911

[edit] 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();
}
}
}

[edit] 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 ());;

[edit] Octave

p = normrnd(1.0, 0.5, 1000, 1);
disp(mean(p));
disp(sqrt(sum((p - mean(p)).^2)/numel(p)));

Output:

1.0209
0.51048


[edit] PARI/GP

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)

[edit] Pascal

See Delphi

[edit] Perl

my $PI = 2 * atan2 1, 0;
 
my @nums = map {
1 + 0.5 * sqrt(-2 * log rand) * cos(2 * $PI * rand)
} 1..1000;

[edit] Perl 6

Works with: Rakudo version #22 "Thousand Oaks"
sub randnorm ($mean, $stddev) {
$mean + $stddev * sqrt(-2 * log rand) * cos(2 * pi * rand)
}
 
my @nums = randnorm(1, 0.5) xx 1000;
 
# Checking
say my $mean = @nums R/ [+] @nums;
say my $stddev = sqrt $mean**2 R- @nums R/ [+] @nums X** 2;
 

[edit] 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);
 
}

[edit] PicoLisp

Translation of: C
(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) " ") ) )

Output:

1.500334 1.212931 1.095283 0.433122 0.459116 1.302446 0.402477

[edit] 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;
 

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   

[edit] 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;

[edit] 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;

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

[edit] Python

Works with: Python version 2.5
import random
values = [random.gauss(1, .5) for i in range(1000)]
# or [ random.normalvariate(1, 0.5) for i in range(1000)]

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

Quick check

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

[edit] R

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

[edit] Racket

 
#lang racket
(for/list ([i 1000])
(add1 (* (sqrt (* -2 (log (random)))) (cos (* 2 pi (random))) 0.5)))
 

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

Quick Check (on linux with code in file rand.rv)

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

[edit] 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.

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

output when using the default input

              old mean= 0.50398978
old standard deviation= 0.292088074

              new mean= 0.999477923
new standard deviation= 0.505691743

[edit] Ruby

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

[edit] Run BASIC

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

[edit] Rust

use std::rand;
use std::rand::distributions::{Normal, IndependentSample};
 
fn main() {
let mut rands = [0.0, ..1000];
let normal = Normal::new(1.0, 0.5);
 
for i in range(0, 1000) {
let v = normal.ind_sample(&mut rand::task_rng());
rands[i] = v;
}
}

[edit] 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;
 

Results:

 The MEANS Procedure

                     Analysis Variable : r

                          Mean         Std Dev
                  ----------------------------
                     0.9907408       0.4844051
                  ----------------------------

[edit] 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;

[edit] Scala

List.fill(1000)(1.0 + 0.5 * scala.util.Random.nextGaussian)

[edit] 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
#|
(-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)

[edit] 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;

[edit] 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.

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

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.

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

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.

[edit] 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()}]
}

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

[edit] TorqueScript

for (%i = 0; %i < 1000; %i++)
%list[%i] = 1 + mSqrt(-2 * mLog(getRandom())) * mCos(2 * $pi * getRandom());

[edit] 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.

#import nat
#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
 
#cast %eWL
 
test =
 
^(mean,stdev)* <
pop_stats(1.,0.5) 1000,
sample_stats(1.,0.5) 1000>

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

[edit] Yorick

Returns array of count random numbers with mean 0 and standard deviation 1.

func random_normal(count) {
return sqrt(-2*log(random(count))) * cos(2*pi*random(count));
}

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

[edit] zkl

fcn mkRand(mean,sd){ //normally distributed random w/mean & standard deviation
pi:=(0.0).pi; // using the Box–Muller transform
rz1:=fcn{1.0-(0.0).random(1)} // from [0,1) to (0,1]
return('wrap(){((-2.0*rz1().log()).sqrt() * (2.0*pi*rz1()).cos())*sd + mean })
}

This creates a new random number generator, now to use it:

var g=mkRand(1,0.5);
ns:=(0).pump(1000,List,g); // 1000 rands with mean==1 & sd==1/2
mean:=(ns.sum(0.0)/1000); //-->1.00379
// calc sd of list of numbers:
(ns.reduce('wrap(p,n){p+(n-mean).pow(2)},0.0)/1000).sqrt() //-->0.494844

[edit] ZX Spectrum Basic

Here we have converted the QBasic code to suit the ZX Spectrum:

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
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox