Random numbers
From Rosetta Code
You are encouraged to solve this task according to the task description, using any language you may know.
Many libraries only generate uniformly distributed random numbers. If so, use this formula to convert them to a normal distribution.
[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
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
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#
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++
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<> dist(1.0, 0.5);
auto rnd = bind(dist, engine);
vector<double> v(1000);
generate(v.begin(), v.end(), rnd);
return 0;
}
#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;
}
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;
// 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();
}
[edit] Alternative Version
(Untested)
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] 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
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
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
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] LabVIEW
[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] 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.
[edit] Modula-3
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] 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
sub randnorm ($mean, $stddev) {
$mean + $stddev * sqrt(-2 * log rand) * cos(2 * pi * rand)
}
my @nums = map { randnorm 1, 0.5 }, ^1000;
[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
(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
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] 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
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] 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] 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
Categories:
- Programming Tasks
- Basic language learning
- Probability and statistics
- Ada
- ALGOL 68
- AutoHotkey
- AWK
- BASIC
- BBC BASIC
- C
- C sharp
- C++
- Boost
- Clojure
- Common Lisp
- D
- Tango
- Delphi
- DWScript
- E
- Eiffel
- Erlang
- Euler Math Toolbox
- Euphoria
- Factor
- Falcon
- Fantom
- Forth
- Fortran
- F Sharp
- Go
- Groovy
- Haskell
- HicEst
- Icon
- Unicon
- IDL
- J
- Java
- JavaScript
- LabVIEW
- Liberty BASIC
- Logo
- Lua
- Mathematica
- MATLAB
- Maxima
- MAXScript
- Metafont
- Mirah
- МК-61/52
- Modula-3
- NetRexx
- NewLISP
- Objeck
- OCaml
- Octave
- PARI/GP
- Pascal
- Perl
- Perl 6
- PHP
- PicoLisp
- PL/I
- PL/SQL
- Pop11
- PureBasic
- Python
- R
- Racket
- Raven
- REXX
- Ruby
- Run BASIC
- Sather
- Scala
- Scheme
- Seed7
- Standard ML
- Tcl
- TI-83 BASIC
- Ursala
- Yorick
- ZX Spectrum Basic
- GUISS/Omit
- UNIX Shell/Omit
- Randomness
