Random numbers: Difference between revisions
m Make together {{omit from}}s. |
|||
Line 739: | Line 739: | ||
}; |
}; |
||
vector(1000,unused,rnormal()/2+1)</lang> |
vector(1000,unused,rnormal()/2+1)</lang> |
||
=={{header|Pascal}}== |
|||
See [[Random_numbers#Delphi | Delphi]] |
|||
=={{header|Perl}}== |
=={{header|Perl}}== |
Revision as of 23:28, 30 October 2011
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.
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
<lang algol68>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]))
)</lang> 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
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
BBC BASIC
<lang bbcbasic> 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</lang>
Output:
Mean = 1.01848064 Standard deviation = 0.503551814
C
<lang 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;
}</lang>
C#
This is the Javascript translation into C#. <lang csharp> private static double randomNormal() { return Math.Cos(2 * Math.PI * tRand.NextDouble()) * Math.Sqrt(-2 * Math.Log(tRand.NextDouble())); } </lang>
Then the methods in Metafont are used to calculate the average and the Standard Deviation: <lang csharp> 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(); } </lang>
An example result:
Average: 1,00510073053613 Standard Deviation: 0,502540443430955
C++
<lang cpp>#include <cstdlib> // for rand
- include <cmath> // for atan, sqrt, log, cos
- include <algorithm> // for generate_n
double const pi = 4*std::atan(1.0);
// simple functor for normal distribution class normal_distribution { public:
normal_distribution(double m, double s): mu(m), sigma(s) {} double operator() // returns a single normally distributed number { double r1 = (std::rand() + 1.0)/(RAND_MAX + 1.0); // gives equal distribution in (0, 1] double r2 = (std::rand() + 1.0)/(RAND_MAX + 1.0); return mu + sigma * std::sqrt(-2*std::log(r1))*std::cos(2*pi*r2); }
private:
double mu, sigma;
};
int main() {
double array[1000]; std::generate_n(array, 1000, normal_distribution(1.0, 0.5)); return 0;
}</lang>
This example used Mersenne Twister generator. It can be changed by changing the typedef.
<lang cpp>
- 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;
} </lang>
C++0x
The new C++ standard looks very similar to the Boost library here.
<lang cpp>#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);
}</lang>
Clojure
<lang lisp>(import '(java.util Random)) (def normals
(let [r (Random.)] (take 1000 (repeatedly #(-> r .nextGaussian (* 0.5) (+ 1.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>
D
<lang d>import std.stdio, std.random, std.math;
struct NormalRandom {
double mean, stdDev;
this(double mean_, double stdDev_) { this.mean = mean_; this.stdDev = stdDev_; }
double opCall() { double r1 = uniform(0.0, 1.0); 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();
}</lang>
<lang d>import tango.math.random.Random; double [1000]list; auto r = new Random(); foreach(ref l; list) {
r.normalSource!(double)()(l); l = 1.0 + 0.5 * l;
}</lang>
Delphi
Delphi has RandG function which generates random numbers with normal distribution using Marsaglia-Bray algorithm:
<lang Delphi>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.</lang> Output:
Mean = 1.0098 Std deviation = 0.5016
E
<lang e>accum [] for _ in 1..1000 { _.with(entropy.nextGaussian()) }</lang>
Euphoria
<lang 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</lang>
Factor
<lang factor>1000 [ 1.0 0.5 normal-random-float ] replicate</lang>
Falcon
<lang 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</lang>
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:
<lang fantom> 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) } }
} </lang>
The second calls out to Java's Gaussian random-number generator:
<lang fantom> 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) } }
} </lang>
Forth
<lang 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</lang>
Fortran
<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#
<lang fsharp>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() ]</lang>
Go
<lang go>package main import "rand" import "time"
func main() {
var list [1000]float64 rand.Seed(time.Nanoseconds()) for i := range list { list[i] = 1.0 + 0.5 * rand.NormFloat64() }
}</lang>
Groovy
<lang groovy>rnd = new Random() result = (1..1000).inject([]) { r, i -> r << rnd.nextGaussian() }</lang>
Haskell
<lang 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</lang>
HicEst
<lang 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 </lang>
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. <lang icon> 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 </lang>
IDL
<lang idl>result = 1.0 + 0.5*randomn(seed,1000)</lang>
J
Solution: <lang j>urand=: ?@$ 0: zrand=: (2 o. 2p1 * urand) * [: %: _2 * [: ^. urand
1 + 0.5 * zrand 100</lang>
Alternative Solution:
Using the normal script from the stats/distribs addon.
<lang j> require 'stats/distribs/normal'
1 0.5 rnorm 1000
1.44868803 1.21548637 0.812460657 1.54295452 1.2470606 ...</lang>
Java
<lang 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();
}</lang>
JavaScript
<lang 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
}</lang>
LabVIEW
Liberty BASIC
<lang lb>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</lang>
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. <lang logo>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 ?] []</lang>
lua
<lang 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 </lang>
Mathematica
Built-in function RandomReal with built-in distribution NormalDistribution as an argument: <lang Mathematica>RandomReal[NormalDistribution[1, 1/2], 1000]</lang>
MATLAB
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.
<lang MATLAB>function randNum = randNorm(mu0,chi2)
radiusSquared = +Inf;
while (radiusSquared >= 1) u = ( 2 * rand(1) ) - 1; v = ( 2 * rand(1) ) - 1;
radiusSquared = u^2 + v^2; end
scaleFactor = sqrt( ( -2*log(radiusSquared) )/ radiusSquared ); randNum = (v * scaleFactor * chi2) + mu0;
end</lang>
Output: <lang MATLAB>>> randNorm(1,.5)
ans =
0.693984121077029</lang>
MAXScript
<lang 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
)</lang>
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 what 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
Mirah
<lang 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 </lang>
Modula-3
<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>
Objeck
<lang 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(); } }
} </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
PARI/GP
<lang parigp>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)</lang>
Pascal
See Delphi
Perl
<lang perl>my $PI = 2 * atan2 1, 0;
my @nums = map {
1 + 0.5 * sqrt(-2 * log rand) * cos(2 * $PI * rand)
} 1..1000;</lang>
Perl 6
<lang perl6>sub randnorm ($mean, $stddev) {
$mean + $stddev * sqrt(-2 * log rand) * cos(2 * pi * rand)
}
my @nums = map { randnorm 1, 0.5 }, ^1000;</lang>
PHP
<lang 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);
}</lang>
PicoLisp
<lang 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) " ") ) )</lang>
Output:
1.500334 1.212931 1.095283 0.433122 0.459116 1.302446 0.402477
PL/I
<lang 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; </lang>
Pop11
<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;</lang>
PureBasic
<lang 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.f(n-1) For i = 0 To n-1
a(i) = 1 + 0.5 * RandomNormal()
Next</lang>
Python
<lang python>import random values = [random.gauss(1, .5) for i in range(1000)]
- 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.
Quick check <lang python>>>> mean = sum(values)/1000 >>> sdeviation = (sum((i - mean)**2 for i in values)/1000)**0.5 >>> mean, sdeviation (1.0127861555468178, 0.5006682783828207) </lang>
R
<lang r>result <- rnorm(1000, mean=1, sd=0.5)</lang>
Ruby
<lang ruby>Array.new(1000) { 1 + Math.sqrt(-2 * Math.log(rand)) * Math.cos(2 * Math::PI * rand) }</lang>
Sather
<lang 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;</lang>
Scala
<lang scala>List.fill(1000)(1.0 + 0.5 * scala.util.Random.nextGaussian)</lang>
Seed7
<lang 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;</lang>
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.
<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 int
s). 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
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. <lang Ursala>#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></lang>
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)>
Yorick
Returns array of count random numbers with mean 0 and standard deviation 1. <lang yorick>func random_normal(count) {
return sqrt(-2*log(random(count))) * cos(2*pi*random(count));
}</lang>
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
ZX Spectrum Basic
Here we have converted the QBasic code to suit the ZX Spectrum:
<lang zxbasic>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</lang>
- 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
- E
- 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
- MAXScript
- Metafont
- Mirah
- Modula-3
- Objeck
- OCaml
- Octave
- PARI/GP
- Pascal
- Perl
- Perl 6
- PHP
- PicoLisp
- PL/I
- Pop11
- PureBasic
- Python
- R
- Ruby
- Sather
- Scala
- Seed7
- Standard ML
- Tcl
- TI-83 BASIC
- Ursala
- Yorick
- ZX Spectrum Basic
- GAP/Omit
- GUISS/Omit
- UNIX Shell/Omit
- Randomness