Random numbers: Difference between revisions
(→{{header|Perl 6}}: style: use of the list repetition operator) |
(→{{header|BASIC}}: Added ANSI BASIC.) |
||
(187 intermediate revisions by 76 users not shown) | |||
Line 1: | Line 1: | ||
{{task|Basic language learning}} |
|||
{{task|Basic language learning}}[[Category:Probability and statistics]]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 [[wp:Standard_deviation|standard deviation]] of 0.5 |
|||
[[Category:Probability and statistics]] |
|||
[[Category:Randomness]] |
|||
{{omit from|GUISS}} |
|||
{{omit from|UNIX Shell|From the shell, we simply invoke the awk solution}} |
|||
;Task: |
|||
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: |
|||
Many libraries only generate uniformly distributed random numbers. If so, use [[wp:Normal_distribution#Generating_values_from_normal_distribution|this formula]] to convert them to a normal distribution. |
|||
* [[Standard deviation]] |
|||
<br><br> |
|||
=={{header|Ada}}== |
=={{header|Ada}}== |
||
< |
<syntaxhighlight lang="ada">with Ada.Numerics; use Ada.Numerics; |
||
with Ada.Numerics.Float_Random; use Ada.Numerics.Float_Random; |
with Ada.Numerics.Float_Random; use Ada.Numerics.Float_Random; |
||
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions; |
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions; |
||
Line 26: | Line 38: | ||
Distribution (I) := Normal_Distribution (Seed); |
Distribution (I) := Normal_Distribution (Seed); |
||
end loop; |
end loop; |
||
end Normal_Random;</ |
end Normal_Random;</syntaxhighlight> |
||
=={{header|ALGOL 68}}== |
=={{header|ALGOL 68}}== |
||
{{trans|C}} |
{{trans|C}} |
||
Line 34: | Line 47: | ||
{{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]}} |
{{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}} |
{{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}} |
||
< |
<syntaxhighlight lang="algol68">PROC random normal = REAL: # normal distribution, centered on 0, std dev 1 # |
||
( |
( |
||
sqrt(-2*log(random)) * cos(2*pi*random) |
sqrt(-2*log(random)) * cos(2*pi*random) |
||
Line 46: | Line 59: | ||
INT limit=10; |
INT limit=10; |
||
printf(($"("n(limit-1)(-d.6d",")-d.5d" ... )"$, rands[:limit])) |
printf(($"("n(limit-1)(-d.6d",")-d.5d" ... )"$, rands[:limit])) |
||
)</ |
)</syntaxhighlight> |
||
{{out}} |
|||
Output sample: |
|||
<pre> |
<pre> |
||
( 0.693461, 0.948424, 0.482261, 1.045939, 0.890818, 1.467935, 0.604153, 0.804811, 0.690227, 0.83462 ... ) |
( 0.693461, 0.948424, 0.482261, 1.045939, 0.890818, 1.467935, 0.604153, 0.804811, 0.690227, 0.83462 ... ) |
||
</pre> |
</pre> |
||
=={{header|Arturo}}== |
|||
<syntaxhighlight lang="rebol">rnd: function []-> (random 0 10000)//10000 |
|||
rands: map 1..1000 'x [ |
|||
1 + (sqrt neg 2 * ln rnd) * (cos 2 * pi * rnd) |
|||
] |
|||
print rands</syntaxhighlight> |
|||
{{out}} |
|||
<pre>0.6219537961087694 1.279396486161406 1.619019280815647 2.119538294228789 -0.1598383851981044 2.67797211803156 0.9304703037226587 1.629254364659528 -0.4171704717398712 0.9082342931486092 -0.5929704390625219 2.117000897984871 -0.1981633787460266 0.01132471973856153 2.102359263212924 0.2408823232884222 2.046195035792376 0.6374831627030295 0.000839808324124558 1.117061838266626 0.7413355299469649 0.4485598815755762 2.999434800016997 2.560580541932842 1.703197984879731 2.889159248353575 -1.800157205708138 1.756810020187321 0.7136708180852145 0.5929151678321705 0.332519993787973 2.660212054362758 0.5835660585480075 0.8527946892567934 1.640573993747053 0.09471843345263908 1.051402997891346 1.116149156137905 -0.7400139019343499 1.782572831979232 2.531779039786426 0.5240268064639871 0.07099232630526586 -0.854892656700071 1.54381929430469 -0.4416899008614745 0.4274356035015117 0.7350027625573482 2.153583935076981 1.461215281535983 -1.041723064151266 2.338060763553139 -0.1492967916030414 0.3799517724040202 0.4577924541353815 0.673317567666373 -2.27731583876462 1.28480355806061 -0.6925811023772748 -0.2642224122781984 0.6590513830891744 2.55537133425143 1.67933335469247 0.8659013395355968 -1.211026941441126 0.9524579534222226 -0.1931750631835656 -0.5119479869237693 0.1814749003063878 3.03139579963414 ...</pre> |
|||
=={{header|AutoHotkey}}== |
=={{header|AutoHotkey}}== |
||
contributed by Laszlo on the ahk |
contributed by Laszlo on the ahk |
||
[http://www.autohotkey.com/forum/post-276261.html#276261 forum] |
[http://www.autohotkey.com/forum/post-276261.html#276261 forum] |
||
< |
<syntaxhighlight lang="autohotkey">Loop 40 |
||
R .= RandN(1,0.5) "`n" ; mean = 1.0, standard deviation = 0.5 |
R .= RandN(1,0.5) "`n" ; mean = 1.0, standard deviation = 0.5 |
||
MsgBox %R% |
MsgBox %R% |
||
Line 69: | Line 96: | ||
} |
} |
||
Return Y |
Return Y |
||
}</syntaxhighlight> |
|||
}</lang> |
|||
=={{header|Avail}}== |
|||
<syntaxhighlight lang="avail">Method "U(_,_)" is |
|||
[ |
|||
lower : number, |
|||
upper : number |
|||
| |
|||
divisor ::= ((1<<32)) ÷ (upper - lower)→double; |
|||
map a pRNG through [i : integer | (i ÷ divisor) + lower] |
|||
]; |
|||
Method "a Marsaglia polar sampler" is |
|||
[ |
|||
generator for |
|||
[ |
|||
yield : [double]→⊤ |
|||
| |
|||
source ::= U(-1, 1); |
|||
Repeat [ |
|||
x ::= take 1 from source[1]; |
|||
y ::= take 1 from source[1]; |
|||
s ::= x^2 + y^2; |
|||
If 0 < s < 1 then |
|||
[ |
|||
factor ::= ((-2 × ln s) ÷ s) ^ 0.5; |
|||
yield(x × factor); |
|||
yield(y × factor); |
|||
]; |
|||
] |
|||
] |
|||
]; |
|||
// 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]; |
|||
values ::= take 1000 from sampler;</syntaxhighlight> |
|||
=={{header|AWK}}== |
=={{header|AWK}}== |
||
'''One-liner:''' |
|||
<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> |
|||
<syntaxhighlight 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}'</syntaxhighlight> |
|||
'''Readable version:''' |
|||
<syntaxhighlight lang="awk"> |
|||
function r() { |
|||
return sqrt( -2*log( rand() ) ) * cos(6.2831853*rand() ) |
|||
} |
|||
BEGIN { |
|||
n=1000 |
|||
for(i=0;i<n;i++) { |
|||
x = 1 + 0.5*r() |
|||
s = s" "x |
|||
} |
|||
print s |
|||
} |
|||
</syntaxhighlight> |
|||
{{out}} first few values only |
|||
<pre> |
|||
0.783753 1.16682 1.17989 1.14975 1.34784 0.29296 0.979227 1.04402 0.567835 1.58812 0.465559 1.27186 0.324533 0.725827 -0.0626549 0.632273 1.0145 1.3387 0.861667 1.04147 1.2576 1.02497 0.58453 0.9619 1.26902 0.851048 -0.126259 0.863256 </pre> |
|||
... |
|||
=={{header|BASIC}}== |
=={{header|BASIC}}== |
||
==={{header|ANSI BASIC}}=== |
|||
{{works with|QuickBasic|4.5}} |
|||
{{trans|FreeBASIC}} |
|||
RANDOMIZE TIMER 'seeds random number generator with the system time |
|||
{{works with|Decimal BASIC}} |
|||
pi = 3.141592653589793# |
|||
<syntaxhighlight lang="basic"> |
|||
DIM a(1 TO 1000) AS DOUBLE |
|||
100 REM Random numbers |
|||
CLS |
|||
110 RANDOMIZE |
|||
FOR i = 1 TO 1000 |
|||
120 DEF RandomNormal = COS(2 * PI * RND) * SQR(-2 * LOG(RND)) |
|||
130 DIM R(0 TO 999) |
|||
NEXT i |
|||
140 LET Sum = 0 |
|||
150 FOR I = 0 TO 999 |
|||
160 LET R(I) = 1 + RandomNormal / 2 |
|||
170 LET Sum = Sum + R(I) |
|||
180 NEXT I |
|||
190 LET Mean = Sum / 1000 |
|||
200 LET Sum = 0 |
|||
210 FOR I = 0 TO 999 |
|||
220 LET Sum = Sum + (R(I) - Mean) ^ 2 |
|||
230 NEXT I |
|||
240 LET SD = SQR(Sum / 1000) |
|||
250 PRINT "Mean is "; Mean |
|||
260 PRINT "Standard Deviation is"; SD |
|||
270 PRINT |
|||
280 END |
|||
</syntaxhighlight> |
|||
{{out}} Two runs. |
|||
<pre> |
|||
Mean is 1.00216454061435 |
|||
Standard Deviation is .504515904812839 |
|||
</pre> |
|||
<pre> |
|||
Mean is .995781408878628 |
|||
Standard Deviation is .499307289407576 |
|||
</pre> |
|||
==={{header|Applesoft BASIC}}=== |
|||
The [[Random_numbers#Commodore_BASIC|Commodore BASIC]] code works in Applesoft BASIC. |
|||
==={{header|BASIC256}}=== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="basic256"># Generates normally distributed random numbers with mean 0 and standard deviation 1 |
|||
function randomNormal() |
|||
return cos(2.0 * pi * rand) * sqr(-2.0 * log(rand)) |
|||
end function |
|||
dim r(1000) |
|||
sum = 0.0 |
|||
# Generate 1000 normally distributed random numbers |
|||
# with mean 1 and standard deviation 0.5 |
|||
# and calculate their sum |
|||
for i = 0 to 999 |
|||
r[i] = 1.0 + randomNormal() / 2.0 |
|||
sum += r[i] |
|||
next i |
|||
mean = sum / 1000.0 |
|||
sum = 0.0 |
|||
# Now calculate their standard deviation |
|||
for i = 0 to 999 |
|||
sum += (r[i] - mean) ^ 2.0 |
|||
next i |
|||
sd = sqr(sum/1000.0) |
|||
print "Mean is "; mean |
|||
print "Standard Deviation is "; sd |
|||
end</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Mean is 1.002092 |
|||
Standard Deviation is 0.4838570687</pre> |
|||
=={{header|BBC BASIC}}== |
==={{header|BBC BASIC}}=== |
||
< |
<syntaxhighlight lang="bbcbasic"> DIM array(999) |
||
FOR number% = 0 TO 999 |
FOR number% = 0 TO 999 |
||
array(number%) = 1.0 + 0.5 * SQR(-2*LN(RND(1))) * COS(2*PI*RND(1)) |
array(number%) = 1.0 + 0.5 * SQR(-2*LN(RND(1))) * COS(2*PI*RND(1)) |
||
Line 95: | Line 238: | ||
PRINT "Mean = " ; mean |
PRINT "Mean = " ; mean |
||
PRINT "Standard deviation = " ; stdev</ |
PRINT "Standard deviation = " ; stdev</syntaxhighlight> |
||
{{out}} |
|||
Output: |
|||
<pre>Mean = 1.01848064 |
<pre>Mean = 1.01848064 |
||
Standard deviation = 0.503551814</pre> |
Standard deviation = 0.503551814</pre> |
||
==={{header|Chipmunk Basic}}=== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="basic"> |
|||
10 ' Random numbers |
|||
20 randomize timer |
|||
30 dim r(999) |
|||
40 sum = 0 |
|||
50 for i = 0 to 999 |
|||
60 r(i) = 1+randomnormal()/2 |
|||
70 sum = sum+r(i) |
|||
80 next |
|||
90 mean = sum/1000 |
|||
100 sum = 0 |
|||
110 for i = 0 to 999 |
|||
120 sum = sum+(r(i)-mean)^2 |
|||
130 next |
|||
140 sd = sqr(sum/1000) |
|||
150 print "Mean is ";mean |
|||
160 print "Standard Deviation is ";sd |
|||
170 print |
|||
180 end |
|||
500 sub randomnormal() |
|||
510 randomnormal = cos(2*pi*rnd(1))*sqr(-2*log(rnd(1))) |
|||
520 end sub |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
Two runs. |
|||
<pre> |
|||
Mean is 1.007087 |
|||
Standard Deviation is 0.496848 |
|||
</pre> |
|||
<pre> |
|||
Mean is 0.9781 |
|||
Standard Deviation is 0.508147 |
|||
</pre> |
|||
==={{header|Commodore BASIC}}=== |
|||
<syntaxhighlight lang="bbcbasic"> |
|||
10 DIM AR(999): DIM DE(999) |
|||
20 FOR N = 0 TO 999 |
|||
30 AR(N)= 0 + SQR(-1.3*LOG(RND(1))) * COS(1.2*PI*RND(1)) |
|||
40 NEXT N |
|||
50 : |
|||
60 REM SUM |
|||
70 LET SU = 0 |
|||
80 FOR N = 0 TO 999 |
|||
90 LET SU = SU + AR(N) |
|||
100 NEXT N |
|||
110 : |
|||
120 REM MEAN |
|||
130 LET ME= 0 |
|||
140 LET ME = SU/1000 |
|||
150 : |
|||
160 REM DEVIATION |
|||
170 FOR N = 0 TO 999 |
|||
180 T = AR(N)-ME: REM SUBTRACT MEAN FROM NUMBER |
|||
190 T = T * T: REM SQUARE THE RESULT |
|||
200 DE(N) = T : REM STORE IN DEVIATION ARRAY |
|||
210 NEXT N |
|||
220 LET DS=0: REM SUM OF DEVIATION ARRAY |
|||
230 FOR N = 0 TO 999 |
|||
240 LET DS = DS + DE(N) |
|||
250 NEXT N |
|||
260 LET DM=0: REM MEAN OF DEVIATION ARRAY |
|||
270 LET DM = DS / 1000 |
|||
280 LET DE = 0: |
|||
290 LET DE = SQR(DM) |
|||
300 : |
|||
310 PRINT "MEAN = "ME |
|||
320 PRINT "STANDARD DEVIATION ="DE |
|||
330 END |
|||
</syntaxhighlight> |
|||
==={{header|FreeBASIC}}=== |
|||
<syntaxhighlight lang="freebasic">' FB 1.05.0 Win64 |
|||
Const pi As Double = 3.141592653589793 |
|||
Randomize |
|||
' Generates normally distributed random numbers with mean 0 and standard deviation 1 |
|||
Function randomNormal() As Double |
|||
Return Cos(2.0 * pi * Rnd) * Sqr(-2.0 * Log(Rnd)) |
|||
End Function |
|||
Dim r(0 To 999) As Double |
|||
Dim sum As Double = 0.0 |
|||
' Generate 1000 normally distributed random numbers |
|||
' with mean 1 and standard deviation 0.5 |
|||
' and calculate their sum |
|||
For i As Integer = 0 To 999 |
|||
r(i) = 1.0 + randomNormal/2.0 |
|||
sum += r(i) |
|||
Next |
|||
Dim mean As Double = sum / 1000.0 |
|||
Dim sd As Double |
|||
sum = 0.0 |
|||
' Now calculate their standard deviation |
|||
For i As Integer = 0 To 999 |
|||
sum += (r(i) - mean) ^ 2.0 |
|||
Next |
|||
sd = Sqr(sum/1000.0) |
|||
Print "Mean is "; mean |
|||
Print "Standard Deviation is"; sd |
|||
Print |
|||
Print "Press any key to quit" |
|||
Sleep</syntaxhighlight> |
|||
Sample result: |
|||
{{out}} |
|||
<pre> |
|||
Mean is 1.000763573902885 |
|||
Standard Deviation is 0.500653063426955 |
|||
</pre> |
|||
==={{header|FutureBasic}}=== |
|||
Note: To generate the random number, rather than using FB's native "rnd" function, this code wraps C code into the RandomZeroToOne function. |
|||
<syntaxhighlight lang="futurebasic">window 1 |
|||
local fn RandomZeroToOne as double |
|||
double result |
|||
cln result = (double)( (rand() % 100000 ) * 0.00001 ); |
|||
end fn = result |
|||
local fn RandomGaussian as double |
|||
double r = fn RandomZeroToOne |
|||
end fn = 1 + .5 * ( sqr( -2 * log(r) ) * cos( 2 * pi * r ) ) |
|||
long i |
|||
double mean, std, a(1000) |
|||
for i = 1 to 1000 |
|||
a(i) = fn RandomGaussian |
|||
mean += a(i) |
|||
next |
|||
mean = mean / 1000 |
|||
for i = 1 to 1000 |
|||
std += ( a(i) - mean )^2 |
|||
next |
|||
std = std / 1000 |
|||
print " Average: "; mean |
|||
print "Standard Deviation: "; std |
|||
HandleEvents</syntaxhighlight> |
|||
{{output}} |
|||
<pre> |
|||
Average: 1.053724951604593 |
|||
Standard Deviation: 0.2897370762627166 |
|||
</pre> |
|||
==={{header|GW-BASIC}}=== |
|||
The [[Random_numbers#Commodore_BASIC|Commodore BASIC]] code works in GW-BASIC. |
|||
==={{header|Liberty BASIC}}=== |
|||
<syntaxhighlight 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</syntaxhighlight> |
|||
==={{header|Minimal BASIC}}=== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="basic"> |
|||
10 REM Random numbers |
|||
20 LET P = 4*ATN(1) |
|||
30 RANDOMIZE |
|||
40 DEF FNN = COS(2*P*RND)*SQR(-2*LOG(RND)) |
|||
50 DIM R(999) |
|||
60 LET S = 0 |
|||
70 FOR I = 0 TO 999 |
|||
80 LET R(I) = 1+FNN/2 |
|||
90 LET S = S+R(I) |
|||
100 NEXT I |
|||
110 LET M = S/1000 |
|||
120 LET S = 0 |
|||
130 FOR I = 0 TO 999 |
|||
140 LET S = S+(R(I)-M)^2 |
|||
150 NEXT I |
|||
160 LET D = SQR(S/1000) |
|||
170 PRINT "Mean is "; M |
|||
180 PRINT "Standard Deviation is"; D |
|||
190 PRINT |
|||
200 END |
|||
</syntaxhighlight> |
|||
==={{header|PureBasic}}=== |
|||
<syntaxhighlight 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.q(n-1) |
|||
For i = 0 To n-1 |
|||
a(i) = 1 + 0.5 * RandomNormal() |
|||
Next</syntaxhighlight> |
|||
==={{header|QuickBASIC}}=== |
|||
{{works with|QuickBasic|4.5}} |
|||
<syntaxhighlight lang="qbasic"> |
|||
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 |
|||
</syntaxhighlight> |
|||
==={{header|Run BASIC}}=== |
|||
<syntaxhighlight lang="runbasic">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</syntaxhighlight> |
|||
==={{header|TI-83 BASIC}}=== |
|||
Built-in function: randNorm() |
|||
randNorm(1,.5) |
|||
Or by a program: |
|||
Calculator symbol translations: |
|||
"STO" arrow: → |
|||
Square root sign: √ |
|||
ClrList L<sub>1</sub> |
|||
Radian |
|||
For(A,1,1000) |
|||
√(-2*ln(rand))*cos(2*π*A)→L<sub>1</sub>(A) |
|||
End |
|||
==={{header|ZX Spectrum Basic}}=== |
|||
Here we have converted the QBasic code to suit the ZX Spectrum: |
|||
<syntaxhighlight 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</syntaxhighlight> |
|||
=={{header|C}}== |
=={{header|C}}== |
||
< |
<syntaxhighlight lang="c">#include <stdlib.h> |
||
#include <math.h> |
#include <math.h> |
||
#ifndef M_PI |
#ifndef M_PI |
||
Line 122: | Line 521: | ||
rands[i] = 1.0 + 0.5*random_normal(); |
rands[i] = 1.0 + 0.5*random_normal(); |
||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
=={{header|C sharp}}== |
=={{header|C sharp}}== |
||
{{trans|JavaScript}} |
{{trans|JavaScript}} |
||
< |
<syntaxhighlight lang="csharp"> |
||
private static double randomNormal() |
private static double randomNormal() |
||
{ |
{ |
||
return Math.Cos(2 * Math.PI * tRand.NextDouble()) * Math.Sqrt(-2 * Math.Log(tRand.NextDouble())); |
return Math.Cos(2 * Math.PI * tRand.NextDouble()) * Math.Sqrt(-2 * Math.Log(tRand.NextDouble())); |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
Then the methods in [[Random numbers#Metafont]] are used to calculate the average and the Standard Deviation: |
Then the methods in [[Random numbers#Metafont]] are used to calculate the average and the Standard Deviation: |
||
< |
<syntaxhighlight lang="csharp"> |
||
static Random tRand = new Random(); |
static Random tRand = new Random(); |
||
Line 162: | Line 561: | ||
Console.ReadLine(); |
Console.ReadLine(); |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
An example result: |
An example result: |
||
Line 175: | Line 574: | ||
The new C++ standard looks very similar to the Boost library example below. |
The new C++ standard looks very similar to the Boost library example below. |
||
< |
<syntaxhighlight lang="cpp">#include <random> |
||
#include <functional> |
#include <functional> |
||
#include <vector> |
#include <vector> |
||
Line 185: | Line 584: | ||
random_device seed; |
random_device seed; |
||
mt19937 engine(seed()); |
mt19937 engine(seed()); |
||
normal_distribution<> dist(1.0, 0.5); |
normal_distribution<double> dist(1.0, 0.5); |
||
auto rnd = bind(dist, engine); |
auto rnd = bind(dist, engine); |
||
Line 191: | Line 590: | ||
generate(v.begin(), v.end(), rnd); |
generate(v.begin(), v.end(), rnd); |
||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{works with|C++03}} |
{{works with|C++03}} |
||
< |
<syntaxhighlight lang="cpp">#include <cstdlib> // for rand |
||
#include <cmath> // for atan, sqrt, log, cos |
#include <cmath> // for atan, sqrt, log, cos |
||
#include <algorithm> // for generate_n |
#include <algorithm> // for generate_n |
||
Line 220: | Line 619: | ||
std::generate_n(array, 1000, normal_distribution(1.0, 0.5)); |
std::generate_n(array, 1000, normal_distribution(1.0, 0.5)); |
||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{libheader|Boost}} |
{{libheader|Boost}} |
||
Line 226: | Line 625: | ||
This example used Mersenne Twister generator. It can be changed by changing the typedef. |
This example used Mersenne Twister generator. It can be changed by changing the typedef. |
||
< |
<syntaxhighlight lang="cpp"> |
||
#include <vector> |
#include <vector> |
||
#include "boost/random.hpp" |
#include "boost/random.hpp" |
||
Line 247: | Line 646: | ||
return 0; |
return 0; |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Clojure}}== |
=={{header|Clojure}}== |
||
< |
<syntaxhighlight lang="lisp">(import '(java.util Random)) |
||
(def normals |
(def normals |
||
(let [r (Random.)] |
(let [r (Random.)] |
||
(take 1000 (repeatedly #(-> r .nextGaussian (* 0.5) (+ 1.0))))))</ |
(take 1000 (repeatedly #(-> r .nextGaussian (* 0.5) (+ 1.0))))))</syntaxhighlight> |
||
=={{header|COBOL}}== |
|||
<syntaxhighlight lang="cobol"> |
|||
IDENTIFICATION DIVISION. |
|||
PROGRAM-ID. RANDOM. |
|||
AUTHOR. Bill Gunshannon |
|||
INSTALLATION. Home. |
|||
DATE-WRITTEN. 14 January 2022. |
|||
************************************************************ |
|||
** Program Abstract: |
|||
** Able to get the Mean to be really close to 1.0 but |
|||
** couldn't get the Standard Deviation any closer than |
|||
** .3 to .4. |
|||
************************************************************ |
|||
DATA DIVISION. |
|||
WORKING-STORAGE SECTION. |
|||
01 Sample-Size PIC 9(5) VALUE 1000. |
|||
01 Total PIC 9(10)V9(5) VALUE 0.0. |
|||
01 Arith-Mean PIC 999V999 VALUE 0.0. |
|||
01 Std-Dev PIC 999V999 VALUE 0.0. |
|||
01 Seed PIC 999V999. |
|||
01 TI PIC 9(8). |
|||
01 Idx PIC 99999 VALUE 0. |
|||
01 Intermediate PIC 9(10)V9(5) VALUE 0.0. |
|||
01 Rnd-Work. |
|||
05 Rnd-Tbl |
|||
OCCURS 1 TO 99999 TIMES DEPENDING ON Sample-Size. |
|||
10 Rnd PIC 9V9999999 VALUE 0.0. |
|||
PROCEDURE DIVISION. |
|||
Main-Program. |
|||
ACCEPT TI FROM TIME. |
|||
MOVE FUNCTION RANDOM(TI) TO Seed. |
|||
PERFORM WITH TEST AFTER VARYING Idx |
|||
FROM 1 BY 1 |
|||
UNTIL Idx = Sample-Size |
|||
COMPUTE Intermediate = |
|||
(FUNCTION RANDOM() * 2.01) |
|||
MOVE Intermediate TO Rnd(Idx) |
|||
END-PERFORM. |
|||
PERFORM WITH TEST AFTER VARYING Idx |
|||
FROM 1 BY 1 |
|||
UNTIL Idx = Sample-Size |
|||
COMPUTE Total = Total + Rnd(Idx) |
|||
END-PERFORM. |
|||
COMPUTE Arith-Mean = Total / Sample-Size. |
|||
DISPLAY "Mean: " Arith-Mean. |
|||
PERFORM WITH TEST AFTER VARYING Idx |
|||
FROM 1 BY 1 |
|||
UNTIL Idx = Sample-Size |
|||
COMPUTE Intermediate = |
|||
Intermediate + (Rnd(Idx) - Arith-Mean) ** 2 |
|||
END-PERFORM. |
|||
COMPUTE Std-Dev = Intermediate / Sample-Size. |
|||
DISPLAY "Std-Dev: " Std-Dev. |
|||
STOP RUN. |
|||
END PROGRAM RANDOM. |
|||
</syntaxhighlight> |
|||
=={{header|Common Lisp}}== |
=={{header|Common Lisp}}== |
||
< |
<syntaxhighlight 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)))</ |
collect (1+ (* (sqrt (* -2 (log (random 1.0)))) (cos (* 2 pi (random 1.0))) 0.5)))</syntaxhighlight> |
||
=={{header|Crystal}}== |
|||
{{trans|Phix}} |
|||
<syntaxhighlight lang="ruby">n, mean, sd, tau = 1000, 1, 0.5, (2 * Math::PI) |
|||
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 ) |
|||
puts "mean = #{mean}, standard deviation = #{standev}"</syntaxhighlight> |
|||
{{out}} |
|||
<pre>mean = 1.0093442539237896, standard deviation = 0.504694489463623</pre> |
|||
=={{header|D}}== |
=={{header|D}}== |
||
< |
<syntaxhighlight lang="d">import std.stdio, std.random, std.math; |
||
struct NormalRandom { |
struct NormalRandom { |
||
double mean, stdDev; |
double mean, stdDev; |
||
// |
// Necessary because it also defines an opCall. |
||
this(in double mean_, in double stdDev_) pure nothrow { |
this(in double mean_, in double stdDev_) pure nothrow { |
||
this.mean = mean_; |
this.mean = mean_; |
||
Line 272: | Line 754: | ||
double opCall() const /*nothrow*/ { |
double opCall() const /*nothrow*/ { |
||
immutable |
immutable r1 = uniform01, r2 = uniform01; // Not nothrow. |
||
return mean + stdDev * sqrt(-2 * r1.log) * cos(2 * PI * r2); |
|||
return mean + stdDev * sqrt(-2 * log(r1)) * cos(2 * PI * r2); |
|||
} |
} |
||
} |
} |
||
Line 280: | Line 761: | ||
void main() { |
void main() { |
||
double[1000] array; |
double[1000] array; |
||
auto |
auto nRnd = NormalRandom(1.0, 0.5); |
||
foreach (ref x; array) |
foreach (ref x; array) |
||
x = |
//x = nRnd; |
||
x = nRnd(); |
|||
}</lang> |
|||
}</syntaxhighlight> |
|||
===Alternative Version=== |
===Alternative Version=== |
||
(Untested) |
(Untested) |
||
{{libheader|tango}} |
{{libheader|tango}} |
||
< |
<syntaxhighlight lang="d">import tango.math.random.Random; |
||
void main() { |
void main() { |
||
Line 297: | Line 780: | ||
l = 1.0 + 0.5 * l; |
l = 1.0 + 0.5 * l; |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
=={{header|Delphi}}== |
=={{header|Delphi}}== |
||
Line 303: | Line 786: | ||
Delphi has RandG function which generates random numbers with normal distribution using Marsaglia-Bray algorithm: |
Delphi has RandG function which generates random numbers with normal distribution using Marsaglia-Bray algorithm: |
||
< |
<syntaxhighlight lang="delphi">program Randoms; |
||
{$APPTYPE CONSOLE} |
{$APPTYPE CONSOLE} |
||
Line 321: | Line 804: | ||
Writeln('Std Deviation = ', StdDev(Values):6:4); |
Writeln('Std Deviation = ', StdDev(Values):6:4); |
||
Readln; |
Readln; |
||
end.</ |
end.</syntaxhighlight> |
||
{{out}} |
|||
Output: |
|||
<pre>Mean = 1.0098 |
<pre>Mean = 1.0098 |
||
Std deviation = 0.5016</pre> |
Std deviation = 0.5016</pre> |
||
=={{header|DWScript}}== |
=={{header|DWScript}}== |
||
< |
<syntaxhighlight lang="delphi">var values : array [0..999] of Float; |
||
var i : Integer; |
var i : Integer; |
||
for i := values.Low to values.High do |
for i := values.Low to values.High do |
||
values := RandG(1, 0.5);</ |
values[i] := RandG(1, 0.5);</syntaxhighlight> |
||
=={{header|E}}== |
=={{header|E}}== |
||
< |
<syntaxhighlight lang="e">accum [] for _ in 1..1000 { _.with(entropy.nextGaussian()) }</syntaxhighlight> |
||
=={{header|EasyLang}}== |
|||
<syntaxhighlight lang="text"> |
|||
numfmt 5 0 |
|||
e = 2.7182818284590452354 |
|||
for i = 1 to 1000 |
|||
a[] &= 1 + 0.5 * sqrt (-2 * log10 randomf / log10 e) * cos (360 * randomf) |
|||
. |
|||
for v in a[] |
|||
avg += v / len a[] |
|||
. |
|||
print "Average: " & avg |
|||
for v in a[] |
|||
s += pow (v - avg) 2 |
|||
. |
|||
s = sqrt (s / len a[]) |
|||
print "Std deviation: " & s |
|||
</syntaxhighlight> |
|||
=={{header|Eiffel}}== |
=={{header|Eiffel}}== |
||
< |
<syntaxhighlight lang="eiffel"> |
||
class |
class |
||
APPLICATION |
APPLICATION |
||
Line 433: | Line 936: | ||
end |
end |
||
end |
end |
||
</syntaxhighlight> |
|||
</lang> |
|||
Example Result |
Example Result |
||
Line 439: | Line 942: | ||
Average: 1.0079398405028137 |
Average: 1.0079398405028137 |
||
Standard Deviation: 0.49042787564453988 |
Standard Deviation: 0.49042787564453988 |
||
</pre> |
|||
=={{header|Elena}}== |
|||
{{trans|C#}} |
|||
ELENA 6.x : |
|||
<syntaxhighlight lang="elena">import extensions; |
|||
import extensions'math; |
|||
randomNormal() |
|||
{ |
|||
^ cos(2 * Pi_value * randomGenerator.nextReal()) |
|||
* sqrt(-2 * ln(randomGenerator.nextReal())) |
|||
} |
|||
public program() |
|||
{ |
|||
real[] a := new real[](1000); |
|||
real tAvg := 0; |
|||
for (int x := 0; x < a.Length; x += 1) |
|||
{ |
|||
a[x] := (randomNormal()) / 2 + 1; |
|||
tAvg += a[x] |
|||
}; |
|||
tAvg /= a.Length; |
|||
console.printLine("Average: ", tAvg); |
|||
real s := 0; |
|||
for (int x := 0; x < a.Length; x += 1) |
|||
{ |
|||
s += power(a[x] - tAvg, 2) |
|||
}; |
|||
s := sqrt(s / 1000); |
|||
console.printLine("Standard Deviation: ", s); |
|||
console.readChar() |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Average: 0.9842420481571 |
|||
Standard Deviation: 0.5109070975558 |
|||
</pre> |
|||
=={{header|Elixir}}== |
|||
<syntaxhighlight lang="elixir">defmodule Random do |
|||
def normal(mean, sd) do |
|||
{a, b} = {:rand.uniform, :rand.uniform} |
|||
mean + sd * (:math.sqrt(-2 * :math.log(a)) * :math.cos(2 * :math.pi * b)) |
|||
end |
|||
end |
|||
std_dev = fn (list) -> |
|||
mean = Enum.sum(list) / length(list) |
|||
sd = Enum.reduce(list, 0, fn x,acc -> acc + (x-mean)*(x-mean) end) / length(list) |
|||
|> :math.sqrt |
|||
IO.puts "Mean: #{mean},\tStdDev: #{sd}" |
|||
end |
|||
xs = for _ <- 1..1000, do: Random.normal(1.0, 0.5) |
|||
std_dev.(xs)</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Mean: 1.009079383094275, StdDev: 0.4991894476975088 |
|||
</pre> |
|||
used Erlang function <code>:rand.normal</code> |
|||
<syntaxhighlight lang="elixir">xs = for _ <- 1..1000, do: 1.0 + :rand.normal * 0.5 |
|||
std_dev.(xs)</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Mean: 0.9955701150615597, StdDev: 0.5036412260426065 |
|||
</pre> |
</pre> |
||
Line 444: | Line 1,021: | ||
{{works with|Erlang}} |
{{works with|Erlang}} |
||
< |
<syntaxhighlight lang="erlang"> |
||
mean(Values) -> |
mean(Values) -> |
||
mean(tl(Values), hd(Values), 1). |
mean(tl(Values), hd(Values), 1). |
||
Line 476: | Line 1,053: | ||
io:format("mean = ~w\n", [mean(X)]), |
io:format("mean = ~w\n", [mean(X)]), |
||
io:format("stddev = ~w\n", [stddev(X)]). |
io:format("stddev = ~w\n", [stddev(X)]). |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
|||
Output: |
|||
<pre> |
<pre> |
||
mean = 1.0118289913718608 |
mean = 1.0118289913718608 |
||
stddev = 0.5021636849524854 |
stddev = 0.5021636849524854 |
||
</pre> |
</pre> |
||
=={{header|ERRE}}== |
|||
<syntaxhighlight lang="text"> |
|||
PROGRAM DISTRIBUTION |
|||
! |
|||
! for rosettacode.org |
|||
! |
|||
! formulas taken from TI-59 Master Library manual |
|||
CONST NUM_ITEM=1000 |
|||
!VAR SUMX#,SUMX2#,R1#,R2#,Z#,I% |
|||
DIM A#[1000] |
|||
BEGIN |
|||
! seeds random number generator with system time |
|||
RANDOMIZE(TIMER) |
|||
PRINT(CHR$(12);) !CLS |
|||
SUMX#=0 SUMX2#=0 |
|||
FOR I%=1 TO NUM_ITEM DO |
|||
R1#=RND(1) R2#=RND(1) |
|||
Z#=SQR(-2*LOG(R1#))*COS(2*π*R2#) |
|||
A#[I%]=Z#/2+1 ! I want a normal distribution with |
|||
! mean=1 and std.dev=0.5 |
|||
SUMX#+=A#[I%] SUMX2#+=A#[I%]*A#[I%] |
|||
END FOR |
|||
Z#=SUMX#/NUM_ITEM |
|||
PRINT("Average is";Z#) |
|||
PRINT("Standard dev. is";SQR(SUMX2#/NUM_ITEM-Z#*Z#)) |
|||
END PROGRAM |
|||
</syntaxhighlight> |
|||
=={{header|Euler Math Toolbox}}== |
=={{header|Euler Math Toolbox}}== |
||
<syntaxhighlight lang="euler math toolbox"> |
|||
<lang Euler Math Toolbox> |
|||
>v=normal(1,1000)*0.5+1; |
>v=normal(1,1000)*0.5+1; |
||
>mean(v), dev(v) |
>mean(v), dev(v) |
||
1.00291801071 |
1.00291801071 |
||
0.498226876528 |
0.498226876528 |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Euphoria}}== |
=={{header|Euphoria}}== |
||
{{trans|PureBasic}} |
{{trans|PureBasic}} |
||
< |
<syntaxhighlight lang="euphoria">include misc.e |
||
function RandomNormal() |
function RandomNormal() |
||
Line 508: | Line 1,124: | ||
for i = 1 to n do |
for i = 1 to n do |
||
s[i] = 1 + 0.5 * RandomNormal() |
s[i] = 1 + 0.5 * RandomNormal() |
||
end for</ |
end for</syntaxhighlight> |
||
=={{header|F_Sharp|F#}}== |
|||
<syntaxhighlight lang="fsharp"> |
|||
let n = MathNet.Numerics.Distributions.Normal(1.0,0.5) |
|||
List.init 1000 (fun _->n.Sample()) |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
[0.734433576; 1.54225304; 0.4407528678; 1.177675412; 0.4318617021; |
|||
0.6026656337; 0.769764924; 1.104693934; 0.6297500925; 0.9594598077; |
|||
1.684736389; 1.160376323; 0.883354356; 0.9513968363; 0.9727698268; |
|||
0.5315570949; 0.9599239266; 1.564976755; 0.7232002879; 1.084139442; |
|||
1.220914517; 0.3553085946; 1.112549824; 1.989443553; 0.5752307543; |
|||
1.156682549; 0.7886670467; 0.02050745923; 1.532060208; 1.18789591; |
|||
1.408946777; 1.038812004; 1.724679503; 1.671565045; 1.266831442; |
|||
1.363611654; 1.705819067; 0.5772366328; 0.4503488498; 1.496891481; |
|||
0.9831877282; 0.3845460366; 0.8253240671; 1.298969969; 0.4265904553; |
|||
0.9303696876; 0.445003361; 0.753175816; 0.6143534043; 1.059982235; |
|||
0.7143206784; 0.2233328038; 1.005178481; 0.7697392436; 0.5904948577; |
|||
0.5127953044; 0.6467346747; 0.7929387604; -0.1501790761; 0.8750780903; |
|||
0.941704369; 1.37941579; 0.4739006145; 1.998886344; 1.219428519; |
|||
0.06270791476; 1.097739804; 0.7584232803; 1.042177217; 1.166561247; |
|||
1.502357164; 1.171525776; 0.1528807432; 0.2289389756; 1.36208422; |
|||
0.3714421124; 1.299571092; 1.171553369; 1.317807265; 1.616662281; |
|||
1.724223246; 1.059580642; 1.270520918; -0.1827677907; 1.938593232; |
|||
1.420362143; 1.888357595; 0.7851629936; 0.7080554899; 0.7747215818; |
|||
1.403719877; 0.5765950249; 1.275206565; 0.6292054813; 1.525562798; |
|||
0.6224640457; 0.8524078517; 0.7646595627; 0.6799834691; 0.773111053; ...] |
|||
</pre> =={{header|F_Sharp|F#}}== |
|||
<syntaxhighlight 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() ]</syntaxhighlight> |
|||
=={{header|Factor}}== |
=={{header|Factor}}== |
||
<lang |
<syntaxhighlight lang="factor">USING: random ; |
||
1000 [ 1.0 0.5 normal-random-float ] replicate</syntaxhighlight> |
|||
=={{header|Falcon |
=={{header|Falcon}}== |
||
< |
<syntaxhighlight lang="falcon">a = [] |
||
for i in [0:1000] : a+= norm_rand_num() |
for i in [0:1000] : a+= norm_rand_num() |
||
Line 520: | Line 1,172: | ||
pi = 2*acos(0) |
pi = 2*acos(0) |
||
return 1 + (cos(2 * pi * random()) * pow(-2 * log(random()) ,1/2)) /2 |
return 1 + (cos(2 * pi * random()) * pow(-2 * log(random()) ,1/2)) /2 |
||
end</ |
end</syntaxhighlight> |
||
=={{header|Fantom}}== |
=={{header|Fantom}}== |
||
Line 526: | Line 1,178: | ||
Two solutions. The first uses Fantom's random-number generator, which produces a uniform distribution. So, convert to a normal distribution using a formula: |
Two solutions. The first uses Fantom's random-number generator, which produces a uniform distribution. So, convert to a normal distribution using a formula: |
||
< |
<syntaxhighlight lang="fantom"> |
||
class Main |
class Main |
||
{ |
{ |
||
Line 544: | Line 1,196: | ||
} |
} |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
The second calls out to Java's Gaussian random-number generator: |
The second calls out to Java's Gaussian random-number generator: |
||
< |
<syntaxhighlight lang="fantom"> |
||
using [java] java.util::Random |
using [java] java.util::Random |
||
Line 569: | Line 1,221: | ||
} |
} |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Forth}}== |
=={{header|Forth}}== |
||
{{works with|gforth|0.6.2}} |
{{works with|gforth|0.6.2}} |
||
< |
<syntaxhighlight lang="forth">require random.fs |
||
here to seed |
here to seed |
||
Line 589: | Line 1,241: | ||
0 do frnd-normal 0.5e f* 1e f+ f, loop ; |
0 do frnd-normal 0.5e f* 1e f+ f, loop ; |
||
create rnd-array 1000 ,normals</ |
create rnd-array 1000 ,normals</syntaxhighlight> |
||
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. |
|||
<syntaxhighlight lang="text">rnd rnd dabs d>f</syntaxhighlight> is necessary, but surprising and definitely not well documented / perhaps not compliant. |
|||
=={{header|Fortran}}== |
=={{header|Fortran}}== |
||
{{works with|Fortran|90 and later}} |
{{works with|Fortran|90 and later}} |
||
< |
<syntaxhighlight lang="fortran">PROGRAM Random |
||
INTEGER, PARAMETER :: n = 1000 |
INTEGER, PARAMETER :: n = 1000 |
||
Line 616: | Line 1,272: | ||
WRITE(*, "(A,F8.6)") "Standard Deviation = ", sd |
WRITE(*, "(A,F8.6)") "Standard Deviation = ", sd |
||
END PROGRAM Random</ |
END PROGRAM Random</syntaxhighlight> |
||
{{out}} |
|||
Output |
|||
<pre> |
|||
Mean = 0.995112 |
Mean = 0.995112 |
||
Standard Deviation = 0.503373 |
Standard Deviation = 0.503373 |
||
</pre> |
|||
=={{header| |
=={{header|Free Pascal}}== |
||
Free Pascal provides the '''randg''' function in the RTL math unit that produces Gaussian-distributed random numbers with the Box-Müller algorithm. |
|||
<lang fsharp>let gaussianRand count = |
|||
let o = new System.Random() |
|||
<syntaxhighlight lang="pascal"> |
|||
let pi = System.Math.PI |
|||
function randg(mean,stddev: float): float; |
|||
let gaussrnd = |
|||
</syntaxhighlight> |
|||
(fun _ -> 1. + 0.5 * sqrt(-2. * log(o.NextDouble())) * cos(2. * pi * o.NextDouble())) |
|||
[ for i in {0 .. (int count)} -> gaussrnd() ]</lang> |
|||
=={{header|Frink}}== |
|||
<syntaxhighlight lang="frink">a = new array[[1000], {|x| randomGaussian[1, 0.5]}]</syntaxhighlight> |
|||
=={{header|Go}}== |
=={{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. |
|||
<lang go>package main |
|||
<syntaxhighlight lang="go">package main |
|||
import ( |
import ( |
||
"fmt" |
|||
"math" |
|||
"math/rand" |
"math/rand" |
||
"strings" |
|||
"time" |
"time" |
||
) |
) |
||
Line 640: | Line 1,304: | ||
const mean = 1.0 |
const mean = 1.0 |
||
const stdv = .5 |
const stdv = .5 |
||
const n = 1000 |
|||
func main() { |
func main() { |
||
var list [ |
var list [n]float64 |
||
rand.Seed(time.Now().UnixNano()) |
rand.Seed(time.Now().UnixNano()) |
||
for i := range list { |
for i := range list { |
||
list[i] = mean + stdv*rand.NormFloat64() |
list[i] = mean + stdv*rand.NormFloat64() |
||
} |
} |
||
// show computed mean and stdv of list |
|||
}</lang> |
|||
var s, sq float64 |
|||
for _, v := range list { |
|||
s += v |
|||
} |
|||
cm := s / n |
|||
for _, v := range list { |
|||
d := v - cm |
|||
sq += d * d |
|||
} |
|||
fmt.Printf("mean %.3f, stdv %.3f\n", cm, math.Sqrt(sq/(n-1))) |
|||
// show histogram by hdiv divisions per stdv over +/-hrange stdv |
|||
const hdiv = 3 |
|||
const hrange = 2 |
|||
var h [1 + 2*hrange*hdiv]int |
|||
for _, v := range list { |
|||
bin := hrange*hdiv + int(math.Floor((v-mean)/stdv*hdiv+.5)) |
|||
if bin >= 0 && bin < len(h) { |
|||
h[bin]++ |
|||
} |
|||
} |
|||
const hscale = 10 |
|||
for _, c := range h { |
|||
fmt.Println(strings.Repeat("*", (c+hscale/2)/hscale)) |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
mean 0.995, stdv 0.503 |
|||
** |
|||
**** |
|||
****** |
|||
******** |
|||
************ |
|||
************ |
|||
************* |
|||
************ |
|||
********** |
|||
******** |
|||
***** |
|||
*** |
|||
** |
|||
</pre> |
|||
=={{header|Groovy}}== |
=={{header|Groovy}}== |
||
< |
<syntaxhighlight lang="groovy">rnd = new Random() |
||
result = (1..1000).inject([]) { r, i -> r << rnd.nextGaussian() }</ |
result = (1..1000).inject([]) { r, i -> r << rnd.nextGaussian() }</syntaxhighlight> |
||
=={{header|Haskell}}== |
=={{header|Haskell}}== |
||
< |
<syntaxhighlight lang="haskell">import System.Random |
||
pairs :: [a] -> [(a,a)] |
pairs :: [a] -> [(a,a)] |
||
Line 668: | Line 1,375: | ||
result :: IO [Double] |
result :: IO [Double] |
||
result = getStdGen >>= \g -> return $ gaussians 1000 g</ |
result = getStdGen >>= \g -> return $ gaussians 1000 g</syntaxhighlight> |
||
Or using Data.Random from random-fu package: |
|||
<syntaxhighlight lang="haskell">replicateM 1000 $ normal 1 0.5</syntaxhighlight> |
|||
To print them: |
|||
<syntaxhighlight lang="haskell">import Data.Random |
|||
import Control.Monad |
|||
thousandRandomNumbers :: RVar [Double] |
|||
thousandRandomNumbers = replicateM 1000 $ normal 1 0.5 |
|||
main = do |
|||
x <- sample thousandRandomNumbers |
|||
print x</syntaxhighlight> |
|||
=={{header|HicEst}}== |
=={{header|HicEst}}== |
||
< |
<syntaxhighlight lang="hicest">REAL :: n=1000, m=1, s=0.5, array(n) |
||
pi = 4 * ATAN(1) |
pi = 4 * ATAN(1) |
||
array = s * (-2*LOG(RAN(1)))^0.5 * COS(2*pi*RAN(1)) + m </ |
array = s * (-2*LOG(RAN(1)))^0.5 * COS(2*pi*RAN(1)) + m </syntaxhighlight> |
||
=={{header|Icon}} and {{header|Unicon}}== |
=={{header|Icon}} and {{header|Unicon}}== |
||
Line 680: | Line 1,400: | ||
Note that Unicon randomly seeds it's generator. |
Note that Unicon randomly seeds it's generator. |
||
< |
<syntaxhighlight lang="icon"> |
||
procedure main() |
procedure main() |
||
local L |
local L |
||
Line 688: | Line 1,408: | ||
every write(!L) |
every write(!L) |
||
end |
end |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|IDL}}== |
=={{header|IDL}}== |
||
< |
<syntaxhighlight lang="idl">result = 1.0 + 0.5*randomn(seed,1000)</syntaxhighlight> |
||
=={{header|J}}== |
=={{header|J}}== |
||
'''Solution:''' |
'''Solution:''' |
||
< |
<syntaxhighlight lang="j">urand=: ?@$ 0: |
||
zrand=: (2 o. 2p1 * urand) * [: %: _2 * [: ^. urand |
zrand=: (2 o. 2p1 * urand) * [: %: _2 * [: ^. urand |
||
1 + 0.5 * zrand 100</ |
1 + 0.5 * zrand 100</syntaxhighlight> |
||
'''Alternative Solution:'''<br> |
'''Alternative Solution:'''<br> |
||
Using the normal script from the [[j:Addons/stats/distribs|stats/distribs addon]]. |
Using the normal script from the [[j:Addons/stats/distribs|stats/distribs addon]]. |
||
< |
<syntaxhighlight lang="j"> require 'stats/distribs/normal' |
||
1 0.5 rnorm 1000 |
1 0.5 rnorm 1000 |
||
1.44868803 1.21548637 0.812460657 1.54295452 1.2470606 ...</ |
1.44868803 1.21548637 0.812460657 1.54295452 1.2470606 ...</syntaxhighlight> |
||
=={{header|Java}}== |
=={{header|Java}}== |
||
< |
<syntaxhighlight lang="java">double[] list = new double[1000]; |
||
double mean = 1.0, std = 0.5; |
double mean = 1.0, std = 0.5; |
||
Random rng = new Random(); |
Random rng = new Random(); |
||
for(int i = 0;i<list.length;i++) { |
for(int i = 0;i<list.length;i++) { |
||
list[i] = mean + std * rng.nextGaussian(); |
list[i] = mean + std * rng.nextGaussian(); |
||
}</ |
}</syntaxhighlight> |
||
=={{header|JavaScript}}== |
=={{header|JavaScript}}== |
||
< |
<syntaxhighlight lang="javascript">function randomNormal() { |
||
return Math.cos(2 * Math.PI * Math.random()) * Math.sqrt(-2 * Math.log(Math.random())) |
return Math.cos(2 * Math.PI * Math.random()) * Math.sqrt(-2 * Math.log(Math.random())) |
||
} |
} |
||
Line 722: | Line 1,442: | ||
for (var i=0; i < 1000; i++){ |
for (var i=0; i < 1000; i++){ |
||
a[i] = randomNormal() / 2 + 1 |
a[i] = randomNormal() / 2 + 1 |
||
}</ |
}</syntaxhighlight> |
||
=={{header|jq}}== |
|||
{{works with|jq|1.4}} |
|||
Since jq is a purely functional language, it is convenient |
|||
to define the pseudo-random number generator functions as filters |
|||
whose inputs and outputs are arrays containing a "seed". |
|||
The following uses the same pseudo-random number generator as the |
|||
Microsoft C Runtime (see [[Linear congruential generator]]). |
|||
''''A Pseudo-Random Number Generator'''' |
|||
<syntaxhighlight lang="jq"># 15-bit integers generated using the same formula as rand() from the Microsoft C Runtime. |
|||
# The random numbers are in [0 -- 32767] inclusive. |
|||
# Input: an array of length at least 2 interpreted as [count, state, ...] |
|||
# Output: [count+1, newstate, r] where r is the next pseudo-random number. |
|||
def next_rand_Microsoft: |
|||
.[0] as $count | .[1] as $state |
|||
| ( (214013 * $state) + 2531011) % 2147483648 # mod 2^31 |
|||
| [$count+1 , ., (. / 65536 | floor) ] ;</syntaxhighlight> |
|||
''''Box-Muller Method'''' |
|||
<syntaxhighlight lang="jq"># Generate a single number following the normal distribution with mean 0, variance 1, |
|||
# 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] |
|||
# Output [n+1, nextstate, r] |
|||
def next_rand_normal: |
|||
def u: next_rand_Microsoft | .[2] /= 32767; |
|||
u as $u1 |
|||
| ($u1 | u) as $u2 |
|||
| ((( (8*(1|atan)) * $u1[2]) | cos) |
|||
* ((-2 * (($u2[2]) | log)) | sqrt)) as $r |
|||
| [ (.[0]+1), $u2[1], $r] ; |
|||
# Generate "count" arrays, each containing a random normal variate with the given mean and standard deviation. |
|||
# Input: [count, state] |
|||
# Output: [updatedcount, updatedstate, rnv] |
|||
# where "state" is a seed and "updatedstate" can be used as a seed. |
|||
def random_normal_variate(mean; sd; count): |
|||
next_rand_normal |
|||
| recurse( if .[0] < count then next_rand_normal else empty end) |
|||
| .[2] = (.[2] * sd) + mean;</syntaxhighlight> |
|||
'''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: |
|||
<syntaxhighlight lang="jq">def summary: |
|||
length as $l | add as $sum | ($sum/$l) as $a |
|||
| reduce .[] as $x (0; . + ( ($x - $a) | .*. )) |
|||
| [ $a, (./$l | sqrt)] ; |
|||
[ [0,1] | random_normal_variate(1; 0.5; 1000) | .[2] ] | summary</syntaxhighlight> |
|||
{{out}} |
|||
$ jq -n -c -f Random_numbers.jq |
|||
[0.9932830741018853,0.4977760644490579] |
|||
=={{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): |
|||
<syntaxhighlight lang="julia">randn(1000) * 0.5 + 1</syntaxhighlight> |
|||
=={{header|Kotlin}}== |
|||
<syntaxhighlight lang="scala">// version 1.0.6 |
|||
import java.util.Random |
|||
fun main(args: Array<String>) { |
|||
val r = Random() |
|||
val da = DoubleArray(1000) |
|||
for (i in 0 until 1000) da[i] = 1.0 + 0.5 * r.nextGaussian() |
|||
// now check actual mean and SD |
|||
val mean = da.average() |
|||
val sd = Math.sqrt(da.map { (it - mean) * (it - mean) }.average()) |
|||
println("Mean is $mean") |
|||
println("S.D. is $sd") |
|||
}</syntaxhighlight> |
|||
Sample output: |
|||
{{out}} |
|||
<pre> |
|||
Mean is 1.0071784073168768 |
|||
S.D. is 0.48567118114896807 |
|||
</pre> |
|||
=={{header|LabVIEW}}== |
=={{header|LabVIEW}}== |
||
Line 728: | Line 1,528: | ||
[[File:LV_array_of_randoms_with_given_mean_and_stdev.png]] |
[[File:LV_array_of_randoms_with_given_mean_and_stdev.png]] |
||
=={{header| |
=={{header|Lingo}}== |
||
<syntaxhighlight lang="lingo">-- Returns a random float value in range 0..1 |
|||
<lang lb>dim a(1000) |
|||
on randf () |
|||
mean =1 |
|||
n = random(the maxinteger)-1 |
|||
sd =0.5 |
|||
return n / float(the maxinteger-1) |
|||
for i = 1 to 1000 ' throw 1000 normal variates |
|||
end</syntaxhighlight> |
|||
a( i) =mean +sd *( sqr( -2 * log( rnd( 0))) * cos( 2 * pi * rnd( 0))) |
|||
next i</lang> |
|||
<syntaxhighlight lang="lingo">normal = [] |
|||
repeat with i = 1 to 1000 |
|||
normal.add(1 + sqrt(-2 * log(randf())) * cos(2 * PI * randf()) / 2) |
|||
end repeat</syntaxhighlight> |
|||
=={{header|Lobster}}== |
|||
Uses built-in <code>rnd_gaussian</code> |
|||
<syntaxhighlight lang="lobster"> |
|||
let mean = 1.0 |
|||
let stdv = 0.5 |
|||
let count = 1000 |
|||
// stats computes a running mean and variance |
|||
// See Knuth TAOCP vol 2, 3rd edition, page 232 |
|||
def stats(xs: [float]) -> float, float: // variance, mean |
|||
var M = xs[0] |
|||
var S = 0.0 |
|||
var n = 1.0 |
|||
for(xs.length - 1) i: |
|||
let x = xs[i + 1] |
|||
n = n + 1.0 |
|||
let mm = (x - M) |
|||
M += mm / n |
|||
S += mm * (x - M) |
|||
return (if n > 0.0: S / n else: 0.0), M |
|||
def test_random_normal() -> [float]: |
|||
rnd_seed(floor(seconds_elapsed() * 1000000)) |
|||
let r = vector_reserve(typeof return, count) |
|||
for (count): |
|||
r.push(rnd_gaussian() * stdv + mean) |
|||
let cvar, cmean = stats(r) |
|||
let cstdv = sqrt(cvar) |
|||
print concat_string(["Mean: ", string(cmean), ", Std.Deviation: ", string(cstdv)], "") |
|||
test_random_normal() |
|||
</syntaxhighlight> |
|||
=={{header|Logo}}== |
=={{header|Logo}}== |
||
{{works with|UCB 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. |
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. |
||
< |
<syntaxhighlight lang="logo">to random.float ; 0..1 |
||
localmake "max.int lshift -1 -1 |
localmake "max.int lshift -1 -1 |
||
output quotient random :max.int :max.int |
output quotient random :max.int :max.int |
||
Line 748: | Line 1,586: | ||
end |
end |
||
make "randoms cascade 1000 [fput random.gaussian / 2 + 1 ?] []</ |
make "randoms cascade 1000 [fput random.gaussian / 2 + 1 ?] []</syntaxhighlight> |
||
=={{header|Lua}}== |
|||
<syntaxhighlight lang="lua">local list = {} |
|||
=={{header|lua}}== |
|||
<lang lua> |
|||
local list = {} |
|||
for i = 1, 1000 do |
for i = 1, 1000 do |
||
list[i] = 1 + math.sqrt(-2 * math.log(math.random())) * math.cos(2 * math.pi * math.random()) / 2 |
list[i] = 1 + math.sqrt(-2 * math.log(math.random())) * math.cos(2 * math.pi * math.random()) / 2 |
||
end</syntaxhighlight> |
|||
end |
|||
</lang> |
|||
=={{header|M2000 Interpreter}}== |
|||
M2000 use a Wichmann - Hill Pseudo Random Number Generator. |
|||
<syntaxhighlight lang="m2000 interpreter"> |
|||
Module CheckIt { |
|||
Function StdDev (A()) { |
|||
\\ A() has a copy of values |
|||
N=Len(A()) |
|||
if N<1 then Error "Empty Array" |
|||
M=Each(A()) |
|||
k=0 |
|||
\\ make sum, dev same type as A(k) |
|||
sum=A(k)-A(k) |
|||
dev=sum |
|||
\\ find mean |
|||
While M { |
|||
sum+=Array(M) |
|||
} |
|||
Mean=sum/N |
|||
\\ make a pointet to A() |
|||
P=A() |
|||
\\ subtruct from each item |
|||
P-=Mean |
|||
M=Each(P) |
|||
While M { |
|||
dev+=Array(M)*Array(M) |
|||
} |
|||
\\ as pointer to arrray |
|||
=(if(dev>0->Sqrt(dev/N), 0), Mean) |
|||
} |
|||
Function randomNormal { |
|||
\\ by default all numbers are double |
|||
\\ cos() get degrees |
|||
=1+Cos(360 * rnd) * Sqrt(-2 * Ln(rnd)) /2 |
|||
} |
|||
\\ fill array calling randomNormal() for each item |
|||
Dim A(1000)<<randomNormal() |
|||
\\ we can pass a pointer to array and place it to stack of values |
|||
DisplayMeanAndStdDeviation(A()) ' mean ~ 1 std deviation ~0.5 |
|||
\\ check M2000 rnd only |
|||
Dim B(1000)<<rnd |
|||
DisplayMeanAndStdDeviation(B()) ' mean ~ 0.5 std deviation ~0.28 |
|||
DisplayMeanAndStdDeviation((0,0,14,14)) ' mean = 7 std deviation = 7 |
|||
DisplayMeanAndStdDeviation((0,6,8,14)) ' mean = 7 std deviation = 5 |
|||
DisplayMeanAndStdDeviation((6,6,8,8)) ' mean = 7 std deviation = 1 |
|||
Sub DisplayMeanAndStdDeviation(A) |
|||
\\ push to stack all items of an array (need an array pointer) |
|||
Push ! StdDev(A) |
|||
\\ read from strack two numbers |
|||
Print "Mean is "; Number |
|||
Print "Standard Deviation is "; Number |
|||
End Sub |
|||
} |
|||
Checkit |
|||
</syntaxhighlight> |
|||
=={{header|Maple}}== |
|||
<syntaxhighlight lang="maple">with(Statistics): |
|||
Sample(Normal(1, 0.5), 1000);</syntaxhighlight> |
|||
'''or''' |
|||
<syntaxhighlight lang="maple">1+0.5*ArrayTools[RandomArray](1000,1,distribution=normal);</syntaxhighlight> |
|||
=={{header|Mathematica}}== |
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
||
Built-in function RandomReal with built-in distribution NormalDistribution as an argument: |
Built-in function RandomReal with built-in distribution NormalDistribution as an argument: |
||
< |
<syntaxhighlight lang="mathematica">RandomReal[NormalDistribution[1, 1/2], 1000]</syntaxhighlight> |
||
=={{header|MATLAB}}== |
=={{header|MATLAB}}== |
||
Native support : |
Native support : |
||
< |
<syntaxhighlight lang="matlab"> mu = 1; sd = 0.5; |
||
x = randn(1000,1) * sd + mu; |
x = randn(1000,1) * sd + mu; |
||
</syntaxhighlight> |
|||
</lang> |
|||
The statistics toolbox provides this function |
The statistics toolbox provides this function |
||
< |
<syntaxhighlight lang="matlab"> x = normrnd(mu, sd, [1000,1]); </syntaxhighlight> |
||
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. |
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. |
||
< |
<syntaxhighlight lang="matlab">function randNum = randNorm(mu0,chi2, sz) |
||
radiusSquared = +Inf; |
radiusSquared = +Inf; |
||
Line 789: | Line 1,691: | ||
randNum = (v .* scaleFactor .* chi2) + mu0; |
randNum = (v .* scaleFactor .* chi2) + mu0; |
||
end</ |
end</syntaxhighlight> |
||
Output: |
Output: |
||
< |
<syntaxhighlight lang="matlab">>> randNorm(1,.5, [1000,1]) |
||
ans = |
ans = |
||
0.693984121077029</ |
0.693984121077029</syntaxhighlight> |
||
=={{header|Maxima}}== |
=={{header|Maxima}}== |
||
< |
<syntaxhighlight lang="maxima">load(distrib)$ |
||
random_normal(1.0, 0.5, 1000);</ |
random_normal(1.0, 0.5, 1000);</syntaxhighlight> |
||
=={{header|MAXScript}}== |
=={{header|MAXScript}}== |
||
< |
<syntaxhighlight lang="maxscript">arr = #() |
||
for i in 1 to 1000 do |
for i in 1 to 1000 do |
||
( |
( |
||
Line 812: | Line 1,714: | ||
c = 1.0 + 0.5 * sqrt (-2*log a) * cos (360*b) -- Maxscript cos takes degrees |
c = 1.0 + 0.5 * sqrt (-2*log a) * cos (360*b) -- Maxscript cos takes degrees |
||
append arr c |
append arr c |
||
)</ |
)</syntaxhighlight> |
||
=={{header|Metafont}}== |
=={{header|Metafont}}== |
||
Line 818: | Line 1,720: | ||
Metafont has <code>normaldeviate</code> which produces pseudorandom normal distributed numbers with mean 0 and variance one. So the following complete the task: |
Metafont has <code>normaldeviate</code> which produces pseudorandom normal distributed numbers with mean 0 and variance one. So the following complete the task: |
||
< |
<syntaxhighlight lang="metafont">numeric col[]; |
||
m := 0; % m holds the mean, for testing purposes |
m := 0; % m holds the mean, for testing purposes |
||
Line 836: | Line 1,738: | ||
show m, s; % and let's show that really they get what we wanted |
show m, s; % and let's show that really they get what we wanted |
||
end</ |
end</syntaxhighlight> |
||
A run gave |
A run gave |
||
Line 845: | Line 1,747: | ||
Assigning a value to the special variable '''randomseed''' will allow to have always |
Assigning a value to the special variable '''randomseed''' will allow to have always |
||
the same sequence of pseudorandom numbers |
the same sequence of pseudorandom numbers |
||
=={{header|MiniScript}}== |
|||
<syntaxhighlight lang="miniscript">randNormal = function(mean=0, stddev=1) |
|||
return mean + sqrt(-2 * log(rnd,2.7182818284)) * cos(2*pi*rnd) * stddev |
|||
end function |
|||
x = [] |
|||
for i in range(1,1000) |
|||
x.push randNormal(1, 0.5) |
|||
end for</syntaxhighlight> |
|||
=={{header|Mirah}}== |
=={{header|Mirah}}== |
||
< |
<syntaxhighlight lang="mirah">import java.util.Random |
||
list = double[999] |
list = double[999] |
||
Line 856: | Line 1,768: | ||
list[i] = mean + std * rng.nextGaussian |
list[i] = mean + std * rng.nextGaussian |
||
end |
end |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|МК-61/52}}== |
=={{header|МК-61/52}}== |
||
<lang>П7 <-> П8 1/x П6 ИП6 П9 СЧ П6 1/x |
<syntaxhighlight lang="text">П7 <-> П8 1/x П6 ИП6 П9 СЧ П6 1/x |
||
ln ИП8 * 2 * КвКор ИП9 2 * пи |
ln ИП8 * 2 * КвКор ИП9 2 * пи |
||
* sin * ИП7 + С/П БП 05</ |
* sin * ИП7 + С/П БП 05</syntaxhighlight> |
||
Input: РY - variance, РX - expectation. |
''Input'': РY - variance, РX - expectation. |
||
Or: |
|||
<syntaxhighlight lang="text">3 10^x П0 ПП 13 2 / 1 + С/П L0 03 С/П |
|||
СЧ lg 2 /-/ * КвКор 2 пи ^ СЧ * * cos * В/О</syntaxhighlight> |
|||
to generate 1000 numbers with a mean of 1.0 and a standard deviation of 0.5. |
|||
=={{header|Modula-3}}== |
=={{header|Modula-3}}== |
||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang="modula3">MODULE Rand EXPORTS Main; |
||
IMPORT Random; |
IMPORT Random; |
||
Line 888: | Line 1,807: | ||
rands[i] := 1.0D0 + 0.5D0 * RandNorm(); |
rands[i] := 1.0D0 + 0.5D0 * RandNorm(); |
||
END; |
END; |
||
END Rand.</ |
END Rand.</syntaxhighlight> |
||
=={{header|Nanoquery}}== |
|||
{{trans|Java}} |
|||
<syntaxhighlight lang="nanoquery">list = {0} * 1000 |
|||
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() |
|||
end</syntaxhighlight> |
|||
=={{header|NetRexx}}== |
=={{header|NetRexx}}== |
||
< |
<syntaxhighlight lang="netrexx">/* NetRexx */ |
||
options replace format comments java crossref symbols nobinary |
options replace format comments java crossref symbols nobinary |
||
Line 946: | Line 1,875: | ||
return |
return |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
|||
'''Output:''' |
|||
<pre> |
<pre> |
||
Mean: 1.0 |
Mean: 1.0 |
||
Line 989: | Line 1,918: | ||
=={{header|NewLISP}}== |
=={{header|NewLISP}}== |
||
<lang |
<syntaxhighlight lang="newlisp">(normal 1 .5 1000)</syntaxhighlight> |
||
=={{header|Nim}}== |
|||
<syntaxhighlight lang="nim">import random, stats, strformat |
|||
var rs: RunningStat |
|||
randomize() |
|||
for _ in 1..5: |
|||
for _ in 1..1000: rs.push gauss(1.0, 0.5) |
|||
echo &"mean: {rs.mean:.5f} stdDev: {rs.standardDeviation:.5f}" |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre>mean: 1.01294 stdDev: 0.49692 |
|||
mean: 1.00262 stdDev: 0.50028 |
|||
mean: 0.99878 stdDev: 0.49662 |
|||
mean: 0.99830 stdDev: 0.49820 |
|||
mean: 1.00658 stdDev: 0.49703</pre> |
|||
=={{header|Objeck}}== |
=={{header|Objeck}}== |
||
< |
<syntaxhighlight lang="objeck">bundle Default { |
||
class RandomNumbers { |
class RandomNumbers { |
||
function : Main(args : String[]) ~ Nil { |
function : Main(args : String[]) ~ Nil { |
||
Line 1,009: | Line 1,956: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
=={{header|OCaml}}== |
=={{header|OCaml}}== |
||
< |
<syntaxhighlight lang="ocaml">let pi = 4. *. atan 1.;; |
||
let random_gaussian () = |
let random_gaussian () = |
||
1. +. sqrt (-2. *. log (Random.float 1.)) *. cos (2. *. pi *. Random.float 1.);; |
1. +. sqrt (-2. *. log (Random.float 1.)) *. cos (2. *. pi *. Random.float 1.);; |
||
let a = Array.init 1000 (fun _ -> random_gaussian ());;</ |
let a = Array.init 1000 (fun _ -> random_gaussian ());;</syntaxhighlight> |
||
=={{header|Octave}}== |
=={{header|Octave}}== |
||
< |
<syntaxhighlight lang="octave">p = normrnd(1.0, 0.5, 1000, 1); |
||
disp(mean(p)); |
disp(mean(p)); |
||
disp(sqrt(sum((p - mean(p)).^2)/numel(p)));</ |
disp(sqrt(sum((p - mean(p)).^2)/numel(p)));</syntaxhighlight> |
||
Output: |
|||
{{out}} |
|||
<pre>1.0209 |
<pre>1.0209 |
||
0.51048</pre> |
0.51048</pre> |
||
=={{header|ooRexx}}== |
|||
{{trans|REXX}} |
|||
===version 1=== |
|||
<syntaxhighlight lang="oorexx">/*REXX pgm gens 1,000 normally distributed #s: mean=1, standard dev.=0.5*/ |
|||
pi=RxCalcPi() /* get value of 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.*/ |
|||
n.g=random(0,1e5)/1e5 /* REXX gens uniform rand integers*/ |
|||
End |
|||
Say ' old mean=' mean() |
|||
Say 'old standard deviation=' stddev() |
|||
Say |
|||
Do j=1 To n-1 By 2 |
|||
m=j+1 |
|||
/*use Box-Muller method */ |
|||
_=sd*RxCalcPower(-2*RxCalcLog(n.j),.5)*RxCalcCos(2*pi*n.m,,'R')+mean |
|||
n.m=sd*RxCalcpower(-2*RxCalcLog(n.j),.5)*RxCalcSin(2*pi*n.m,,'R')+, |
|||
mean /* rand # must be 0???1. */ |
|||
n.j=_ |
|||
End /* j */ |
|||
Say ' new mean=' mean() |
|||
Say 'new standard deviation=' stddev() |
|||
Exit |
|||
mean: |
|||
_=0 |
|||
Do k=1 For n |
|||
_=_+n.k |
|||
End |
|||
Return _/n |
|||
stddev: |
|||
_avg=mean() |
|||
_=0 |
|||
Do k=1 For n |
|||
_=_+(n.k-_avg)**2 |
|||
End |
|||
Return RxCalcPower(_/n,.5) |
|||
:: requires rxmath library</syntaxhighlight> |
|||
{{out}} |
|||
<pre> old mean= 0.49830002 |
|||
old standard deviation= 0.283199568 |
|||
new mean= 1.00377404 |
|||
new standard deviation= 0.501444536 </pre> |
|||
===version 2=== |
|||
Using the nice function names in the algorithm. |
|||
<syntaxhighlight lang="oorexx">/*REXX pgm gens 1,000 normally distributed #s: mean=1, standard dev.=0.5*/ |
|||
pi=RxCalcPi() /* get value of 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.*/ |
|||
n.g=random(0,1e5)/1e5 /* REXX gens uniform rand integers*/ |
|||
End |
|||
Say ' old mean=' mean() |
|||
Say 'old standard deviation=' stddev() |
|||
Say |
|||
Do j=1 To n-1 By 2 |
|||
m=j+1 |
|||
/*use Box-Muller method */ |
|||
_=sd*sqrt(-2*ln(n.j))*cos(2*pi*n.m)+mean |
|||
n.m=sd*sqrt(-2*ln(n.j))*sin(2*pi*n.m)+mean |
|||
n.j=_ |
|||
End |
|||
Say ' new mean=' mean() |
|||
Say 'new standard deviation=' stddev() |
|||
Exit |
|||
mean: |
|||
_=0 |
|||
Do k=1 For n |
|||
_=_+n.k |
|||
End |
|||
Return _/n |
|||
stddev: |
|||
_avg=mean() |
|||
_=0 |
|||
Do k=1 For n |
|||
_=_+(n.k-_avg)**2 |
|||
End |
|||
Return sqrt(_/n) |
|||
sqrt: Return RxCalcSqrt(arg(1)) |
|||
ln: Return RxCalcLog(arg(1)) |
|||
cos: Return RxCalcCos(arg(1),,'R') |
|||
sin: Return RxCalcSin(arg(1),,'R') |
|||
:: requires rxmath library</syntaxhighlight> |
|||
=={{header|PARI/GP}}== |
=={{header|PARI/GP}}== |
||
< |
<syntaxhighlight 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); |
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* |
sqrt(-2*log(u1))*cos(2*Pi*u2) \\ in previous version "u1" instead of "u2" was used --> has given crap distribution |
||
\\ Could easily be extended with a second normal at very little cost. |
\\ Could easily be extended with a second normal at very little cost. |
||
}; |
}; |
||
vector(1000,unused,rnormal()/2+1)</ |
vector(1000,unused,rnormal()/2+1)</syntaxhighlight> |
||
=={{header|Pascal}}== |
=={{header|Pascal}}== |
||
See [[Random_numbers#Delphi | Delphi]] |
|||
The following function calculates Gaussian-distributed random numbers with the Box-Müller algorithm: |
|||
<syntaxhighlight lang="pascal"> |
|||
function rnorm (mean, sd: real): real; |
|||
{Calculates Gaussian random numbers according to the Box-Müller approach} |
|||
var |
|||
u1, u2: real; |
|||
begin |
|||
u1 := random; |
|||
u2 := random; |
|||
rnorm := mean * abs(1 + sqrt(-2 * (ln(u1))) * cos(2 * pi * u2) * sd); |
|||
/* error !?! Shouldn't it be "mean +" instead of "mean *" ? */ |
|||
end; |
|||
</syntaxhighlight> |
|||
[[#Delphi | Delphi]] and [[#Free Pascal|Free Pascal]] support implement a '''randg''' function that delivers Gaussian-distributed random numbers. |
|||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
< |
<syntaxhighlight lang="perl">my $PI = 2 * atan2 1, 0; |
||
my @nums = map { |
my @nums = map { |
||
1 + 0.5 * sqrt(-2 * log rand) * cos(2 * $PI * rand) |
1 + 0.5 * sqrt(-2 * log rand) * cos(2 * $PI * rand) |
||
} 1..1000;</ |
} 1..1000;</syntaxhighlight> |
||
=={{header| |
=={{header|Phix}}== |
||
{{Trans|Euphoria}} |
|||
{{works with|Rakudo|#22 "Thousand Oaks"}} |
|||
<!--<syntaxhighlight lang="phix">--> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">RandomNormal</span><span style="color: #0000FF;">()</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #7060A8;">log</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">rnd</span><span style="color: #0000FF;">()))</span> <span style="color: #0000FF;">*</span> <span style="color: #7060A8;">cos</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">PI</span><span style="color: #0000FF;">*</span><span style="color: #7060A8;">rnd</span><span style="color: #0000FF;">())</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1000</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">0.5</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">RandomNormal</span><span style="color: #0000FF;">()</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<!--</syntaxhighlight>--> |
|||
=={{header|Phixmonti}}== |
|||
<lang perl6>sub randnorm ($mean, $stddev) { |
|||
<syntaxhighlight lang="phixmonti">include ..\Utilitys.pmt |
|||
$mean + $stddev * sqrt(-2 * log rand) * cos(2 * pi * rand) |
|||
} |
|||
def RandomNormal |
|||
my @nums = randnorm(1, 0.5) xx 1000;</lang> |
|||
drop rand log -2 * sqrt 2 pi * rand * cos * 0.5 * 1 + |
|||
enddef |
|||
1000 var n |
|||
0 n repeat |
|||
getid RandomNormal map |
|||
dup |
|||
sum n / var mean |
|||
"Mean: " print mean print nl |
|||
0 swap n for |
|||
get mean - 2 power rot + swap |
|||
endfor |
|||
swap n / sqrt "Standard deviation: " print print</syntaxhighlight> |
|||
=={{header|PHP}}== |
=={{header|PHP}}== |
||
< |
<syntaxhighlight lang="php">function random() { |
||
return mt_rand() / mt_getrandmax(); |
return mt_rand() / mt_getrandmax(); |
||
} |
} |
||
Line 1,067: | Line 2,152: | ||
$a[$i] = 1.0 + ((sqrt(-2 * log(random())) * cos(2 * $pi * random())) * 0.5); |
$a[$i] = 1.0 + ((sqrt(-2 * log(random())) * cos(2 * $pi * random())) * 0.5); |
||
}</ |
}</syntaxhighlight> |
||
=={{header|Picat}}== |
|||
<syntaxhighlight lang="picat">main => |
|||
_ = random2(), % random seed |
|||
G = [gaussian_dist(1,0.5) : _ in 1..1000], |
|||
println(first_10=G[1..10]), |
|||
println([mean=avg(G),stdev=stdev(G)]), |
|||
nl. |
|||
% Gaussian (Normal) distribution, Box-Muller algorithm |
|||
gaussian01() = Y => |
|||
U = frand(0,1), |
|||
V = frand(0,1), |
|||
Y = sqrt(-2*log(U))*sin(2*math.pi*V). |
|||
gaussian_dist(Mean,Stdev) = Mean + (gaussian01() * Stdev). |
|||
% Variance of Xs |
|||
variance(Xs) = Variance => |
|||
Mu = avg(Xs), |
|||
N = Xs.len, |
|||
Variance = sum([ (X-Mu)**2 : X in Xs ]) / N. |
|||
% Standard deviation |
|||
stdev(Xs) = sqrt(variance(Xs)).</syntaxhighlight> |
|||
{{out}} |
|||
<pre>first_10 = [1.639965415776091,0.705425965005482,0.981532402477848,0.309148743347499,1.252800181962738,0.098829881195179,0.74888084504147,0.181494956495445,1.304931340021904,0.595939453660087] |
|||
[mean = 0.99223677282248,stdev = 0.510336641737154]</pre> |
|||
=={{header|PicoLisp}}== |
=={{header|PicoLisp}}== |
||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang="picolisp">(load "@lib/math.l") |
||
(de randomNormal () # Normal distribution, centered on 0, std dev 1 |
(de randomNormal () # Normal distribution, centered on 0, std dev 1 |
||
Line 1,086: | Line 2,201: | ||
(link (+ 1.0 (/ (randomNormal) 2))) ) ) |
(link (+ 1.0 (/ (randomNormal) 2))) ) ) |
||
(for N (head 7 Result) # Print first 7 results |
(for N (head 7 Result) # Print first 7 results |
||
(prin (format N *Scl) " ") ) )</ |
(prin (format N *Scl) " ") ) )</syntaxhighlight> |
||
{{out}} |
|||
Output: |
|||
<pre>1.500334 1.212931 1.095283 0.433122 0.459116 1.302446 0.402477</pre> |
<pre>1.500334 1.212931 1.095283 0.433122 0.459116 1.302446 0.402477</pre> |
||
=={{header|PL/I}}== |
=={{header|PL/I}}== |
||
<syntaxhighlight lang="pl/i"> |
|||
<lang PL/I> |
|||
/* CONVERTED FROM WIKI FORTRAN */ |
/* CONVERTED FROM WIKI FORTRAN */ |
||
Normal_Random: procedure options (main); |
Normal_Random: procedure options (main); |
||
Line 1,113: | Line 2,228: | ||
put skip edit ( "Standard Deviation = ", sd) (a, F(18,16)); |
put skip edit ( "Standard Deviation = ", sd) (a, F(18,16)); |
||
END Normal_Random; |
END Normal_Random; |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
|||
Output: |
|||
<pre> |
<pre> |
||
Mean = 1.0125630677913652 Standard Deviation = 0.5067289784535284 |
Mean = 1.0125630677913652 Standard Deviation = 0.5067289784535284 |
||
Line 1,124: | Line 2,239: | ||
=={{header|PL/SQL}}== |
=={{header|PL/SQL}}== |
||
<syntaxhighlight lang="pl/sql"> |
|||
<lang PL/SQL>create or replace |
|||
DECLARE |
|||
PROCEDURE PROCEDURE1 AS |
|||
--The desired collection |
|||
TYPE numsColl is TABLE OF NUMBER; |
|||
type t_coll is table of number index by binary_integer; |
|||
nums numsColl; |
|||
l_coll t_coll; |
|||
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; |
|||
c_max pls_integer := 1000; |
|||
BEGIN |
BEGIN |
||
FOR l_counter IN 1 .. c_max |
|||
nums := GenNums(10); |
|||
LOOP |
|||
-- dbms_random.normal delivers normal distributed random numbers with a mean of 0 and a variance of 1 |
|||
DBMS_OUTPUT.PUT_LINE(nums(i)); |
|||
-- We just adjust the values and get the desired result: |
|||
END LOOP; |
|||
l_coll(l_counter) := DBMS_RANDOM.normal * 0.5 + 1; |
|||
END PROCEDURE1;</lang> |
|||
DBMS_OUTPUT.put_line (l_coll(l_counter)); |
|||
END LOOP; |
|||
END; |
|||
</syntaxhighlight> |
|||
=={{header|Pop11}}== |
=={{header|Pop11}}== |
||
< |
<syntaxhighlight lang="pop11">;;; Choose radians as arguments to trigonometic functions |
||
true -> popradians; |
true -> popradians; |
||
Line 1,164: | Line 2,272: | ||
for i from 1 to 1000 do 1.0+0.5*random_normal() endfor; |
for i from 1 to 1000 do 1.0+0.5*random_normal() endfor; |
||
;;; collect them into array |
;;; collect them into array |
||
consvector(1000) -> array;</ |
consvector(1000) -> array;</syntaxhighlight> |
||
=={{header| |
=={{header|PowerShell}}== |
||
Equation adapted from Liberty BASIC |
|||
<lang PureBasic>Procedure.f RandomNormal() |
|||
<syntaxhighlight lang="powershell">function Get-RandomNormal |
|||
; This procedure can return any real number. |
|||
{ |
|||
Protected.f x1, x2 |
|||
[CmdletBinding()] |
|||
Param ( [double]$Mean, [double]$StandardDeviation ) |
|||
$RandomNormal = $Mean + $StandardDeviation * [math]::Sqrt( -2 * [math]::Log( ( Get-Random -Minimum 0.0 -Maximum 1.0 ) ) ) * [math]::Cos( 2 * [math]::PI * ( Get-Random -Minimum 0.0 -Maximum 1.0 ) ) |
|||
return $RandomNormal |
|||
} |
|||
# Standard deviation function for testing |
|||
function Get-StandardDeviation |
|||
{ |
|||
[CmdletBinding()] |
|||
param ( [double[]]$Numbers ) |
|||
$Measure = $Numbers | Measure-Object -Average |
|||
$PopulationDeviation = 0 |
|||
ForEach ($Number in $Numbers) { $PopulationDeviation += [math]::Pow( ( $Number - $Measure.Average ), 2 ) } |
|||
$StandardDeviation = [math]::Sqrt( $PopulationDeviation / ( $Measure.Count - 1 ) ) |
|||
return $StandardDeviation |
|||
} |
|||
# Test |
|||
$RandomNormalNumbers = 1..1000 | ForEach { Get-RandomNormal -Mean 1 -StandardDeviation 0.5 } |
|||
$Measure = $RandomNormalNumbers | Measure-Object -Average |
|||
$Stats = [PSCustomObject]@{ |
|||
Count = $Measure.Count |
|||
Average = $Measure.Average |
|||
StandardDeviation = Get-StandardDeviation -Numbers $RandomNormalNumbers |
|||
} |
|||
$Stats | Format-List |
|||
; random numbers from the open interval ]0, 1[ |
|||
</syntaxhighlight> |
|||
x1 = (Random(999998)+1) / 1000000 ; must be > 0 because of Log(x1) |
|||
{{out}} |
|||
x2 = (Random(999998)+1) / 1000000 |
|||
<pre> |
|||
Count : 1000 |
|||
Average : 1.01206560135809 |
|||
StandardDeviation : 0.489099623426272 |
|||
</pre> |
|||
=={{header|Python}}== |
|||
ProcedureReturn Sqr(-2*Log(x1)) * Cos(2*#PI*x2) |
|||
;Using random.gauss: |
|||
EndProcedure |
|||
<syntaxhighlight lang="python">>>> import random |
|||
>>> values = [random.gauss(1, .5) for i in range(1000)] |
|||
>>> </syntaxhighlight> |
|||
;Quick check of distribution: |
|||
<syntaxhighlight lang="python">>>> def quick_check(numbers): |
|||
count = len(numbers) |
|||
mean = sum(numbers) / count |
|||
sdeviation = (sum((i - mean)**2 for i in numbers) / count)**0.5 |
|||
return mean, sdeviation |
|||
>>> quick_check(values) |
|||
Define i, n=1000 |
|||
(1.0140373306786599, 0.49943411329234066) |
|||
>>> </syntaxhighlight> |
|||
Dim a.q(n-1) |
|||
For i = 0 To n-1 |
|||
a(i) = 1 + 0.5 * RandomNormal() |
|||
Next</lang> |
|||
=={{header|Python}}== |
|||
{{works with|Python|2.5}} |
|||
<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. |
Note that the ''random'' module in the Python standard library supports a number of statistical distribution methods. |
||
;Alternatively using random.normalvariate: |
|||
'''Quick check''' |
|||
< |
<syntaxhighlight lang="python">>>> values = [ random.normalvariate(1, 0.5) for i in range(1000)] |
||
>>> quick_check(values) |
|||
>>> sdeviation = (sum((i - mean)**2 for i in values)/1000)**0.5 |
|||
(0.990099111944864, 0.5029847005836282) |
|||
>>> mean, sdeviation |
|||
>>> </syntaxhighlight> |
|||
(1.0127861555468178, 0.5006682783828207) |
|||
</lang> |
|||
=={{header|R}}== |
=={{header|R}}== |
||
<syntaxhighlight lang="rsplus"># For reproducibility, set the seed: |
|||
<lang r>result <- rnorm(1000, mean=1, sd=0.5)</lang> |
|||
set.seed(12345L) |
|||
result <- rnorm(1000, mean = 1, sd = 0.5)</syntaxhighlight> |
|||
=={{header|Racket}}== |
=={{header|Racket}}== |
||
<syntaxhighlight lang="racket"> |
|||
<lang Racket> |
|||
#lang racket |
#lang racket |
||
(for/list ([i 1000]) |
(for/list ([i 1000]) |
||
(add1 (* (sqrt (* -2 (log (random)))) (cos (* 2 pi (random))) 0.5))) |
(add1 (* (sqrt (* -2 (log (random)))) (cos (* 2 pi (random))) 0.5))) |
||
</syntaxhighlight> |
|||
</lang> |
|||
Alternative: |
|||
<syntaxhighlight lang="racket"> |
|||
#lang racket |
|||
(require math/distributions) |
|||
(sample (normal-dist 1.0 0.5) 1000) |
|||
</syntaxhighlight> |
|||
=={{header|Raku}}== |
|||
(formerly Perl 6) |
|||
{{works with|Rakudo|#22 "Thousand Oaks"}} |
|||
<syntaxhighlight lang="raku" line>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; |
|||
</syntaxhighlight> |
|||
=={{header|Raven}}== |
=={{header|Raven}}== |
||
< |
<syntaxhighlight lang="raven">define PI |
||
-1 acos |
-1 acos |
||
Line 1,224: | Line 2,392: | ||
2 / 1 + |
2 / 1 + |
||
1000 each drop randNormal "%f\n" print</ |
1000 each drop randNormal "%f\n" print</syntaxhighlight> |
||
Quick Check (on linux with code in file rand.rv) |
Quick Check (on linux with code in file rand.rv) |
||
< |
<syntaxhighlight lang="raven">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 |
stdev = 0.497773 |
||
mean = 1.01497</ |
mean = 1.01497</syntaxhighlight> |
||
=={{header|ReScript}}== |
|||
{{trans|OCaml}} |
|||
<syntaxhighlight lang="rescript">let pi = 4.0 *. atan(1.0) |
|||
let random_gaussian = () => { |
|||
1.0 +. |
|||
sqrt(-2.0 *. log(Random.float(1.0))) *. |
|||
cos(2.0 *. pi *. Random.float(1.0)) |
|||
} |
|||
let a = Belt.Array.makeBy(1000, (_) => random_gaussian ()) |
|||
for i in 0 to 10 { |
|||
Js.log(a[i]) |
|||
}</syntaxhighlight> |
|||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
The REXX language doesn't have any "higher math" |
The REXX language doesn't have any "higher math" functions like SQRT/SIN/COS/LN/LOG/EXP/POW/etc., |
||
<br>so we ''hoi polloi'' REXX programmers have to roll our own. |
|||
functions like SIN/COS/LN/LOG/SQRT/POW/etc, |
|||
<br>so we hoi polloi programmers have roll our own. |
|||
<lang rexx>/*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*/ |
|||
Programming note: note the range of the random numbers: (0,1] |
|||
say ' old mean=' mean() |
|||
<br>(that is, random numbers from zero──►unity, excluding zero, including unity). |
|||
say 'old standard deviation=' stddev() |
|||
<syntaxhighlight lang="rexx">/*REXX pgm generates 1,000 normally distributed numbers: mean=1, standard deviation=½.*/ |
|||
numeric digits 20 /*the default decimal digit precision=9*/ |
|||
parse arg n seed . /*allow specification of N and the seed*/ |
|||
if n=='' | n=="," then n=1000 /*N: is the size of the array. */ |
|||
if datatype(seed,'W') then call random ,,seed /*SEED: for repeatable random numbers. */ |
|||
newMean=1 /*the desired new mean (arithmetic avg)*/ |
|||
sd=1/2 /*the desired new standard deviation. */ |
|||
do g=1 for n /*generate N uniform random #'s (0,1].*/ |
|||
#.g = random(1, 1e5) / 1e5 /*REXX's RANDOM BIF generates integers.*/ |
|||
end /*g*/ /* [↑] random integers ──► fractions. */ |
|||
say ' old mean=' mean() |
|||
say 'old standard deviation=' stdDev() |
|||
call pi; pi2=pi * 2 /*define pi and also 2 * pi. */ |
|||
say |
say |
||
do j=1 to n-1 by 2 |
do j=1 to n-1 by 2; m=j+1 /*step through the iterations by two. */ |
||
_=sd * sqrt(ln(#.j) * -2) /*calculate the used-twice expression.*/ |
|||
m=j+1 |
|||
#.j=_ * cos(pi2 * #.m) + newMean /*utilize the Box─Muller method. */ |
|||
#.m= |
#.m=_ * sin(pi2 * #.m) + newMean /*random number must be: (0,1] */ |
||
end /*j*/ |
|||
say ' new mean=' mean() |
|||
say ' |
say 'new standard deviation=' stdDev() |
||
exit /*stick a fork in it, we're all done. */ |
|||
say 'new standard deviation=' stddev() |
|||
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ |
|||
exit /*stick a fork in it, we're done.*/ |
|||
mean: _=0; do k=1 for n; _=_ + #.k; end; return _/n |
|||
/*──────────────────────────────────subroutines─────────────────────────*/ |
|||
stdDev: _avg=mean(); _=0; do k=1 for n; _=_ + (#.k - _avg)**2; end; return sqrt(_/n) |
|||
e: e =2.7182818284590452353602874713526624977572470936999595749669676277240766303535; return e /*digs overkill*/ |
|||
stddev: _avg=mean(); _=0; do k=1 for n; _=_+(#.k-_avg)**2; end; return sqrt(_/n) |
|||
pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862; return pi /* " " */ |
|||
e: e=2.7182818284590452353602874713526624977572470936999595749669676277240766303535; return e |
|||
r2r: return arg(1) // (pi() * 2) /*normalize ang*/ |
|||
pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862; return pi |
|||
sin: procedure; parse arg x;x=r2r(x);numeric fuzz min(5,digits()-3);if abs(x)=pi then return 0;return .sincos(x,x,1) |
|||
r2r: return arg(1)//(2*pi()) |
|||
.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 |
|||
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ |
|||
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 |
|||
ln: procedure; parse arg x,f; call e; ig= x>1.5; is=1 - 2 * (ig\==1); ii=0; xx=x |
|||
g=.5*(g+x/g); end; numeric digits d; return g/1 |
|||
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 |
|||
.sqrtGuess: numeric form; m.=11; p=d+d%4+2 |
|||
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 |
|||
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 |
|||
cos: procedure; parse arg x; x=r2r(x); a=abs(x); hpi=pi * .5 |
|||
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) |
|||
numeric fuzz min(6, digits() - 3); if a=pi then return -1 |
|||
if a=hpi | a=hpi*3 then return 0; if a=pi/3 then return .5 |
|||
if a=pi * 2/3 then return -.5; return .sinCos(1,1,-1) |
|||
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ |
|||
.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 |
|||
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits; h=d+6 |
|||
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</lang> |
|||
numeric form; parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g * .5'e'_ %2 |
|||
'''output''' when using the default input |
|||
m.=9; do j=0 while h>9; m.j=h; h=h%2 + 1; end /*j*/ |
|||
<pre style="overflow:scroll"> |
|||
do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/ |
|||
numeric digits d; return g/1</syntaxhighlight> |
|||
old standard deviation= 0.292088074 |
|||
'''output''' when using the default inputs: |
|||
<pre> |
|||
old mean= 0.5015724 |
|||
old standard deviation= 0.28652466389342471402 |
|||
new mean= 0. |
new mean= 0.98807025356443262689 |
||
new standard deviation= 0. |
new standard deviation= 0.50002924192766720838 |
||
</pre> |
|||
=={{header|Ring}}== |
|||
<syntaxhighlight lang="ring"> |
|||
for i = 1 to 10 |
|||
see random(i) + nl |
|||
next i |
|||
</syntaxhighlight> |
|||
=={{header|RPL}}== |
|||
≪ RAND LN NEG 2 * √ |
|||
RAND 2 * π * COS * |
|||
→NUM 2 / 1 + |
|||
≫ '<span style="color:blue>RANDN</span>' STO |
|||
≪ CL∑ |
|||
1 1000 '''START''' <span style="color:blue>RANDN</span> ∑+ '''NEXT''' |
|||
MEAN PSDEV |
|||
≫ '<span style="color:blue>TASK</span>' STO |
|||
{{out}} |
|||
<pre> |
|||
1: .990779804949 |
|||
2: .487204045227 |
|||
</pre> |
</pre> |
||
The collection is stored in a predefined array named <code>∑DAT</code>, which is automatically created/updated when using the <code>∑+</code> instruction and remains available until the user decides to purge it, typically by calling the <code>CL∑</code> command. |
|||
=={{header|Ruby}}== |
=={{header|Ruby}}== |
||
< |
<syntaxhighlight lang="ruby">Array.new(1000) { 1 + Math.sqrt(-2 * Math.log(rand)) * Math.cos(2 * Math::PI * rand) }</syntaxhighlight> |
||
=={{header| |
=={{header|Rust}}== |
||
{{libheader|rand}} |
|||
<lang runbasic>dim a(1000) |
|||
'''Using a for-loop:''' |
|||
pi = 22/7 |
|||
<syntaxhighlight lang="rust">extern crate rand; |
|||
for i = 1 to 1000 |
|||
use rand::distributions::{Normal, IndependentSample}; |
|||
a( i) = 1 + .5 * (sqr(-2 * log(rnd(0))) * cos(2 * pi * rnd(0))) |
|||
next i</lang> |
|||
fn main() { |
|||
let mut rands = [0.0; 1000]; |
|||
let normal = Normal::new(1.0, 0.5); |
|||
let mut rng = rand::thread_rng(); |
|||
for num in rands.iter_mut() { |
|||
*num = normal.ind_sample(&mut rng); |
|||
} |
|||
}</syntaxhighlight> |
|||
'''Using iterators:''' |
|||
<syntaxhighlight lang="rust">extern crate rand; |
|||
use rand::distributions::{Normal, IndependentSample}; |
|||
fn main() { |
|||
let rands: Vec<_> = { |
|||
let normal = Normal::new(1.0, 0.5); |
|||
let mut rng = rand::thread_rng(); |
|||
(0..1000).map(|_| normal.ind_sample(&mut rng)).collect() |
|||
}; |
|||
}</syntaxhighlight> |
|||
=={{header|SAS}}== |
=={{header|SAS}}== |
||
<syntaxhighlight lang="sas"> |
|||
<lang SAS> |
|||
/* Generate 1000 random numbers with mean 1 and standard deviation 0.5. |
/* Generate 1000 random numbers with mean 1 and standard deviation 0.5. |
||
SAS version 9.2 was used to create this code.*/ |
SAS version 9.2 was used to create this code.*/ |
||
Line 1,308: | Line 2,542: | ||
end; |
end; |
||
run; |
run; |
||
</syntaxhighlight> |
|||
</lang> |
|||
Results: |
Results: |
||
<pre> |
<pre> |
||
Line 1,322: | Line 2,556: | ||
=={{header|Sather}}== |
=={{header|Sather}}== |
||
< |
<syntaxhighlight lang="sather">class MAIN is |
||
main is |
main is |
||
a:ARRAY{FLTD} := #(1000); |
a:ARRAY{FLTD} := #(1000); |
||
Line 1,340: | Line 2,574: | ||
#OUT + "dev " + dev + "\n"; |
#OUT + "dev " + dev + "\n"; |
||
end; |
end; |
||
end;</ |
end;</syntaxhighlight> |
||
=={{header|Scala}}== |
=={{header|Scala}}== |
||
===One liner=== |
|||
<lang scala>List.fill(1000)(1.0 + 0.5 * scala.util.Random.nextGaussian)</lang> |
|||
<syntaxhighlight lang="scala">List.fill(1000)(1.0 + 0.5 * scala.util.Random.nextGaussian)</syntaxhighlight> |
|||
===Academic=== |
|||
<syntaxhighlight lang="scala"> |
|||
object RandomNumbers extends App { |
|||
val distribution: LazyList[Double] = { |
|||
def randomNormal: Double = 1.0 + 0.5 * scala.util.Random.nextGaussian |
|||
def normalDistribution(a: Double): LazyList[Double] = a #:: normalDistribution(randomNormal) |
|||
normalDistribution(randomNormal) |
|||
} |
|||
/* |
|||
* Let's test it |
|||
*/ |
|||
def calcAvgAndStddev[T](ts: Iterable[T])(implicit num: Fractional[T]): (T, Double) = { |
|||
val mean: T = |
|||
num.div(ts.sum, num.fromInt(ts.size)) // Leaving with type of function T |
|||
// Root of mean diffs |
|||
val stdDev = Math.sqrt(ts.map { x => |
|||
val diff = num.toDouble(num.minus(x, mean)) |
|||
diff * diff |
|||
}.sum / ts.size) |
|||
(mean, stdDev) |
|||
} |
|||
println(calcAvgAndStddev(distribution.take(1000))) // e.g. (1.0061433267806525,0.5291834867560893) |
|||
} |
|||
</syntaxhighlight> |
|||
=={{header|Scheme}}== |
=={{header|Scheme}}== |
||
< |
<syntaxhighlight lang="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)) |
(define ((c-rand seed)) (set! seed (remainder (+ (* 1103515245 seed) 12345) 2147483648)) (quotient seed 65536)) |
||
Line 1,392: | Line 2,658: | ||
(mean-sdev v) |
(mean-sdev v) |
||
; (0.9562156817697293 0.5097087109575911)</ |
; (0.9562156817697293 0.5097087109575911)</syntaxhighlight> |
||
=={{header|Seed7}}== |
=={{header|Seed7}}== |
||
< |
<syntaxhighlight lang="seed7">$ include "seed7_05.s7i"; |
||
include "float.s7i"; |
include "float.s7i"; |
||
include "math.s7i"; |
include "math.s7i"; |
||
Line 1,419: | Line 2,685: | ||
rands[i] := 1.0 + 0.5 * randomNormal; |
rands[i] := 1.0 + 0.5 * randomNormal; |
||
end for; |
end for; |
||
end func;</ |
end func;</syntaxhighlight> |
||
=={{header|Sidef}}== |
|||
<syntaxhighlight lang="ruby">var arr = 1000.of { 1 + (0.5 * sqrt(-2 * 1.rand.log) * cos(Num.tau * 1.rand)) } |
|||
arr.each { .say }</syntaxhighlight> |
|||
=={{header|Standard ML}}== |
=={{header|Standard ML}}== |
||
Line 1,425: | Line 2,695: | ||
SML/NJ has two structures for random numbers: |
SML/NJ has two structures for random numbers: |
||
1) Rand (a linear congruential generator). |
|||
1) Rand (a linear congruential generator). You create the generator by calling <code>Rand.mkRandom</code> with a seed (of <code>word</code> type). 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. |
|||
You create the generator by calling <code>Rand.mkRandom</code> with a seed (of <code>word</code> type). |
|||
<lang sml>val seed = 0w42; |
|||
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. |
|||
<syntaxhighlight lang="sml">val seed = 0w42; |
|||
val gen = Rand.mkRandom seed; |
val gen = Rand.mkRandom seed; |
||
fun random_gaussian () = |
fun random_gaussian () = |
||
1.0 + Math.sqrt (~2.0 * Math.ln (Rand.norm (gen ()))) * Math.cos (2.0 * Math.pi * Rand.norm (gen ())); |
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 ());</ |
val a = List.tabulate (1000, fn _ => random_gaussian ());</syntaxhighlight> |
||
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. |
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. |
||
< |
<syntaxhighlight lang="sml">val seed = (47,42); |
||
val gen = Random.rand seed; |
val gen = Random.rand seed; |
||
fun random_gaussian () = |
fun random_gaussian () = |
||
1.0 + Math.sqrt (~2.0 * Math.ln (Random.randReal gen)) * Math.cos (2.0 * Math.pi * Random.randReal gen); |
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 ());</ |
val a = List.tabulate (1000, fn _ => random_gaussian ());</syntaxhighlight> |
||
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. |
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. |
||
{{works with|Poly/ML}} |
|||
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: |
|||
<syntaxhighlight lang="sml"> |
|||
val urandomlist = fn seed => fn n => |
|||
let |
|||
val uniformdeviate = fn seed => |
|||
let |
|||
val in31m = (Real.fromInt o Int32.toInt ) (getOpt (Int32.maxInt,0) ); |
|||
val in31 = in31m +1.0; |
|||
val s1 = 41160.0; |
|||
val s2 = 950665216.0; |
|||
val v = Real.realFloor seed; |
|||
val val1 = v*s1; |
|||
val val2 = v*s2; |
|||
val next1 = Real.fromLargeInt (Real.toLargeInt IEEEReal.TO_NEGINF (val1/in31)) ; |
|||
val next2 = Real.rem(Real.realFloor(val2/in31) , in31m ); |
|||
val valt = val1+val2 - (next1+next2)*in31m; |
|||
val nextt = Real.realFloor(valt/in31m); |
|||
val valt = valt - nextt*in31m; |
|||
in |
|||
(valt/in31m,valt) |
|||
end; |
|||
val store = ref (0.0,0.0); |
|||
val rec u = fn S => fn 0 => [] | n=> (store:=uniformdeviate S; (#1 (!store)):: (u (#2 (!store)) (n-1))) ; |
|||
in |
|||
u seed n |
|||
end; |
|||
local |
|||
open Math |
|||
in |
|||
val bmconv = fn urand => fn vrand => 1.0+0.5*(sqrt(~2.0*ln urand)*cos (2.0*pi*vrand) ) |
|||
end; |
|||
val rec makeNormals = fn once => fn u::v::[] => [once u v] | |
|||
u::v::rm => (once u v )::(makeNormals once rm ); |
|||
val anyrealseed=1009.0 ; |
|||
makeNormals bmconv (urandomlist anyrealseed 2000); |
|||
</syntaxhighlight> |
|||
=={{header|Stata}}== |
|||
<syntaxhighlight lang="stata">clear all |
|||
set obs 1000 |
|||
gen x=rnormal(1,0.5)</syntaxhighlight> |
|||
=== Mata === |
|||
<syntaxhighlight lang="stata">a = rnormal(1000,1,1,0.5)</syntaxhighlight> |
|||
=={{header|Tcl}}== |
=={{header|Tcl}}== |
||
< |
<syntaxhighlight lang="tcl">package require Tcl 8.5 |
||
variable ::pi [expr acos(0)] |
variable ::pi [expr acos(0)] |
||
proc ::tcl::mathfunc::nrand {} { |
proc ::tcl::mathfunc::nrand {} { |
||
Line 1,452: | Line 2,773: | ||
for {set i 0} {$i < 1000} {incr i} { |
for {set i 0} {$i < 1000} {incr i} { |
||
lappend result [expr {$mean + $stddev*nrand()}] |
lappend result [expr {$mean + $stddev*nrand()}] |
||
}</ |
}</syntaxhighlight> |
||
=={{header| |
=={{header|TorqueScript}}== |
||
<syntaxhighlight lang="tqs">for (%i = 0; %i < 1000; %i++) |
|||
Calculator symbol translations: |
|||
%list[%i] = 1 + mSqrt(-2 * mLog(getRandom())) * mCos(2 * $pi * getRandom());</syntaxhighlight> |
|||
"STO" arrow: → |
|||
Square root sign: √ |
|||
ClrList L<sub>1</sub> |
|||
Radian |
|||
For(A,1,1000) |
|||
√(-2*ln(rand))*cos(2*π*A)→L<sub>1</sub>(A) |
|||
End |
|||
=={{header|Ursala}}== |
=={{header|Ursala}}== |
||
Line 1,477: | Line 2,789: | ||
a standard normal distribution. Mean and standard deviation |
a standard normal distribution. Mean and standard deviation |
||
library functions are also used in this example. |
library functions are also used in this example. |
||
< |
<syntaxhighlight lang="ursala">#import nat |
||
#import flo |
#import flo |
||
Line 1,490: | Line 2,802: | ||
^(mean,stdev)* < |
^(mean,stdev)* < |
||
pop_stats(1.,0.5) 1000, |
pop_stats(1.,0.5) 1000, |
||
sample_stats(1.,0.5) 1000></ |
sample_stats(1.,0.5) 1000></syntaxhighlight> |
||
The output shows the mean and standard deviation for both sample vectors, |
The output shows the mean and standard deviation for both sample vectors, |
||
the latter being exact by construction. |
the latter being exact by construction. |
||
Line 1,496: | Line 2,808: | ||
(1.004504e+00,4.915525e-01), |
(1.004504e+00,4.915525e-01), |
||
(1.000000e+00,5.000000e-01)></pre> |
(1.000000e+00,5.000000e-01)></pre> |
||
=={{header|Visual FoxPro}}== |
|||
<syntaxhighlight lang="vfp"> |
|||
LOCAL i As Integer, m As Double, n As Integer, sd As Double |
|||
py = PI() |
|||
SET TALK OFF |
|||
SET DECIMALS TO 6 |
|||
CREATE CURSOR gdev (deviate B(6)) |
|||
RAND(-1) |
|||
n = 1000 |
|||
m = 1 |
|||
sd = 0.5 |
|||
CLEAR |
|||
FOR i = 1 TO n |
|||
INSERT INTO gdev VALUES (GaussDev(m, 1/sd)) |
|||
ENDFOR |
|||
CALCULATE AVG(deviate), STD(deviate) TO m, sd |
|||
? "Mean", m, "Std Dev", sd |
|||
SET TALK ON |
|||
SET DECIMALS TO |
|||
FUNCTION GaussDev(mean As Double, sdev As Double) As Double |
|||
LOCAL z As Double |
|||
z = SQRT(-2*LOG(RAND()))*COS(py*RAND()) |
|||
IF sdev # 0 |
|||
z = mean + z/sdev |
|||
ENDIF |
|||
RETURN z |
|||
ENDFUNC |
|||
</syntaxhighlight> |
|||
=={{header|V (Vlang)}}== |
|||
<syntaxhighlight lang="Vlang"> |
|||
import crypto.rand |
|||
fn main() { |
|||
mut nums := []u64{} |
|||
for _ in 0..1000 { |
|||
nums << rand.int_u64(10000) or {0} // returns random unsigned 64-bit integer from real OS source of entropy |
|||
} |
|||
println(nums) |
|||
} |
|||
</syntaxhighlight> |
|||
=={{header|Wren}}== |
|||
<syntaxhighlight lang="wren">import "random" for Random |
|||
var rand = Random.new() |
|||
var randNormal = Fn.new { (-2 * rand.float().log).sqrt * (2 * Num.pi * rand.float()).cos } |
|||
var stdDev = Fn.new { |a, m| |
|||
var c = a.count |
|||
return ((a.reduce(0) { |acc, x| acc + x*x } - m*m*c) / c).sqrt |
|||
} |
|||
var n = 1000 |
|||
var numbers = List.filled(n, 0) |
|||
var mu = 1 |
|||
var sigma = 0.5 |
|||
var sum = 0 |
|||
for (i in 0...n) { |
|||
numbers[i] = mu + sigma*randNormal.call() |
|||
sum = sum + numbers[i] |
|||
} |
|||
var mean = sum / n |
|||
System.print("Actual mean : %(mean)") |
|||
System.print("Actual std dev: %(stdDev.call(numbers, mean))")</syntaxhighlight> |
|||
{{out}} |
|||
Sample run: |
|||
<pre> |
|||
Actual mean : 1.0053988699746 |
|||
Actual std dev: 0.4961645117026 |
|||
</pre> |
|||
=={{header|XPL0}}== |
|||
{{trans|C}} |
|||
<syntaxhighlight lang "XPL0">define PI = 3.14159265358979323846; |
|||
func real DRand; \Uniform distribution, [0..1] |
|||
return float(Ran(1_000_000)) / 1e6; |
|||
func real RandomNormal; \Normal distribution, centered on 0, std dev 1 |
|||
return sqrt(-2.*Log(DRand)) * Cos(2.*PI*DRand); |
|||
int I; |
|||
real Rands(1000); |
|||
for I:= 0 to 1000-1 do |
|||
Rands(I):= 1.0 + 0.5*RandomNormal</syntaxhighlight> |
|||
=={{header|Yorick}}== |
=={{header|Yorick}}== |
||
Returns array of ''count'' random numbers with mean 0 and standard deviation 1. |
Returns array of ''count'' random numbers with mean 0 and standard deviation 1. |
||
< |
<syntaxhighlight lang="yorick">func random_normal(count) { |
||
return sqrt(-2*log(random(count))) * cos(2*pi*random(count)); |
return sqrt(-2*log(random(count))) * cos(2*pi*random(count)); |
||
}</ |
}</syntaxhighlight> |
||
Example of basic use: |
Example of basic use: |
||
Line 1,517: | Line 2,919: | ||
0.496853</pre> |
0.496853</pre> |
||
=={{header| |
=={{header|zkl}}== |
||
<syntaxhighlight lang="zkl">fcn mkRand(mean,sd){ //normally distributed random w/mean & standard deviation |
|||
pi:=(0.0).pi; // using the Box–Muller transform |
|||
Here we have converted the QBasic code to suit the ZX Spectrum: |
|||
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 }) |
|||
<lang zxbasic>10 RANDOMIZE 0 : REM seeds random number generator based on uptime |
|||
}</syntaxhighlight> |
|||
20 DIM a(1000) |
|||
This creates a new random number generator, now to use it: |
|||
30 CLS |
|||
<syntaxhighlight lang="zkl">var g=mkRand(1,0.5); |
|||
40 FOR i = 1 TO 1000 |
|||
ns:=(0).pump(1000,List,g); // 1000 rands with mean==1 & sd==1/2 |
|||
50 LET a(i) = 1 + SQR(-2 * LN(RND)) * COS(2 * PI * RND) |
|||
mean:=(ns.sum(0.0)/1000); //-->1.00379 |
|||
60 NEXT i</lang> |
|||
// calc sd of list of numbers: |
|||
(ns.reduce('wrap(p,n){p+(n-mean).pow(2)},0.0)/1000).sqrt() //-->0.494844</syntaxhighlight> |
|||
{{omit from|GUISS}} |
|||
{{omit from|UNIX Shell|From the shell, we simply invoke the awk solution}} |
|||
[[Category:Randomness]] |
Latest revision as of 00:23, 1 May 2024
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
You are encouraged to solve this task according to the task description, using any language you may know.
- Task
Generate a collection filled with 1000 normally distributed random (or pseudo-random) numbers with a mean of 1.0 and a standard deviation of 0.5
Many libraries only generate uniformly distributed random numbers. If so, you may use one of these algorithms.
- Related task
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;
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:
( 0.693461, 0.948424, 0.482261, 1.045939, 0.890818, 1.467935, 0.604153, 0.804811, 0.690227, 0.83462 ... )
Arturo
rnd: function []-> (random 0 10000)//10000
rands: map 1..1000 'x [
1 + (sqrt neg 2 * ln rnd) * (cos 2 * pi * rnd)
]
print rands
- Output:
0.6219537961087694 1.279396486161406 1.619019280815647 2.119538294228789 -0.1598383851981044 2.67797211803156 0.9304703037226587 1.629254364659528 -0.4171704717398712 0.9082342931486092 -0.5929704390625219 2.117000897984871 -0.1981633787460266 0.01132471973856153 2.102359263212924 0.2408823232884222 2.046195035792376 0.6374831627030295 0.000839808324124558 1.117061838266626 0.7413355299469649 0.4485598815755762 2.999434800016997 2.560580541932842 1.703197984879731 2.889159248353575 -1.800157205708138 1.756810020187321 0.7136708180852145 0.5929151678321705 0.332519993787973 2.660212054362758 0.5835660585480075 0.8527946892567934 1.640573993747053 0.09471843345263908 1.051402997891346 1.116149156137905 -0.7400139019343499 1.782572831979232 2.531779039786426 0.5240268064639871 0.07099232630526586 -0.854892656700071 1.54381929430469 -0.4416899008614745 0.4274356035015117 0.7350027625573482 2.153583935076981 1.461215281535983 -1.041723064151266 2.338060763553139 -0.1492967916030414 0.3799517724040202 0.4577924541353815 0.673317567666373 -2.27731583876462 1.28480355806061 -0.6925811023772748 -0.2642224122781984 0.6590513830891744 2.55537133425143 1.67933335469247 0.8659013395355968 -1.211026941441126 0.9524579534222226 -0.1931750631835656 -0.5119479869237693 0.1814749003063878 3.03139579963414 ...
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
}
Avail
Method "U(_,_)" is
[
lower : number,
upper : number
|
divisor ::= ((1<<32)) ÷ (upper - lower)→double;
map a pRNG through [i : integer | (i ÷ divisor) + lower]
];
Method "a Marsaglia polar sampler" is
[
generator for
[
yield : [double]→⊤
|
source ::= U(-1, 1);
Repeat [
x ::= take 1 from source[1];
y ::= take 1 from source[1];
s ::= x^2 + y^2;
If 0 < s < 1 then
[
factor ::= ((-2 × ln s) ÷ s) ^ 0.5;
yield(x × factor);
yield(y × factor);
];
]
]
];
// 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];
values ::= take 1000 from sampler;
AWK
One-liner:
$ 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}'
Readable version:
function r() {
return sqrt( -2*log( rand() ) ) * cos(6.2831853*rand() )
}
BEGIN {
n=1000
for(i=0;i<n;i++) {
x = 1 + 0.5*r()
s = s" "x
}
print s
}
- Output:
first few values only
0.783753 1.16682 1.17989 1.14975 1.34784 0.29296 0.979227 1.04402 0.567835 1.58812 0.465559 1.27186 0.324533 0.725827 -0.0626549 0.632273 1.0145 1.3387 0.861667 1.04147 1.2576 1.02497 0.58453 0.9619 1.26902 0.851048 -0.126259 0.863256
...
BASIC
ANSI BASIC
100 REM Random numbers
110 RANDOMIZE
120 DEF RandomNormal = COS(2 * PI * RND) * SQR(-2 * LOG(RND))
130 DIM R(0 TO 999)
140 LET Sum = 0
150 FOR I = 0 TO 999
160 LET R(I) = 1 + RandomNormal / 2
170 LET Sum = Sum + R(I)
180 NEXT I
190 LET Mean = Sum / 1000
200 LET Sum = 0
210 FOR I = 0 TO 999
220 LET Sum = Sum + (R(I) - Mean) ^ 2
230 NEXT I
240 LET SD = SQR(Sum / 1000)
250 PRINT "Mean is "; Mean
260 PRINT "Standard Deviation is"; SD
270 PRINT
280 END
- Output:
Two runs.
Mean is 1.00216454061435 Standard Deviation is .504515904812839
Mean is .995781408878628 Standard Deviation is .499307289407576
Applesoft BASIC
The Commodore BASIC code works in Applesoft BASIC.
BASIC256
# Generates normally distributed random numbers with mean 0 and standard deviation 1
function randomNormal()
return cos(2.0 * pi * rand) * sqr(-2.0 * log(rand))
end function
dim r(1000)
sum = 0.0
# Generate 1000 normally distributed random numbers
# with mean 1 and standard deviation 0.5
# and calculate their sum
for i = 0 to 999
r[i] = 1.0 + randomNormal() / 2.0
sum += r[i]
next i
mean = sum / 1000.0
sum = 0.0
# Now calculate their standard deviation
for i = 0 to 999
sum += (r[i] - mean) ^ 2.0
next i
sd = sqr(sum/1000.0)
print "Mean is "; mean
print "Standard Deviation is "; sd
end
- Output:
Mean is 1.002092 Standard Deviation is 0.4838570687
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
Chipmunk Basic
10 ' Random numbers
20 randomize timer
30 dim r(999)
40 sum = 0
50 for i = 0 to 999
60 r(i) = 1+randomnormal()/2
70 sum = sum+r(i)
80 next
90 mean = sum/1000
100 sum = 0
110 for i = 0 to 999
120 sum = sum+(r(i)-mean)^2
130 next
140 sd = sqr(sum/1000)
150 print "Mean is ";mean
160 print "Standard Deviation is ";sd
170 print
180 end
500 sub randomnormal()
510 randomnormal = cos(2*pi*rnd(1))*sqr(-2*log(rnd(1)))
520 end sub
- Output:
Two runs.
Mean is 1.007087 Standard Deviation is 0.496848
Mean is 0.9781 Standard Deviation is 0.508147
Commodore BASIC
10 DIM AR(999): DIM DE(999)
20 FOR N = 0 TO 999
30 AR(N)= 0 + SQR(-1.3*LOG(RND(1))) * COS(1.2*PI*RND(1))
40 NEXT N
50 :
60 REM SUM
70 LET SU = 0
80 FOR N = 0 TO 999
90 LET SU = SU + AR(N)
100 NEXT N
110 :
120 REM MEAN
130 LET ME= 0
140 LET ME = SU/1000
150 :
160 REM DEVIATION
170 FOR N = 0 TO 999
180 T = AR(N)-ME: REM SUBTRACT MEAN FROM NUMBER
190 T = T * T: REM SQUARE THE RESULT
200 DE(N) = T : REM STORE IN DEVIATION ARRAY
210 NEXT N
220 LET DS=0: REM SUM OF DEVIATION ARRAY
230 FOR N = 0 TO 999
240 LET DS = DS + DE(N)
250 NEXT N
260 LET DM=0: REM MEAN OF DEVIATION ARRAY
270 LET DM = DS / 1000
280 LET DE = 0:
290 LET DE = SQR(DM)
300 :
310 PRINT "MEAN = "ME
320 PRINT "STANDARD DEVIATION ="DE
330 END
FreeBASIC
' FB 1.05.0 Win64
Const pi As Double = 3.141592653589793
Randomize
' Generates normally distributed random numbers with mean 0 and standard deviation 1
Function randomNormal() As Double
Return Cos(2.0 * pi * Rnd) * Sqr(-2.0 * Log(Rnd))
End Function
Dim r(0 To 999) As Double
Dim sum As Double = 0.0
' Generate 1000 normally distributed random numbers
' with mean 1 and standard deviation 0.5
' and calculate their sum
For i As Integer = 0 To 999
r(i) = 1.0 + randomNormal/2.0
sum += r(i)
Next
Dim mean As Double = sum / 1000.0
Dim sd As Double
sum = 0.0
' Now calculate their standard deviation
For i As Integer = 0 To 999
sum += (r(i) - mean) ^ 2.0
Next
sd = Sqr(sum/1000.0)
Print "Mean is "; mean
Print "Standard Deviation is"; sd
Print
Print "Press any key to quit"
Sleep
Sample result:
- Output:
Mean is 1.000763573902885 Standard Deviation is 0.500653063426955
FutureBasic
Note: To generate the random number, rather than using FB's native "rnd" function, this code wraps C code into the RandomZeroToOne function.
window 1
local fn RandomZeroToOne as double
double result
cln result = (double)( (rand() % 100000 ) * 0.00001 );
end fn = result
local fn RandomGaussian as double
double r = fn RandomZeroToOne
end fn = 1 + .5 * ( sqr( -2 * log(r) ) * cos( 2 * pi * r ) )
long i
double mean, std, a(1000)
for i = 1 to 1000
a(i) = fn RandomGaussian
mean += a(i)
next
mean = mean / 1000
for i = 1 to 1000
std += ( a(i) - mean )^2
next
std = std / 1000
print " Average: "; mean
print "Standard Deviation: "; std
HandleEvents
- Output:
Average: 1.053724951604593 Standard Deviation: 0.2897370762627166
GW-BASIC
The Commodore BASIC code works in GW-BASIC.
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
Minimal BASIC
10 REM Random numbers
20 LET P = 4*ATN(1)
30 RANDOMIZE
40 DEF FNN = COS(2*P*RND)*SQR(-2*LOG(RND))
50 DIM R(999)
60 LET S = 0
70 FOR I = 0 TO 999
80 LET R(I) = 1+FNN/2
90 LET S = S+R(I)
100 NEXT I
110 LET M = S/1000
120 LET S = 0
130 FOR I = 0 TO 999
140 LET S = S+(R(I)-M)^2
150 NEXT I
160 LET D = SQR(S/1000)
170 PRINT "Mean is "; M
180 PRINT "Standard Deviation is"; D
190 PRINT
200 END
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
QuickBASIC
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
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
TI-83 BASIC
Built-in function: randNorm()
randNorm(1,.5)
Or by a program:
Calculator symbol translations:
"STO" arrow: →
Square root sign: √
ClrList L1 Radian For(A,1,1000) √(-2*ln(rand))*cos(2*π*A)→L1(A) End
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
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;
}
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
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<double> 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;
}
Clojure
(import '(java.util Random))
(def normals
(let [r (Random.)]
(take 1000 (repeatedly #(-> r .nextGaussian (* 0.5) (+ 1.0))))))
COBOL
IDENTIFICATION DIVISION.
PROGRAM-ID. RANDOM.
AUTHOR. Bill Gunshannon
INSTALLATION. Home.
DATE-WRITTEN. 14 January 2022.
************************************************************
** Program Abstract:
** Able to get the Mean to be really close to 1.0 but
** couldn't get the Standard Deviation any closer than
** .3 to .4.
************************************************************
DATA DIVISION.
WORKING-STORAGE SECTION.
01 Sample-Size PIC 9(5) VALUE 1000.
01 Total PIC 9(10)V9(5) VALUE 0.0.
01 Arith-Mean PIC 999V999 VALUE 0.0.
01 Std-Dev PIC 999V999 VALUE 0.0.
01 Seed PIC 999V999.
01 TI PIC 9(8).
01 Idx PIC 99999 VALUE 0.
01 Intermediate PIC 9(10)V9(5) VALUE 0.0.
01 Rnd-Work.
05 Rnd-Tbl
OCCURS 1 TO 99999 TIMES DEPENDING ON Sample-Size.
10 Rnd PIC 9V9999999 VALUE 0.0.
PROCEDURE DIVISION.
Main-Program.
ACCEPT TI FROM TIME.
MOVE FUNCTION RANDOM(TI) TO Seed.
PERFORM WITH TEST AFTER VARYING Idx
FROM 1 BY 1
UNTIL Idx = Sample-Size
COMPUTE Intermediate =
(FUNCTION RANDOM() * 2.01)
MOVE Intermediate TO Rnd(Idx)
END-PERFORM.
PERFORM WITH TEST AFTER VARYING Idx
FROM 1 BY 1
UNTIL Idx = Sample-Size
COMPUTE Total = Total + Rnd(Idx)
END-PERFORM.
COMPUTE Arith-Mean = Total / Sample-Size.
DISPLAY "Mean: " Arith-Mean.
PERFORM WITH TEST AFTER VARYING Idx
FROM 1 BY 1
UNTIL Idx = Sample-Size
COMPUTE Intermediate =
Intermediate + (Rnd(Idx) - Arith-Mean) ** 2
END-PERFORM.
COMPUTE Std-Dev = Intermediate / Sample-Size.
DISPLAY "Std-Dev: " Std-Dev.
STOP RUN.
END PROGRAM RANDOM.
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)))
Crystal
n, mean, sd, tau = 1000, 1, 0.5, (2 * Math::PI)
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 )
puts "mean = #{mean}, standard deviation = #{standev}"
- Output:
mean = 1.0093442539237896, standard deviation = 0.504694489463623
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();
}
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;
}
}
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
DWScript
var values : array [0..999] of Float;
var i : Integer;
for i := values.Low to values.High do
values[i] := RandG(1, 0.5);
E
accum [] for _ in 1..1000 { _.with(entropy.nextGaussian()) }
EasyLang
numfmt 5 0
e = 2.7182818284590452354
for i = 1 to 1000
a[] &= 1 + 0.5 * sqrt (-2 * log10 randomf / log10 e) * cos (360 * randomf)
.
for v in a[]
avg += v / len a[]
.
print "Average: " & avg
for v in a[]
s += pow (v - avg) 2
.
s = sqrt (s / len a[])
print "Std deviation: " & s
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
Elena
ELENA 6.x :
import extensions;
import extensions'math;
randomNormal()
{
^ cos(2 * Pi_value * randomGenerator.nextReal())
* sqrt(-2 * ln(randomGenerator.nextReal()))
}
public program()
{
real[] a := new real[](1000);
real tAvg := 0;
for (int x := 0; x < a.Length; x += 1)
{
a[x] := (randomNormal()) / 2 + 1;
tAvg += a[x]
};
tAvg /= a.Length;
console.printLine("Average: ", tAvg);
real s := 0;
for (int x := 0; x < a.Length; x += 1)
{
s += power(a[x] - tAvg, 2)
};
s := sqrt(s / 1000);
console.printLine("Standard Deviation: ", s);
console.readChar()
}
- Output:
Average: 0.9842420481571 Standard Deviation: 0.5109070975558
Elixir
defmodule Random do
def normal(mean, sd) do
{a, b} = {:rand.uniform, :rand.uniform}
mean + sd * (:math.sqrt(-2 * :math.log(a)) * :math.cos(2 * :math.pi * b))
end
end
std_dev = fn (list) ->
mean = Enum.sum(list) / length(list)
sd = Enum.reduce(list, 0, fn x,acc -> acc + (x-mean)*(x-mean) end) / length(list)
|> :math.sqrt
IO.puts "Mean: #{mean},\tStdDev: #{sd}"
end
xs = for _ <- 1..1000, do: Random.normal(1.0, 0.5)
std_dev.(xs)
- Output:
Mean: 1.009079383094275, StdDev: 0.4991894476975088
used Erlang function :rand.normal
xs = for _ <- 1..1000, do: 1.0 + :rand.normal * 0.5
std_dev.(xs)
- Output:
Mean: 0.9955701150615597, StdDev: 0.5036412260426065
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
ERRE
PROGRAM DISTRIBUTION
!
! for rosettacode.org
!
! formulas taken from TI-59 Master Library manual
CONST NUM_ITEM=1000
!VAR SUMX#,SUMX2#,R1#,R2#,Z#,I%
DIM A#[1000]
BEGIN
! seeds random number generator with system time
RANDOMIZE(TIMER)
PRINT(CHR$(12);) !CLS
SUMX#=0 SUMX2#=0
FOR I%=1 TO NUM_ITEM DO
R1#=RND(1) R2#=RND(1)
Z#=SQR(-2*LOG(R1#))*COS(2*π*R2#)
A#[I%]=Z#/2+1 ! I want a normal distribution with
! mean=1 and std.dev=0.5
SUMX#+=A#[I%] SUMX2#+=A#[I%]*A#[I%]
END FOR
Z#=SUMX#/NUM_ITEM
PRINT("Average is";Z#)
PRINT("Standard dev. is";SQR(SUMX2#/NUM_ITEM-Z#*Z#))
END PROGRAM
Euler Math Toolbox
>v=normal(1,1000)*0.5+1;
>mean(v), dev(v)
1.00291801071
0.498226876528
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
F#
let n = MathNet.Numerics.Distributions.Normal(1.0,0.5)
List.init 1000 (fun _->n.Sample())
- Output:
[0.734433576; 1.54225304; 0.4407528678; 1.177675412; 0.4318617021; 0.6026656337; 0.769764924; 1.104693934; 0.6297500925; 0.9594598077; 1.684736389; 1.160376323; 0.883354356; 0.9513968363; 0.9727698268; 0.5315570949; 0.9599239266; 1.564976755; 0.7232002879; 1.084139442; 1.220914517; 0.3553085946; 1.112549824; 1.989443553; 0.5752307543; 1.156682549; 0.7886670467; 0.02050745923; 1.532060208; 1.18789591; 1.408946777; 1.038812004; 1.724679503; 1.671565045; 1.266831442; 1.363611654; 1.705819067; 0.5772366328; 0.4503488498; 1.496891481; 0.9831877282; 0.3845460366; 0.8253240671; 1.298969969; 0.4265904553; 0.9303696876; 0.445003361; 0.753175816; 0.6143534043; 1.059982235; 0.7143206784; 0.2233328038; 1.005178481; 0.7697392436; 0.5904948577; 0.5127953044; 0.6467346747; 0.7929387604; -0.1501790761; 0.8750780903; 0.941704369; 1.37941579; 0.4739006145; 1.998886344; 1.219428519; 0.06270791476; 1.097739804; 0.7584232803; 1.042177217; 1.166561247; 1.502357164; 1.171525776; 0.1528807432; 0.2289389756; 1.36208422; 0.3714421124; 1.299571092; 1.171553369; 1.317807265; 1.616662281; 1.724223246; 1.059580642; 1.270520918; -0.1827677907; 1.938593232; 1.420362143; 1.888357595; 0.7851629936; 0.7080554899; 0.7747215818; 1.403719877; 0.5765950249; 1.275206565; 0.6292054813; 1.525562798; 0.6224640457; 0.8524078517; 0.7646595627; 0.6799834691; 0.773111053; ...]
==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() ]
Factor
USING: random ;
1000 [ 1.0 0.5 normal-random-float ] replicate
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
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) }
}
}
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
For newer versions of gforth (tested on 0.7.3), it seems you need to use HERE SEED ! instead of HERE TO SEED, because SEED has been made a variable instead of a value.
rnd rnd dabs d>f
is necessary, but surprising and definitely not well documented / perhaps not compliant.
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
Free Pascal
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;
Frink
a = new array[[1000], {|x| randomGaussian[1, 0.5]}]
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.
package main
import (
"fmt"
"math"
"math/rand"
"strings"
"time"
)
const mean = 1.0
const stdv = .5
const n = 1000
func main() {
var list [n]float64
rand.Seed(time.Now().UnixNano())
for i := range list {
list[i] = mean + stdv*rand.NormFloat64()
}
// show computed mean and stdv of list
var s, sq float64
for _, v := range list {
s += v
}
cm := s / n
for _, v := range list {
d := v - cm
sq += d * d
}
fmt.Printf("mean %.3f, stdv %.3f\n", cm, math.Sqrt(sq/(n-1)))
// show histogram by hdiv divisions per stdv over +/-hrange stdv
const hdiv = 3
const hrange = 2
var h [1 + 2*hrange*hdiv]int
for _, v := range list {
bin := hrange*hdiv + int(math.Floor((v-mean)/stdv*hdiv+.5))
if bin >= 0 && bin < len(h) {
h[bin]++
}
}
const hscale = 10
for _, c := range h {
fmt.Println(strings.Repeat("*", (c+hscale/2)/hscale))
}
}
- Output:
mean 0.995, stdv 0.503 ** **** ****** ******** ************ ************ ************* ************ ********** ******** ***** *** **
Groovy
rnd = new Random()
result = (1..1000).inject([]) { r, i -> r << rnd.nextGaussian() }
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
Or using Data.Random from random-fu package:
replicateM 1000 $ normal 1 0.5
To print them:
import Data.Random
import Control.Monad
thousandRandomNumbers :: RVar [Double]
thousandRandomNumbers = replicateM 1000 $ normal 1 0.5
main = do
x <- sample thousandRandomNumbers
print x
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
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.
IDL
result = 1.0 + 0.5*randomn(seed,1000)
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 ...
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();
}
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
}
jq
Since jq is a purely functional language, it is convenient to define the pseudo-random number generator functions as filters whose inputs and outputs are arrays containing a "seed".
The following uses the same pseudo-random number generator as the Microsoft C Runtime (see Linear congruential generator).
'A Pseudo-Random Number Generator'
# 15-bit integers generated using the same formula as rand() from the Microsoft C Runtime.
# The random numbers are in [0 -- 32767] inclusive.
# Input: an array of length at least 2 interpreted as [count, state, ...]
# Output: [count+1, newstate, r] where r is the next pseudo-random number.
def next_rand_Microsoft:
.[0] as $count | .[1] as $state
| ( (214013 * $state) + 2531011) % 2147483648 # mod 2^31
| [$count+1 , ., (. / 65536 | floor) ] ;
'Box-Muller Method'
# Generate a single number following the normal distribution with mean 0, variance 1,
# 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]
# Output [n+1, nextstate, r]
def next_rand_normal:
def u: next_rand_Microsoft | .[2] /= 32767;
u as $u1
| ($u1 | u) as $u2
| ((( (8*(1|atan)) * $u1[2]) | cos)
* ((-2 * (($u2[2]) | log)) | sqrt)) as $r
| [ (.[0]+1), $u2[1], $r] ;
# Generate "count" arrays, each containing a random normal variate with the given mean and standard deviation.
# Input: [count, state]
# Output: [updatedcount, updatedstate, rnv]
# where "state" is a seed and "updatedstate" can be used as a seed.
def random_normal_variate(mean; sd; count):
next_rand_normal
| recurse( if .[0] < count then next_rand_normal else empty end)
| .[2] = (.[2] * sd) + mean;
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:
def summary:
length as $l | add as $sum | ($sum/$l) as $a
| reduce .[] as $x (0; . + ( ($x - $a) | .*. ))
| [ $a, (./$l | sqrt)] ;
[ [0,1] | random_normal_variate(1; 0.5; 1000) | .[2] ] | summary
- Output:
$ jq -n -c -f Random_numbers.jq [0.9932830741018853,0.4977760644490579]
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
Kotlin
// version 1.0.6
import java.util.Random
fun main(args: Array<String>) {
val r = Random()
val da = DoubleArray(1000)
for (i in 0 until 1000) da[i] = 1.0 + 0.5 * r.nextGaussian()
// now check actual mean and SD
val mean = da.average()
val sd = Math.sqrt(da.map { (it - mean) * (it - mean) }.average())
println("Mean is $mean")
println("S.D. is $sd")
}
Sample output:
- Output:
Mean is 1.0071784073168768 S.D. is 0.48567118114896807
LabVIEW
Lingo
-- Returns a random float value in range 0..1
on randf ()
n = random(the maxinteger)-1
return n / float(the maxinteger-1)
end
normal = []
repeat with i = 1 to 1000
normal.add(1 + sqrt(-2 * log(randf())) * cos(2 * PI * randf()) / 2)
end repeat
Lobster
Uses built-in rnd_gaussian
let mean = 1.0
let stdv = 0.5
let count = 1000
// stats computes a running mean and variance
// See Knuth TAOCP vol 2, 3rd edition, page 232
def stats(xs: [float]) -> float, float: // variance, mean
var M = xs[0]
var S = 0.0
var n = 1.0
for(xs.length - 1) i:
let x = xs[i + 1]
n = n + 1.0
let mm = (x - M)
M += mm / n
S += mm * (x - M)
return (if n > 0.0: S / n else: 0.0), M
def test_random_normal() -> [float]:
rnd_seed(floor(seconds_elapsed() * 1000000))
let r = vector_reserve(typeof return, count)
for (count):
r.push(rnd_gaussian() * stdv + mean)
let cvar, cmean = stats(r)
let cstdv = sqrt(cvar)
print concat_string(["Mean: ", string(cmean), ", Std.Deviation: ", string(cstdv)], "")
test_random_normal()
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 ?] []
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
M2000 Interpreter
M2000 use a Wichmann - Hill Pseudo Random Number Generator.
Module CheckIt {
Function StdDev (A()) {
\\ A() has a copy of values
N=Len(A())
if N<1 then Error "Empty Array"
M=Each(A())
k=0
\\ make sum, dev same type as A(k)
sum=A(k)-A(k)
dev=sum
\\ find mean
While M {
sum+=Array(M)
}
Mean=sum/N
\\ make a pointet to A()
P=A()
\\ subtruct from each item
P-=Mean
M=Each(P)
While M {
dev+=Array(M)*Array(M)
}
\\ as pointer to arrray
=(if(dev>0->Sqrt(dev/N), 0), Mean)
}
Function randomNormal {
\\ by default all numbers are double
\\ cos() get degrees
=1+Cos(360 * rnd) * Sqrt(-2 * Ln(rnd)) /2
}
\\ fill array calling randomNormal() for each item
Dim A(1000)<<randomNormal()
\\ we can pass a pointer to array and place it to stack of values
DisplayMeanAndStdDeviation(A()) ' mean ~ 1 std deviation ~0.5
\\ check M2000 rnd only
Dim B(1000)<<rnd
DisplayMeanAndStdDeviation(B()) ' mean ~ 0.5 std deviation ~0.28
DisplayMeanAndStdDeviation((0,0,14,14)) ' mean = 7 std deviation = 7
DisplayMeanAndStdDeviation((0,6,8,14)) ' mean = 7 std deviation = 5
DisplayMeanAndStdDeviation((6,6,8,8)) ' mean = 7 std deviation = 1
Sub DisplayMeanAndStdDeviation(A)
\\ push to stack all items of an array (need an array pointer)
Push ! StdDev(A)
\\ read from strack two numbers
Print "Mean is "; Number
Print "Standard Deviation is "; Number
End Sub
}
Checkit
Maple
with(Statistics):
Sample(Normal(1, 0.5), 1000);
or
1+0.5*ArrayTools[RandomArray](1000,1,distribution=normal);
Mathematica/Wolfram Language
Built-in function RandomReal with built-in distribution NormalDistribution as an argument:
RandomReal[NormalDistribution[1, 1/2], 1000]
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
Maxima
load(distrib)$
random_normal(1.0, 0.5, 1000);
MAXScript
arr = #()
for i in 1 to 1000 do
(
a = random 0.0 1.0
b = random 0.0 1.0
c = 1.0 + 0.5 * sqrt (-2*log a) * cos (360*b) -- Maxscript cos takes degrees
append arr c
)
Metafont
Metafont has normaldeviate
which produces pseudorandom normal distributed numbers with mean 0 and variance one. So the following complete the task:
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
MiniScript
randNormal = function(mean=0, stddev=1)
return mean + sqrt(-2 * log(rnd,2.7182818284)) * cos(2*pi*rnd) * stddev
end function
x = []
for i in range(1,1000)
x.push randNormal(1, 0.5)
end for
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
МК-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.
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.
Nanoquery
list = {0} * 1000
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()
end
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
NewLISP
(normal 1 .5 1000)
Nim
import random, stats, strformat
var rs: RunningStat
randomize()
for _ in 1..5:
for _ in 1..1000: rs.push gauss(1.0, 0.5)
echo &"mean: {rs.mean:.5f} stdDev: {rs.standardDeviation:.5f}"
- Output:
mean: 1.01294 stdDev: 0.49692 mean: 1.00262 stdDev: 0.50028 mean: 0.99878 stdDev: 0.49662 mean: 0.99830 stdDev: 0.49820 mean: 1.00658 stdDev: 0.49703
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();
}
}
}
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 ());;
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
ooRexx
version 1
/*REXX pgm gens 1,000 normally distributed #s: mean=1, standard dev.=0.5*/
pi=RxCalcPi() /* get value of 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.*/
n.g=random(0,1e5)/1e5 /* REXX gens uniform rand integers*/
End
Say ' old mean=' mean()
Say 'old standard deviation=' stddev()
Say
Do j=1 To n-1 By 2
m=j+1
/*use Box-Muller method */
_=sd*RxCalcPower(-2*RxCalcLog(n.j),.5)*RxCalcCos(2*pi*n.m,,'R')+mean
n.m=sd*RxCalcpower(-2*RxCalcLog(n.j),.5)*RxCalcSin(2*pi*n.m,,'R')+,
mean /* rand # must be 0???1. */
n.j=_
End /* j */
Say ' new mean=' mean()
Say 'new standard deviation=' stddev()
Exit
mean:
_=0
Do k=1 For n
_=_+n.k
End
Return _/n
stddev:
_avg=mean()
_=0
Do k=1 For n
_=_+(n.k-_avg)**2
End
Return RxCalcPower(_/n,.5)
:: requires rxmath library
- Output:
old mean= 0.49830002 old standard deviation= 0.283199568 new mean= 1.00377404 new standard deviation= 0.501444536
version 2
Using the nice function names in the algorithm.
/*REXX pgm gens 1,000 normally distributed #s: mean=1, standard dev.=0.5*/
pi=RxCalcPi() /* get value of 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.*/
n.g=random(0,1e5)/1e5 /* REXX gens uniform rand integers*/
End
Say ' old mean=' mean()
Say 'old standard deviation=' stddev()
Say
Do j=1 To n-1 By 2
m=j+1
/*use Box-Muller method */
_=sd*sqrt(-2*ln(n.j))*cos(2*pi*n.m)+mean
n.m=sd*sqrt(-2*ln(n.j))*sin(2*pi*n.m)+mean
n.j=_
End
Say ' new mean=' mean()
Say 'new standard deviation=' stddev()
Exit
mean:
_=0
Do k=1 For n
_=_+n.k
End
Return _/n
stddev:
_avg=mean()
_=0
Do k=1 For n
_=_+(n.k-_avg)**2
End
Return sqrt(_/n)
sqrt: Return RxCalcSqrt(arg(1))
ln: Return RxCalcLog(arg(1))
cos: Return RxCalcCos(arg(1),,'R')
sin: Return RxCalcSin(arg(1),,'R')
:: requires rxmath library
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*u2) \\ in previous version "u1" instead of "u2" was used --> has given crap distribution
\\ Could easily be extended with a second normal at very little cost.
};
vector(1000,unused,rnormal()/2+1)
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}
var
u1, u2: real;
begin
u1 := random;
u2 := random;
rnorm := mean * abs(1 + sqrt(-2 * (ln(u1))) * cos(2 * pi * u2) * sd);
/* error !?! Shouldn't it be "mean +" instead of "mean *" ? */
end;
Delphi and Free Pascal support implement a randg function that delivers Gaussian-distributed random numbers.
Perl
my $PI = 2 * atan2 1, 0;
my @nums = map {
1 + 0.5 * sqrt(-2 * log rand) * cos(2 * $PI * rand)
} 1..1000;
Phix
function RandomNormal() return sqrt(-2*log(rnd())) * cos(2*PI*rnd()) end function sequence s = repeat(0,1000) for i=1 to length(s) do s[i] = 1 + 0.5 * RandomNormal() end for
Phixmonti
include ..\Utilitys.pmt
def RandomNormal
drop rand log -2 * sqrt 2 pi * rand * cos * 0.5 * 1 +
enddef
1000 var n
0 n repeat
getid RandomNormal map
dup
sum n / var mean
"Mean: " print mean print nl
0 swap n for
get mean - 2 power rot + swap
endfor
swap n / sqrt "Standard deviation: " print print
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);
}
Picat
main =>
_ = random2(), % random seed
G = [gaussian_dist(1,0.5) : _ in 1..1000],
println(first_10=G[1..10]),
println([mean=avg(G),stdev=stdev(G)]),
nl.
% Gaussian (Normal) distribution, Box-Muller algorithm
gaussian01() = Y =>
U = frand(0,1),
V = frand(0,1),
Y = sqrt(-2*log(U))*sin(2*math.pi*V).
gaussian_dist(Mean,Stdev) = Mean + (gaussian01() * Stdev).
% Variance of Xs
variance(Xs) = Variance =>
Mu = avg(Xs),
N = Xs.len,
Variance = sum([ (X-Mu)**2 : X in Xs ]) / N.
% Standard deviation
stdev(Xs) = sqrt(variance(Xs)).
- Output:
first_10 = [1.639965415776091,0.705425965005482,0.981532402477848,0.309148743347499,1.252800181962738,0.098829881195179,0.74888084504147,0.181494956495445,1.304931340021904,0.595939453660087] [mean = 0.99223677282248,stdev = 0.510336641737154]
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
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
PL/SQL
DECLARE
--The desired collection
type t_coll is table of number index by binary_integer;
l_coll t_coll;
c_max pls_integer := 1000;
BEGIN
FOR l_counter IN 1 .. c_max
LOOP
-- dbms_random.normal delivers normal distributed random numbers with a mean of 0 and a variance of 1
-- We just adjust the values and get the desired result:
l_coll(l_counter) := DBMS_RANDOM.normal * 0.5 + 1;
DBMS_OUTPUT.put_line (l_coll(l_counter));
END LOOP;
END;
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;
PowerShell
Equation adapted from Liberty BASIC
function Get-RandomNormal
{
[CmdletBinding()]
Param ( [double]$Mean, [double]$StandardDeviation )
$RandomNormal = $Mean + $StandardDeviation * [math]::Sqrt( -2 * [math]::Log( ( Get-Random -Minimum 0.0 -Maximum 1.0 ) ) ) * [math]::Cos( 2 * [math]::PI * ( Get-Random -Minimum 0.0 -Maximum 1.0 ) )
return $RandomNormal
}
# Standard deviation function for testing
function Get-StandardDeviation
{
[CmdletBinding()]
param ( [double[]]$Numbers )
$Measure = $Numbers | Measure-Object -Average
$PopulationDeviation = 0
ForEach ($Number in $Numbers) { $PopulationDeviation += [math]::Pow( ( $Number - $Measure.Average ), 2 ) }
$StandardDeviation = [math]::Sqrt( $PopulationDeviation / ( $Measure.Count - 1 ) )
return $StandardDeviation
}
# Test
$RandomNormalNumbers = 1..1000 | ForEach { Get-RandomNormal -Mean 1 -StandardDeviation 0.5 }
$Measure = $RandomNormalNumbers | Measure-Object -Average
$Stats = [PSCustomObject]@{
Count = $Measure.Count
Average = $Measure.Average
StandardDeviation = Get-StandardDeviation -Numbers $RandomNormalNumbers
}
$Stats | Format-List
- Output:
Count : 1000 Average : 1.01206560135809 StandardDeviation : 0.489099623426272
Python
- Using random.gauss
>>> import random
>>> values = [random.gauss(1, .5) for i in range(1000)]
>>>
- Quick check of distribution
>>> def quick_check(numbers):
count = len(numbers)
mean = sum(numbers) / count
sdeviation = (sum((i - mean)**2 for i in numbers) / count)**0.5
return mean, sdeviation
>>> 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
>>> values = [ random.normalvariate(1, 0.5) for i in range(1000)]
>>> quick_check(values)
(0.990099111944864, 0.5029847005836282)
>>>
R
# For reproducibility, set the seed:
set.seed(12345L)
result <- rnorm(1000, mean = 1, sd = 0.5)
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)
Raku
(formerly Perl 6)
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;
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
ReScript
let pi = 4.0 *. atan(1.0)
let random_gaussian = () => {
1.0 +.
sqrt(-2.0 *. log(Random.float(1.0))) *.
cos(2.0 *. pi *. Random.float(1.0))
}
let a = Belt.Array.makeBy(1000, (_) => random_gaussian ())
for i in 0 to 10 {
Js.log(a[i])
}
REXX
The REXX language doesn't have any "higher math" functions like SQRT/SIN/COS/LN/LOG/EXP/POW/etc.,
so we hoi polloi REXX programmers have to roll our own.
Programming note: note the range of the random numbers: (0,1]
(that is, random numbers from zero──►unity, excluding zero, including unity).
/*REXX pgm generates 1,000 normally distributed numbers: mean=1, standard deviation=½.*/
numeric digits 20 /*the default decimal digit precision=9*/
parse arg n seed . /*allow specification of N and the seed*/
if n=='' | n=="," then n=1000 /*N: is the size of the array. */
if datatype(seed,'W') then call random ,,seed /*SEED: for repeatable random numbers. */
newMean=1 /*the desired new mean (arithmetic avg)*/
sd=1/2 /*the desired new standard deviation. */
do g=1 for n /*generate N uniform random #'s (0,1].*/
#.g = random(1, 1e5) / 1e5 /*REXX's RANDOM BIF generates integers.*/
end /*g*/ /* [↑] random integers ──► fractions. */
say ' old mean=' mean()
say 'old standard deviation=' stdDev()
call pi; pi2=pi * 2 /*define pi and also 2 * pi. */
say
do j=1 to n-1 by 2; m=j+1 /*step through the iterations by two. */
_=sd * sqrt(ln(#.j) * -2) /*calculate the used-twice expression.*/
#.j=_ * cos(pi2 * #.m) + newMean /*utilize the Box─Muller method. */
#.m=_ * sin(pi2 * #.m) + newMean /*random number must be: (0,1] */
end /*j*/
say ' new mean=' mean()
say 'new standard deviation=' stdDev()
exit /*stick a fork in it, we're all done. */
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/
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 /*digs overkill*/
pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862; return pi /* " " */
r2r: return arg(1) // (pi() * 2) /*normalize ang*/
sin: procedure; parse 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
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
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/
cos: procedure; parse arg x; x=r2r(x); a=abs(x); hpi=pi * .5
numeric fuzz min(6, digits() - 3); if a=pi then return -1
if a=hpi | a=hpi*3 then return 0; if a=pi/3 then return .5
if a=pi * 2/3 then return -.5; return .sinCos(1,1,-1)
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits; h=d+6
numeric form; parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g * .5'e'_ %2
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*/
numeric digits d; return g/1
output when using the default inputs:
old mean= 0.5015724 old standard deviation= 0.28652466389342471402 new mean= 0.98807025356443262689 new standard deviation= 0.50002924192766720838
Ring
for i = 1 to 10
see random(i) + nl
next i
RPL
≪ RAND LN NEG 2 * √ RAND 2 * π * COS * →NUM 2 / 1 + ≫ 'RANDN' STO ≪ CL∑ 1 1000 START RANDN ∑+ NEXT MEAN PSDEV ≫ 'TASK' STO
- Output:
1: .990779804949 2: .487204045227
The collection is stored in a predefined array named ∑DAT
, which is automatically created/updated when using the ∑+
instruction and remains available until the user decides to purge it, typically by calling the CL∑
command.
Ruby
Array.new(1000) { 1 + Math.sqrt(-2 * Math.log(rand)) * Math.cos(2 * Math::PI * rand) }
Rust
Using a for-loop:
extern crate rand;
use rand::distributions::{Normal, IndependentSample};
fn main() {
let mut rands = [0.0; 1000];
let normal = Normal::new(1.0, 0.5);
let mut rng = rand::thread_rng();
for num in rands.iter_mut() {
*num = normal.ind_sample(&mut rng);
}
}
Using iterators:
extern crate rand;
use rand::distributions::{Normal, IndependentSample};
fn main() {
let rands: Vec<_> = {
let normal = Normal::new(1.0, 0.5);
let mut rng = rand::thread_rng();
(0..1000).map(|_| normal.ind_sample(&mut rng)).collect()
};
}
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 ----------------------------
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;
Scala
One liner
List.fill(1000)(1.0 + 0.5 * scala.util.Random.nextGaussian)
Academic
object RandomNumbers extends App {
val distribution: LazyList[Double] = {
def randomNormal: Double = 1.0 + 0.5 * scala.util.Random.nextGaussian
def normalDistribution(a: Double): LazyList[Double] = a #:: normalDistribution(randomNormal)
normalDistribution(randomNormal)
}
/*
* Let's test it
*/
def calcAvgAndStddev[T](ts: Iterable[T])(implicit num: Fractional[T]): (T, Double) = {
val mean: T =
num.div(ts.sum, num.fromInt(ts.size)) // Leaving with type of function T
// Root of mean diffs
val stdDev = Math.sqrt(ts.map { x =>
val diff = num.toDouble(num.minus(x, mean))
diff * diff
}.sum / ts.size)
(mean, stdDev)
}
println(calcAvgAndStddev(distribution.take(1000))) // e.g. (1.0061433267806525,0.5291834867560893)
}
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)
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;
Sidef
var arr = 1000.of { 1 + (0.5 * sqrt(-2 * 1.rand.log) * cos(Num.tau * 1.rand)) }
arr.each { .say }
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 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.
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.
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 uniformdeviate = fn seed =>
let
val in31m = (Real.fromInt o Int32.toInt ) (getOpt (Int32.maxInt,0) );
val in31 = in31m +1.0;
val s1 = 41160.0;
val s2 = 950665216.0;
val v = Real.realFloor seed;
val val1 = v*s1;
val val2 = v*s2;
val next1 = Real.fromLargeInt (Real.toLargeInt IEEEReal.TO_NEGINF (val1/in31)) ;
val next2 = Real.rem(Real.realFloor(val2/in31) , in31m );
val valt = val1+val2 - (next1+next2)*in31m;
val nextt = Real.realFloor(valt/in31m);
val valt = valt - nextt*in31m;
in
(valt/in31m,valt)
end;
val store = ref (0.0,0.0);
val rec u = fn S => fn 0 => [] | n=> (store:=uniformdeviate S; (#1 (!store)):: (u (#2 (!store)) (n-1))) ;
in
u seed n
end;
local
open Math
in
val bmconv = fn urand => fn vrand => 1.0+0.5*(sqrt(~2.0*ln urand)*cos (2.0*pi*vrand) )
end;
val rec makeNormals = fn once => fn u::v::[] => [once u v] |
u::v::rm => (once u v )::(makeNormals once rm );
val anyrealseed=1009.0 ;
makeNormals bmconv (urandomlist anyrealseed 2000);
Stata
clear all
set obs 1000
gen x=rnormal(1,0.5)
Mata
a = rnormal(1000,1,1,0.5)
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()}]
}
TorqueScript
for (%i = 0; %i < 1000; %i++)
%list[%i] = 1 + mSqrt(-2 * mLog(getRandom())) * mCos(2 * $pi * getRandom());
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)>
Visual FoxPro
LOCAL i As Integer, m As Double, n As Integer, sd As Double
py = PI()
SET TALK OFF
SET DECIMALS TO 6
CREATE CURSOR gdev (deviate B(6))
RAND(-1)
n = 1000
m = 1
sd = 0.5
CLEAR
FOR i = 1 TO n
INSERT INTO gdev VALUES (GaussDev(m, 1/sd))
ENDFOR
CALCULATE AVG(deviate), STD(deviate) TO m, sd
? "Mean", m, "Std Dev", sd
SET TALK ON
SET DECIMALS TO
FUNCTION GaussDev(mean As Double, sdev As Double) As Double
LOCAL z As Double
z = SQRT(-2*LOG(RAND()))*COS(py*RAND())
IF sdev # 0
z = mean + z/sdev
ENDIF
RETURN z
ENDFUNC
V (Vlang)
import crypto.rand
fn main() {
mut nums := []u64{}
for _ in 0..1000 {
nums << rand.int_u64(10000) or {0} // returns random unsigned 64-bit integer from real OS source of entropy
}
println(nums)
}
Wren
import "random" for Random
var rand = Random.new()
var randNormal = Fn.new { (-2 * rand.float().log).sqrt * (2 * Num.pi * rand.float()).cos }
var stdDev = Fn.new { |a, m|
var c = a.count
return ((a.reduce(0) { |acc, x| acc + x*x } - m*m*c) / c).sqrt
}
var n = 1000
var numbers = List.filled(n, 0)
var mu = 1
var sigma = 0.5
var sum = 0
for (i in 0...n) {
numbers[i] = mu + sigma*randNormal.call()
sum = sum + numbers[i]
}
var mean = sum / n
System.print("Actual mean : %(mean)")
System.print("Actual std dev: %(stdDev.call(numbers, mean))")
- Output:
Sample run:
Actual mean : 1.0053988699746 Actual std dev: 0.4961645117026
XPL0
define PI = 3.14159265358979323846;
func real DRand; \Uniform distribution, [0..1]
return float(Ran(1_000_000)) / 1e6;
func real RandomNormal; \Normal distribution, centered on 0, std dev 1
return sqrt(-2.*Log(DRand)) * Cos(2.*PI*DRand);
int I;
real Rands(1000);
for I:= 0 to 1000-1 do
Rands(I):= 1.0 + 0.5*RandomNormal
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
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
- Programming Tasks
- Basic language learning
- Probability and statistics
- Randomness
- GUISS/Omit
- UNIX Shell/Omit
- Ada
- ALGOL 68
- Arturo
- AutoHotkey
- Avail
- AWK
- BASIC
- ANSI BASIC
- Applesoft BASIC
- BASIC256
- BBC BASIC
- Chipmunk Basic
- Commodore BASIC
- FreeBASIC
- FutureBasic
- GW-BASIC
- Liberty BASIC
- Minimal BASIC
- PureBasic
- QuickBASIC
- Run BASIC
- TI-83 BASIC
- ZX Spectrum Basic
- C
- C sharp
- C++
- Boost
- Clojure
- COBOL
- Common Lisp
- Crystal
- D
- Tango
- Delphi
- DWScript
- E
- EasyLang
- Eiffel
- Elena
- Elixir
- Erlang
- ERRE
- Euler Math Toolbox
- Euphoria
- F Sharp
- Factor
- Falcon
- Fantom
- Forth
- Fortran
- Free Pascal
- Frink
- Go
- Groovy
- Haskell
- HicEst
- Icon
- Unicon
- IDL
- J
- Java
- JavaScript
- Jq
- Julia
- Kotlin
- LabVIEW
- Lingo
- Lobster
- Logo
- Lua
- M2000 Interpreter
- Maple
- Mathematica
- Wolfram Language
- MATLAB
- Maxima
- MAXScript
- Metafont
- MiniScript
- Mirah
- МК-61/52
- Modula-3
- Nanoquery
- NetRexx
- NewLISP
- Nim
- Objeck
- OCaml
- Octave
- OoRexx
- PARI/GP
- Pascal
- Perl
- Phix
- Phixmonti
- PHP
- Picat
- PicoLisp
- PL/I
- PL/SQL
- Pop11
- PowerShell
- Python
- R
- Racket
- Raku
- Raven
- ReScript
- REXX
- Ring
- RPL
- Ruby
- Rust
- Rand
- SAS
- Sather
- Scala
- Scheme
- Seed7
- Sidef
- Standard ML
- Stata
- Tcl
- TorqueScript
- Ursala
- Visual FoxPro
- V (Vlang)
- Wren
- XPL0
- Yorick
- Zkl