Random numbers: Difference between revisions

Generate a collection filled with   '''1000'''   normally distributed random (or pseudo-random) numbers

with a mean of   '''1.0'''   and a   [[wp:Standard_deviation|standard deviation]]   of   '''0.5'''

 

Many libraries only generate uniformly distributed random numbers. If so, you may use [[wp:Normal_distribution#Generating_values_from_normal_distribution|one of these algorithms]].

 

;Related task:

 

=={{header|Ada}}==

with Ada.Numerics.Float_Random; use Ada.Numerics.Float_Random;

with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;

Distribution (I) := Normal_Distribution (Seed);

end loop;

 

=={{header|ALGOL 68}}==

{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-1.18.0/algol68g-1.18.0-9h.tiny.el5.centos.fc11.i386.rpm/download 1.18.0-9h.tiny]}}

{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of FORMATted transput}}

(

sqrt(-2*log(random)) * cos(2*pi*random)

INT limit=10;

printf(($"("n(limit-1)(-d.6d",")-d.5d" ... )"$, rands[:limit]))

{{out}}

<pre>

=={{header|Arturo}}==

 

 

rands: map 1..1000 'x [

]

 

 

{{out}}

contributed by Laszlo on the ahk

[http://www.autohotkey.com/forum/post-276261.html#276261 forum]

R .= RandN(1,0.5) "`n" ; mean = 1.0, standard deviation = 0.5

MsgBox %R%

}

Return Y

 

=={{header|Avail}}==

[

lower : number,

// the default distribution has mean 0 and std dev 1.0, so we scale the values

sampler ::= map a Marsaglia polar sampler through [d : double | d ÷ 2.0 + 1.0];

 

=={{header|AWK}}==

'''One-liner:'''

 

'''Readable version:'''

function r() {

return sqrt( -2*log( rand() ) ) * cos(6.2831853*rand() )

print s

}

{{out}} first few values only

<pre>

 

=={{header|BASIC}}==

 

FOR number% = 0 TO 999

array(number%) = 1.0 + 0.5 * SQR(-2*LN(RND(1))) * COS(2*PI*RND(1))

PRINT "Mean = " ; mean

{{out}}

<pre>Mean = 1.01848064

Standard deviation = 0.503551814</pre>

 

=={{header|C}}==

#include <math.h>

#ifndef M_PI

rands[i] = 1.0 + 0.5*random_normal();

return 0;

 

=={{header|C sharp}}==

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

 

Console.ReadLine();

}

 

An example result:

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

 

#include <functional>

#include <vector>

generate(v.begin(), v.end(), rnd);

return 0;

 

{{works with|C++03}}

#include <cmath> // for atan, sqrt, log, cos

#include <algorithm> // for generate_n

std::generate_n(array, 1000, normal_distribution(1.0, 0.5));

return 0;

 

{{libheader|Boost}}

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

 

#include <vector>

#include "boost/random.hpp"

return 0;

}

 

=={{header|Clojure}}==

(def normals

(let [r (Random.)]

 

=={{header|Common Lisp}}==

 

=={{header|Crystal}}==

{{trans|Phix}}

array = Array.new(n) { mean + sd * Math.sqrt(-2 * Math.log(rand)) * Math.cos(tau * rand) }

 

mean = array.sum / array.size

standev = Math.sqrt( array.sum{ |x| (x - mean) ** 2 } / array.size )

{{out}}

<pre>mean = 1.0093442539237896, standard deviation = 0.504694489463623</pre>

 

=={{header|D}}==

 

struct NormalRandom {

//x = nRnd;

x = nRnd();

 

===Alternative Version===

{{libheader|tango}}

 

 

void main() {

l = 1.0 + 0.5 * l;

}

 

=={{header|Delphi}}==

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

 

 

{$APPTYPE CONSOLE}

Writeln('Std Deviation = ', StdDev(Values):6:4);

Readln;

{{out}}

<pre>Mean = 1.0098

 

=={{header|DWScript}}==

var i : Integer;

 

for i := values.Low to values.High do

 

=={{header|E}}==

 

=={{header|EasyLang}}==

 

.

 

=={{header|Eiffel}}==

class

APPLICATION

end

end

 

Example Result

=={{header|Elena}}==

{{trans|C#}}

import extensions'math;

real tAvg := 0;

{

a[x] := (randomNormal()) / 2 + 1;

real s := 0;

{

s += power(a[x] - tAvg, 2)

console.readChar()

{{out}}

<pre>

 

=={{header|Elixir}}==

def normal(mean, sd) do

{a, b} = {:rand.uniform, :rand.uniform}

 

xs = for _ <- 1..1000, do: Random.normal(1.0, 0.5)

{{out}}

<pre>

 

used Erlang function <code>:rand.normal</code>

{{out}}

<pre>

{{works with|Erlang}}

 

mean(Values) ->

mean(tl(Values), hd(Values), 1).

io:format("mean = ~w\n", [mean(X)]),

io:format("stddev = ~w\n", [stddev(X)]).

{{out}}

<pre>

 

=={{header|ERRE}}==

PROGRAM DISTRIBUTION

 

 

END PROGRAM

 

=={{header|Euler Math Toolbox}}==

 

>v=normal(1,1000)*0.5+1;

>mean(v), dev(v)

1.00291801071

0.498226876528

 

=={{header|Euphoria}}==

{{trans|PureBasic}}

 

function RandomNormal()

for i = 1 to n do

s[i] = 1 + 0.5 * RandomNormal()

 

=={{header|F_Sharp|F#}}==

let n = MathNet.Numerics.Distributions.Normal(1.0,0.5)

List.init 1000 (fun _->n.Sample())

{{out}}

<pre>

0.6224640457; 0.8524078517; 0.7646595627; 0.6799834691; 0.773111053; ...]

</pre> =={{header|F_Sharp|F#}}==

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

 

=={{header|Factor}}==

 

for i in [0:1000] : a+= norm_rand_num()

 

pi = 2*acos(0)

return 1 + (cos(2 * pi * random()) * pow(-2 * log(random()) ,1/2)) /2

 

=={{header|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

{

}

}

 

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

 

using [java] java.util::Random

 

}

}

 

=={{header|Forth}}==

{{works with|gforth|0.6.2}}

 

here to seed

 

0 do frnd-normal 0.5e f* 1e f+ f, loop ;

 

 

For newer versions of gforth (tested on 0.7.3), it seems you need to use <tt>HERE SEED !</tt> instead of <tt>HERE TO SEED</tt>, because <tt>SEED</tt> has been made a variable instead of a value.

 

 

=={{header|Fortran}}==

{{works with|Fortran|90 and later}}

 

INTEGER, PARAMETER :: n = 1000

WRITE(*, "(A,F8.6)") "Standard Deviation = ", sd

 

 

{{out}}

Free Pascal provides the '''randg''' function in the RTL math unit that produces Gaussian-distributed random numbers with the Box-Müller algorithm.

 

function randg(mean,stddev: float): float;

 

=={{header|Frink}}==

 

=={{header|Go}}==

This solution uses math/rand package in the standard library. See also though the subrepository rand package at https://godoc.org/golang.org/x/exp/rand, which also has a NormFloat64 and has a rand source with a number of advantages over the one in standard library.

 

import (

fmt.Println(strings.Repeat("*", (c+hscale/2)/hscale))

}

{{out}}

<pre>

 

=={{header|Groovy}}==

 

=={{header|Haskell}}==

 

 

pairs :: [a] -> [(a,a)]

 

result :: IO [Double]

 

Or using Data.Random from random-fu package:

To print them:

import Control.Monad

 

main = do

x <- sample thousandRandomNumbers

 

=={{header|HicEst}}==

 

pi = 4 * ATAN(1)

 

=={{header|Icon}} and {{header|Unicon}}==

 

Note that Unicon randomly seeds it's generator.

procedure main()

local L

every write(!L)

end

 

=={{header|IDL}}==

 

=={{header|J}}==

'''Solution:'''

zrand=: (2 o. 2p1 * urand) * [: %: _2 * [: ^. urand

 

 

'''Alternative Solution:'''<br>

Using the normal script from the [[j:Addons/stats/distribs|stats/distribs addon]].

1 0.5 rnorm 1000

 

=={{header|Java}}==

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

 

=={{header|JavaScript}}==

return Math.cos(2 * Math.PI * Math.random()) * Math.sqrt(-2 * Math.log(Math.random()))

}

for (var i=0; i < 1000; i++){

a[i] = randomNormal() / 2 + 1

 

=={{header|jq}}==

 

''''A Pseudo-Random Number Generator''''

# The random numbers are in [0 -- 32767] inclusive.

# Input: an array of length at least 2 interpreted as [count, state, ...]

.[0] as $count | .[1] as $state

| ( (214013 * $state) + 2531011) % 2147483648 # mod 2^31

''''Box-Muller Method''''

# using the Box-Muller method: X = sqrt(-2 ln U) * cos(2 pi V) where U and V are uniform on [0,1].

# Input: [n, state]

next_rand_normal

| recurse( if .[0] < count then next_rand_normal else empty end)

'''Example'''

The task can be completed using: [0,1] | random_normal_variate(1; 0.5; 1000) | .[2]

 

We show just the sample average and standard deviation:

length as $l | add as $sum | ($sum/$l) as $a

| reduce .[] as $x (0; . + ( ($x - $a) | .*. ))

| [ $a, (./$l | sqrt)] ;

 

{{out}}

$ jq -n -c -f Random_numbers.jq

=={{header|Julia}}==

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

 

=={{header|Kotlin}}==

 

import java.util.Random

println("Mean is $mean")

println("S.D. is $sd")

Sample output:

{{out}}

{{works with|LabVIEW|8.6}}

[[File:LV_array_of_randoms_with_given_mean_and_stdev.png]]

 

=={{header|Lingo}}==

on randf ()

n = random(the maxinteger)-1

return n / float(the maxinteger-1)

 

repeat with i = 1 to 1000

normal.add(1 + sqrt(-2 * log(randf())) * cos(2 * PI * randf()) / 2)

 

=={{header|Lobster}}==

Uses built-in <code>rnd_gaussian</code>

let mean = 1.0

let stdv = 0.5

 

test_random_normal()

 

=={{header|Logo}}==

{{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.

localmake "max.int lshift -1 -1

output quotient random :max.int :max.int

end

 

 

=={{header|Lua}}==

for i = 1, 1000 do

list[i] = 1 + math.sqrt(-2 * math.log(math.random())) * math.cos(2 * math.pi * math.random()) / 2

 

=={{header|M2000 Interpreter}}==

M2000 use a Wichmann - Hill Pseudo Random Number Generator.

Module CheckIt {

Function StdDev (A()) {

}

Checkit

 

=={{header|Maple}}==

 

'''or'''

 

 

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

 

=={{header|MATLAB}}==

 

Native support :

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

 

The statistics toolbox provides this function

 

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.

 

radiusSquared = +Inf;

randNum = (v .* scaleFactor .* chi2) + mu0;

 

 

Output:

 

ans =

 

 

=={{header|Maxima}}==

 

 

=={{header|MAXScript}}==

 

for i in 1 to 1000 do

(

c = 1.0 + 0.5 * sqrt (-2*log a) * cos (360*b) -- Maxscript cos takes degrees

append arr c

 

=={{header|Metafont}}==

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

 

 

m := 0; % m holds the mean, for testing purposes

 

show m, s; % and let's show that really they get what we wanted

 

A run gave

 

=={{header|MiniScript}}==

return mean + sqrt(-2 * log(rnd,2.7182818284)) * cos(2*pi*rnd) * stddev

end function

for i in range(1,1000)

x.push randNormal(1, 0.5)

 

=={{header|Mirah}}==

 

list = double[999]

list[i] = mean + std * rng.nextGaussian

end

 

=={{header|МК-61/52}}==

ln ИП8 * 2 * КвКор ИП9 2 * пи

 

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

Or:

 

 

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

{{trans|C}}

 

 

IMPORT Random;

rands[i] := 1.0D0 + 0.5D0 * RandNorm();

END;

 

=={{header|Nanoquery}}==

{{trans|Java}}

mean = 1.0; std = 0.5

rng = new(Nanoquery.Util.Random)

for i in range(0, len(list) - 1)

list[i] = mean + std * rng.getGaussian()

 

=={{header|NetRexx}}==

options replace format comments java crossref symbols nobinary

 

 

return

{{out}}

<pre>

 

=={{header|NewLISP}}==

 

=={{header|Nim}}==

 

var rs: RunningStat

 

for _ in 1..5:

{{out}}

 

=={{header|Objeck}}==

class RandomNumbers {

function : Main(args : String[]) ~ Nil {

}

}

 

=={{header|OCaml}}==

let random_gaussian () =

1. +. sqrt (-2. *. log (Random.float 1.)) *. cos (2. *. pi *. Random.float 1.);;

 

=={{header|Octave}}==

disp(mean(p));

 

{{out}}

{{trans|REXX}}

===version 1===

pi=RxCalcPi() /* get value of pi */

Parse Arg n seed . /* allow specification of N & seed*/

Return RxCalcPower(_/n,.5)

 

{{out}}

<pre> old mean= 0.49830002

===version 2===

Using the nice function names in the algorithm.

pi=RxCalcPi() /* get value of pi */

Parse Arg n seed . /* allow specification of N & seed*/

sin: Return RxCalcSin(arg(1),,'R')

 

 

=={{header|PARI/GP}}==

my(pr=32*ceil(default(realprecision)*log(10)/log(4294967296)),u1=random(2^pr)*1.>>pr,u2=random(2^pr)*1.>>pr);

\\ Could easily be extended with a second normal at very little cost.

};

 

=={{header|Pascal}}==

 

The following function calculates Gaussian-distributed random numbers with the Box-Müller algorithm:

function rnorm (mean, sd: real): real;

{Calculates Gaussian random numbers according to the Box-Müller approach}

u1 := random;

u2 := random;

end;

 

[[#Delphi | Delphi]] and [[#Free Pascal|Free Pascal]] support implement a '''randg''' function that delivers Gaussian-distributed random numbers.

 

=={{header|Perl}}==

 

my @nums = map {

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

 

=={{header|Phix}}==

{{Trans|Euphoria}}

 

=={{header|Phixmonti}}==

 

def RandomNormal

get mean - 2 power rot + swap

endfor

 

=={{header|PHP}}==

 

return mt_rand() / mt_getrandmax();

}

$a[$i] = 1.0 + ((sqrt(-2 * log(random())) * cos(2 * $pi * random())) * 0.5);

 

=={{header|PicoLisp}}==

{{trans|C}}

 

(de randomNormal () # Normal distribution, centered on 0, std dev 1

(link (+ 1.0 (/ (randomNormal) 2))) ) )

(for N (head 7 Result) # Print first 7 results

{{out}}

<pre>1.500334 1.212931 1.095283 0.433122 0.459116 1.302446 0.402477</pre>

 

=={{header|PL/I}}==

/* CONVERTED FROM WIKI FORTRAN */

Normal_Random: procedure options (main);

put skip edit ( "Standard Deviation = ", sd) (a, F(18,16));

END Normal_Random;

{{out}}

<pre>

 

=={{header|PL/SQL}}==

DECLARE

--The desired collection

END LOOP;

END;

 

=={{header|Pop11}}==

true -> popradians;

 

for i from 1 to 1000 do 1.0+0.5*random_normal() endfor;

;;; collect them into array

 

=={{header|PowerShell}}==

Equation adapted from Liberty BASIC

{

[CmdletBinding()]

 

$Stats | Format-List

{{out}}

<pre>

StandardDeviation : 0.489099623426272

</pre>

 

=={{header|Python}}==

;Using random.gauss:

>>> values = [random.gauss(1, .5) for i in range(1000)]

 

;Quick check of distribution:

count = len(numbers)

mean = sum(numbers) / count

>>> quick_check(values)

(1.0140373306786599, 0.49943411329234066)

 

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

 

;Alternatively using random.normalvariate:

>>> quick_check(values)

(0.990099111944864, 0.5029847005836282)

 

=={{header|R}}==

 

=={{header|Racket}}==

#lang racket

(for/list ([i 1000])

(add1 (* (sqrt (* -2 (log (random)))) (cos (* 2 pi (random))) 0.5)))

 

Alternative:

#lang racket

(require math/distributions)

(sample (normal-dist 1.0 0.5) 1000)

 

=={{header|Raku}}==

{{works with|Rakudo|#22 "Thousand Oaks"}}

 

$mean + $stddev * sqrt(-2 * log rand) * cos(2 * pi * rand)

}

say my $mean = @nums R/ [+] @nums;

say my $stddev = sqrt $mean**2 R- @nums R/ [+] @nums X** 2;

 

=={{header|Raven}}==

-1 acos

2 / 1 +

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

stdev = 0.497773

 

=={{header|REXX}}==

Programming note: &nbsp; note the range of the random numbers: &nbsp; (0,1]

<br>(that is, random numbers from &nbsp; zero──►unity, &nbsp; excluding zero, including unity).

numeric digits 20 /*the default decimal digit precision=9*/

parse arg n seed . /*allow specification of N and the seed*/

m.=9; do j=0 while h>9; m.j=h; h=h%2 + 1; end /*j*/

do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/

'''output''' &nbsp; when using the default inputs:

<pre>

 

=={{header|Ring}}==

for i = 1 to 10

see random(i) + nl

next i

 

 

 

=={{header|Rust}}==

{{libheader|rand}}

'''Using a for-loop:'''

use rand::distributions::{Normal, IndependentSample};

 

*num = normal.ind_sample(&mut rng);

}

 

'''Using iterators:'''

use rand::distributions::{Normal, IndependentSample};

 

(0..1000).map(|_| normal.ind_sample(&mut rng)).collect()

};

 

=={{header|SAS}}==

/* Generate 1000 random numbers with mean 1 and standard deviation 0.5.

SAS version 9.2 was used to create this code.*/

end;

run;

Results:

<pre>

 

=={{header|Sather}}==

main is

a:ARRAY{FLTD} := #(1000);

#OUT + "dev " + dev + "\n";

end;

 

=={{header|Scala}}==

===One liner===

===Academic===

object RandomNumbers extends App {

 

println(calcAvgAndStddev(distribution.take(1000))) // e.g. (1.0061433267806525,0.5291834867560893)

}

 

=={{header|Scheme}}==

(define ((c-rand seed)) (set! seed (remainder (+ (* 1103515245 seed) 12345) 2147483648)) (quotient seed 65536))

 

 

(mean-sdev v)

 

=={{header|Seed7}}==

include "float.s7i";

include "math.s7i";

rands[i] := 1.0 + 0.5 * randomNormal;

end for;

 

=={{header|Sidef}}==

 

=={{header|Standard ML}}==

You can call the generator with <code>()</code> repeatedly to get a word in the range <code>[Rand.randMin, Rand.randMax]</code>.

You can use the <code>Rand.norm</code> function to transform the output into a <code>real</code> from 0 to 1, or use the <code>Rand.range (i,j)</code> function to transform the output into an <code>int</code> of the given range.

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

 

2) Random (a subtract-with-borrow generator). You create the generator by calling <code>Random.rand</code> with a seed (of a pair of <code>int</code>s). You can use the <code>Random.randInt</code> function to generate a random int over its whole range; <code>Random.randNat</code> to generate a non-negative random int; <code>Random.randReal</code> to generate a <code>real</code> between 0 and 1; or <code>Random.randRange (i,j)</code> to generate an <code>int</code> in the given range.

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

 

Other implementations of Standard ML have their own random number generators. For example, Moscow ML has a <code>Random</code> structure that is different from the one from SML/NJ.

The SML Basis Library does not provide a routine for uniform deviate generation, and PolyML does not have one. Using a routine from "Monte Carlo" by Fishman (Springer), in the function uniformdeviate, and avoiding the slow IntInf's:

val urandomlist = fn seed => fn n =>

let

val anyrealseed=1009.0 ;

makeNormals bmconv (urandomlist anyrealseed 2000);

 

=={{header|Stata}}==

set obs 1000

 

=== Mata ===

 

=={{header|Tcl}}==

variable ::pi [expr acos(0)]

proc ::tcl::mathfunc::nrand {} {

for {set i 0} {$i < 1000} {incr i} {

lappend result [expr {$mean + $stddev*nrand()}]

 

=={{header|TorqueScript}}==

 

=={{header|Ursala}}==

a standard normal distribution. Mean and standard deviation

library functions are also used in this example.

#import flo

 

^(mean,stdev)* <

pop_stats(1.,0.5) 1000,

The output shows the mean and standard deviation for both sample vectors,

the latter being exact by construction.

 

=={{header|Visual FoxPro}}==

LOCAL i As Integer, m As Double, n As Integer, sd As Double

py = PI()

RETURN z

ENDFUNC

 

=={{header|Wren}}==

 

var rand = Random.new()

var mean = sum / n

System.print("Actual mean : %(mean)")

 

{{out}}

Actual std dev: 0.4961645117026

</pre>

 

=={{header|Yorick}}==

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

return sqrt(-2*log(random(count))) * cos(2*pi*random(count));

 

Example of basic use:

 

=={{header|zkl}}==

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:

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: