Statistics/Normal distribution: Difference between revisions
SqrtNegInf (talk | contribs) m (→{{header|Perl}}: Fix link: Perl 6 --> Raku) |
m (→{{header|Wren}}: Minor tidy) |
||
(16 intermediate revisions by 14 users not shown) | |||
Line 13: | Line 13: | ||
=={{header|C}}== |
=={{header|C}}== |
||
<syntaxhighlight lang="c">/* |
|||
<lang C>/* |
|||
* RosettaCode example: Statistics/Normal distribution in C |
* RosettaCode example: Statistics/Normal distribution in C |
||
* |
* |
||
Line 139: | Line 139: | ||
} |
} |
||
return EXIT_FAILURE; |
return EXIT_FAILURE; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>mean = 0.000477941, stddev = 0.999945 |
<pre>mean = 0.000477941, stddev = 0.999945 |
||
Line 208: | Line 208: | ||
=={{header|C sharp|C#}}== |
=={{header|C sharp|C#}}== |
||
{{libheader|Math.Net}} |
{{libheader|Math.Net}} |
||
< |
<syntaxhighlight lang="csharp">using System; |
||
using MathNet.Numerics.Distributions; |
using MathNet.Numerics.Distributions; |
||
using MathNet.Numerics.Statistics; |
using MathNet.Numerics.Statistics; |
||
Line 239: | Line 239: | ||
RunNormal(10000); |
RunNormal(10000); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Sample size: 100 |
<pre>Sample size: 100 |
||
Line 285: | Line 285: | ||
=={{header|C++}}== |
=={{header|C++}}== |
||
showing features of C++11 here |
showing features of C++11 here |
||
< |
<syntaxhighlight lang="cpp">#include <random> |
||
#include <map> |
#include <map> |
||
#include <string> |
#include <string> |
||
Line 315: | Line 315: | ||
} |
} |
||
return 0 ; |
return 0 ; |
||
}</ |
}</syntaxhighlight> |
||
Output: |
Output: |
||
<pre>The mean of the distribution is 1 , the standard deviation 1 ! |
<pre>The mean of the distribution is 1 , the standard deviation 1 ! |
||
Line 347: | Line 347: | ||
=={{header|D}}== |
=={{header|D}}== |
||
This uses the Box-Muller method as in the Go entry, and the module from the Statistics/Basic. A ziggurat-based normal generator for the Phobos standard library is in the works. |
This uses the Box-Muller method as in the Go entry, and the module from the Statistics/Basic. A ziggurat-based normal generator for the Phobos standard library is in the works. |
||
< |
<syntaxhighlight lang="d">import std.stdio, std.random, std.math, std.range, std.algorithm, |
||
statistics_basic; |
statistics_basic; |
||
Line 373: | Line 373: | ||
writefln("Mean: %8.6f, SD: %8.6f\n", data.meanStdDev[]); |
writefln("Mean: %8.6f, SD: %8.6f\n", data.meanStdDev[]); |
||
data.map!q{ max(0.0, min(0.9999, a / 3 + 0.5)) }.showHistogram01; |
data.map!q{ max(0.0, min(0.9999, a / 3 + 0.5)) }.showHistogram01; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Mean: 0.000528, SD: 0.502245 |
<pre>Mean: 0.000528, SD: 0.502245 |
||
Line 387: | Line 387: | ||
0.8: ****** |
0.8: ****** |
||
0.9: *</pre> |
0.9: *</pre> |
||
=={{header|EDSAC order code}}== |
|||
<syntaxhighlight lang="edsac"> |
|||
[Normal distribution for Rosetta Code |
|||
EDSAC program, Initial Orders 2] |
|||
[================================================================================== |
|||
Uses an accept-reject method, which requires only logarithms (no trig functions). |
|||
Let u, v be independent uniform variates in (0, 1). Let x = -ln(u), y = -ln(v). |
|||
Accept x iff (x - 1)^2 <= 2y. If x is accepted, negate it with probability 1/2. |
|||
Then x is normally distributed with mean 0 and standard deviation 1. |
|||
The algorithm is modified for this EDSAC version: |
|||
(1) Uses EDSAC library subroutine L1 to calculate (1/32)log_2() rather than ln() |
|||
(2) Since real numbers on EDSAC are restricted to the interval [-1, 1), scales so |
|||
that the standard deviation is 1/4, and reject values >= 4 s.d. from the mean. |
|||
In the histogram, counts are divided by 16, with rounding. |
|||
On the EDSAC PC simulator, takes 2.75 EDSAC hours to find 4096 normal variates. |
|||
[=================================================================================] |
|||
[Arrange the storage] |
|||
T46K P70F [N parameter: library subroutine P1 to print +ve fraction] |
|||
T47K P200F [M parameter: main routine and dependent subroutines] |
|||
T49K P92F [L parameter: library subroutine L1 to calculate log_2] |
|||
T51K P130F [G parameter: generator for pseudo-random numbers] |
|||
T52K P168F [A parameter: library subroutine S2 for square root] |
|||
T54K P190F [C parameter: constants read in by library subroutine R9] |
|||
[Library subroutine R9, reads non-negative integers at load time. |
|||
Fractions can be read by converting to integers (multiply by 2^34). |
|||
15 locations, must be loaded at location 56.] |
|||
T56KGKT20FVDL8FA40DUDTFI40FA40FS39FG@S2FG23FA5@T5@E4@ |
|||
[Library subroutine M3, prints header at load time. |
|||
M3 and header are then overwritten.] |
|||
PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF |
|||
*MEAN#!K0BL!!*SD#!K0M25BL@&#..PZ |
|||
[========================== C parameter ==========================] |
|||
[Tell R9 where to store integers read from tape] |
|||
E69K T#C ['T m D' in WWG, but this also works] |
|||
[List of integers; separated by F; end list with #TZ] |
|||
123456789F5F1549082005F11908177887F858993#TZ |
|||
[0] [seed for random number generator] |
|||
[2] [minimum argument for library subroutine L1] |
|||
[4] [1/(16*ln(2)) for accept-reject] |
|||
[6] [ln(2) for scaling standard deviation] |
|||
[8] [0.00005 for rounding in print routine] |
|||
[========================== M parameter ==========================] |
|||
E25K TM GK |
|||
[0] P4096F [number of data, in address field; code below assumes 4096] |
|||
[1] PF [negative count of data] |
|||
[2] PF [worlspace, low word] |
|||
[3] PF [workspace, high word] |
|||
[4] PF PF [auxiliary for accept-reject algorithm] |
|||
[6] PF PF [sum of value/count, for mean] |
|||
[8] PF PF [sum of value^2/count, for variance] |
|||
[10] PFPFPFPFPFPFPFPFPFPFPFPFPFPFPFPF [16 bins for histogram] |
|||
[26] CF [11110...0 binary, to isolate bin index; also prints colon] |
|||
[27] A10@ [A order for bin{0}] |
|||
[28] AF [A order for bin{16}] |
|||
[29] P8F [8 in address field] |
|||
[30] MF [subtract from A order to make T order; also prints dot] |
|||
[31] WF [1/8] |
|||
[32] FF [-15/16, mid value of histogram bin{0}] |
|||
[33] PF [mid value of current bin] |
|||
[Teleprinter] |
|||
[34] #F [set figures mode] |
|||
[35] PF |
|||
[36] AF [minus (in figures mode)] |
|||
[37] !F [space] |
|||
[38] @F [carriage return] |
|||
[39] &F [line feed] |
|||
[Enter with acc = 0] |
|||
[40] |
|||
S@ T1@ [store negated data count in address field] |
|||
A#C T4D [copy seed to 4D for PRNG] |
|||
[44] A44@ GG [initialize PRNG] |
|||
T6#@ T8#@ [clear sum and sum of squares] |
|||
O34@ [set teleprinter to figures] |
|||
[Start of loop to generate normal variates] |
|||
[49] TF [clear acc] |
|||
[Calculate and store logarithms of uniform variates] |
|||
[50] A50@ G176@ SD T2#@ [corresponds to x = -ln(u)] |
|||
[54] A54@ G176@ SD T4#@ [corresponds to y = -ln(v)] |
|||
[Accept or reject] |
|||
A4#C RD [acc := (1/32*ln(2))] |
|||
S2#@ [subtract x] |
|||
TD HD ND [acc := negated square] |
|||
H4#C V4#@ [add 2*(auxiliary value)] |
|||
G49@ [reject if result < 0] |
|||
[First log is accepted; multiply by 8*ln(2) so that s.d. becomes 0.25] |
|||
[Reject values outside (-1,1) (that is >= 4 s.d. from mean, unlikely)] |
|||
T6F [clear acc] |
|||
H6#C V2#@ [times ln(2)] |
|||
S31@ E49@ [if product >= 1/8 (unlikely) reject and try again] |
|||
A31@ [restore acc after test] |
|||
L2F [shift 3 left to complete scaling] |
|||
YF T2#@ [round and store back] |
|||
[Finally change sign with probability 1/2] |
|||
T6F T4D A2F T4F [pass range = 2 to PRNG] |
|||
[80] A80@ G1G [call PRNG, 0 or 1 returned in 0D] |
|||
SD E87@ [don't change sign if 0D = 0] |
|||
T6F S2#@ T2#@ [change sign, so value < 0] |
|||
[87] [Here with acc = 0.] |
|||
[To print the values, replace the following jump with a no-op (X F)] |
|||
E94@ |
|||
A2#@ TD [pass abs(value) to print routine] |
|||
[90] A90@ G191@ |
|||
O38@ O39@ [print CR, LF] |
|||
[94] A2#@ R1024F YF [load value, divide by count, round] |
|||
A6#@ T6#@ [update sum] |
|||
H2#@ V2#@ R1024F YF |
|||
A8#@ T8#@ [update sum of squares] |
|||
[Inc count in appropriate bin] |
|||
H26@ [mult reg := mask to isolate bin index] |
|||
C3@ [acc := top 4 bits of value] |
|||
R1024F [12 right, get bin -8..7 in address field] |
|||
A29@ [add 8 to address, bin index now 0..15] |
|||
A27@ [make A order for bin] |
|||
U113@ [plant in code] |
|||
S30@ [convert to T order] |
|||
T115@ [plant in code] |
|||
[113] AF [(planted) acc := count in bin] |
|||
A2F [inc by 1 in address field] |
|||
[115] TF [(planted) store updated bin count] |
|||
A1@ A2F U1@ [inc negative count of variates] |
|||
G49@ [loop till got required number] |
|||
[Print mean and standard deviation] |
|||
A6#@ TD [pass mean to print subroutine] |
|||
[122] A122@ G191@ |
|||
O37@ O37@ O37@ [print spaces] |
|||
A8#@ H6#@ |
|||
N6#@ T4D [calc variance, pass to square root s/r] |
|||
[131] A131@ GA [4D := standard deviation] |
|||
A4D U8#@ TD [pass to print subroutine] |
|||
[136] A136@ G191@ [print s.d.] |
|||
O38@ O39@ O38@ O39@ [print CR, LF twice] |
|||
[Print histogram] |
|||
A29@ LD A27@ [make A order for exclusive end bin] |
|||
T28@ [store as test for end] |
|||
A32@ T33@ [initialize mid-value of bin] |
|||
A27@ [A order for bin{0}] |
|||
[149] T158@ [loop: plant A order for current bin] |
|||
TD A33@ U1F [0D = mid-value of bin, extended to 35 bits] |
|||
A31@ T33@ [update mid-value for next time] |
|||
[155] A155@ G191@ [print mid-value] |
|||
O26@ [print colon] |
|||
[158] AF [(planted) load number of hits in address field] |
|||
A29@ [add 8 for rounding] |
|||
R4F [divide by 16] |
|||
E163@ [jump to middle of loop] |
|||
[162] O175@ [print plus sign] |
|||
[163] S2F E162@ [loop till printed enough plus signs] |
|||
O38@ O39@ [print CR, LF] |
|||
TF [clear acc] |
|||
A158@ A2F [make A order for next bin] |
|||
S28@ [compare with end order] |
|||
E174@ [exit if no more bins] |
|||
A28@ [restore acc after test] |
|||
G149@ [loop back (A order is negative)] |
|||
[174] O34@ [dummy character to flush print buffer] |
|||
[175] ZF [halt program; also serves as plus sign] |
|||
[Subroutine of main routine. Sets 0D := (1/32)log_2 of uniform variate] |
|||
[176] A3F T188@ [plant return link as usual] |
|||
[178] T4D [pass range = 0 to PRNG] |
|||
[179] A179@ G1G [0D := uniform variate] |
|||
AD S2#C [test for too small (unlikely)] |
|||
G189@ [jump to try again if so] |
|||
A2#C T6D [OK, pass uniform variate to logarithm subroutine] |
|||
[186] A186@ GL [0D := logarithm] |
|||
[188] ZF [(planted) return to caller] |
|||
[189] TF E178@ [if failed, clear acc and try again] |
|||
[Subroutine of main routine; prints signed fraction in 0D to 4 decimals.] |
|||
[Wrapper for library subroutine P1; adds sign, rounding, layout.] |
|||
[191] A3F T207@ [plant return link as usual] |
|||
AD G197@ [acc := value, jump if < 0] |
|||
O37@ E200@ [print space, jump to common code] |
|||
[197] O36@ [value < 0, print minus] |
|||
TD SD [acc := abs(valus)] |
|||
[200] A8#C TD [add 0.00005 for rounding, pass to P1 in 0D] |
|||
O35@ O30@ [print '0.'] |
|||
[204] A204@ GN P4F [call P1 to print 5 decimals] |
|||
[207] ZF [(planted) jump back to caller] |
|||
[========================== G parameter ==========================] |
|||
[Linear congruential generator, same algorithm as Delphi 7 LCG.] |
|||
[38 storage locations.] |
|||
[Initialize: Call 0G with 35-bit seed in 4D.] |
|||
[Input: Call 1G with 35-bit range in 4D.] |
|||
[Output: if range > 0, 0D := random 35-bit integer in 0..(range - 1).] |
|||
[if range = 0, 0D := random 35-bit real in [0, 1)] |
|||
E25K TG |
|||
GKG10@G15@T2#ZPFT2ZI514DP257FT4#ZPFT4ZPDPFT6#ZPFT6ZPFRFA6#@S4#@ |
|||
T6#@E25FE8ZPFT8ZPFPFA3FT14@A4DT8#@ZFA3FT37@H2#@V8#@L512FL512FL1024F |
|||
A4#@T8#@H6#@C8#@T8#@S4DG32@TDA8#@E35@H4DTDV8#@L1FTDZF |
|||
[==================== Library subroutines ============================] |
|||
E25K TL |
|||
[L1: Logarithm to base 2.] |
|||
[0D := (1/32)*log_2(6D), provided 6D > 2^(-32)] |
|||
[38 storage locations; working positions 4D and 8F.] |
|||
GKA3FT33@E11@IFP1024FP512FA3@LDT6DADS4@TDS3@A6DG6@T8FS5@T4D |
|||
H6DV6DS3@E34@A3@LDYFT6DA4DADTDA4DRDYFG17@EFA3@YFT6DE29@ |
|||
E25K TN |
|||
[P1: Print non-negative fraction in 0D, without layout or rounding] |
|||
[21 storage locations.] |
|||
GKA18@U17@S20@T5@H19@PFT5@VDUFOFFFSFL4FTDA5@A2FG6@EFU3FJFM1F |
|||
E25K TA |
|||
[S2: square root.] |
|||
[Closed: 22 storage locations, working position 0D.] |
|||
[Forms sqrt( C(4D)) where C(4D) > 0, and places result in 4D.] |
|||
GKA3FT20@A4DS9@A6@UDHDR1FS21@TDN4DRDA4DYFT4DVDTDVDYFG5@EFSF |
|||
[======================= M parameter again ======================] |
|||
E25K TM GK |
|||
E40Z [define entry point] |
|||
PF [acc = 0 on entry] |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
MEAN (0?) SD (0.25?) |
|||
-0.0015 0.2488 |
|||
-0.9375: |
|||
-0.8125:+ |
|||
-0.6875:++ |
|||
-0.5625:++++ |
|||
-0.4375:+++++++++++ |
|||
-0.3125:++++++++++++++++++++++ |
|||
-0.1875:++++++++++++++++++++++++++++++++++++++++ |
|||
-0.0625:++++++++++++++++++++++++++++++++++++++++++++++++ |
|||
0.0625:++++++++++++++++++++++++++++++++++++++++++++++++++ |
|||
0.1875:++++++++++++++++++++++++++++++++++++++++ |
|||
0.3125:+++++++++++++++++++++++ |
|||
0.4375:+++++++++++ |
|||
0.5625:++++ |
|||
0.6875:+ |
|||
0.8125: |
|||
0.9375: |
|||
</pre> |
|||
=={{header|Elixir}}== |
=={{header|Elixir}}== |
||
< |
<syntaxhighlight lang="elixir">defmodule Statistics do |
||
def normal_distribution(n, w\\5) do |
def normal_distribution(n, w\\5) do |
||
{sum, sum2, hist} = generate(n, w) |
{sum, sum2, hist} = generate(n, w) |
||
Line 418: | Line 663: | ||
Enum.each([100,1000,10000], fn n -> |
Enum.each([100,1000,10000], fn n -> |
||
Statistics.normal_distribution(n) |
Statistics.normal_distribution(n) |
||
end)</ |
end)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 518: | Line 763: | ||
2.4: == |
2.4: == |
||
2.6: = |
2.6: = |
||
</pre> |
|||
=={{header|Factor}}== |
|||
{{works with|Factor|0.99 2020-01-23}} |
|||
<syntaxhighlight lang="factor">USING: assocs formatting kernel math math.functions |
|||
math.statistics random sequences sorting ; |
|||
2,000,000 [ 0 1 normal-random-float ] replicate ! make data set |
|||
dup [ mean ] [ population-std ] bi ! calculate and show |
|||
"Mean: %f\nStdev: %f\n\n" printf ! mean and stddev |
|||
[ 10 * floor 10 / ] map ! map data to buckets |
|||
histogram >alist [ first ] sort-with ! create histogram sorted by bucket (key) |
|||
dup values supremum ! find maximum count |
|||
[ |
|||
[ /f 100 * >integer ] keepd ! how big should this histogram bar be? |
|||
[ [ CHAR: * ] "" replicate-as ] dip ! make the bar |
|||
"% 5.2f: %s %d\n" printf ! print a line of the histogram |
|||
] curry assoc-each</syntaxhighlight> |
|||
{{out}} |
|||
<pre style="font-size:80%; height: 120ex; overflow: scroll"> |
|||
Mean: 0.000798 |
|||
Stdev: 1.000549 |
|||
-4.90: 2 |
|||
-4.80: 1 |
|||
-4.70: 1 |
|||
-4.60: 3 |
|||
-4.50: 3 |
|||
-4.40: 6 |
|||
-4.30: 15 |
|||
-4.20: 13 |
|||
-4.10: 16 |
|||
-4.00: 42 |
|||
-3.90: 62 |
|||
-3.80: 68 |
|||
-3.70: 98 |
|||
-3.60: 145 |
|||
-3.50: 205 |
|||
-3.40: 311 |
|||
-3.30: 379 |
|||
-3.20: 580 |
|||
-3.10: 739 |
|||
-3.00: * 1002 |
|||
-2.90: * 1349 |
|||
-2.80: ** 1893 |
|||
-2.70: *** 2499 |
|||
-2.60: **** 3211 |
|||
-2.50: ***** 4035 |
|||
-2.40: ****** 5141 |
|||
-2.30: ******* 6392 |
|||
-2.20: ********* 7869 |
|||
-2.10: ************ 9780 |
|||
-2.00: ************** 11787 |
|||
-1.90: ****************** 14483 |
|||
-1.80: ********************* 17183 |
|||
-1.70: ************************* 20387 |
|||
-1.60: ****************************** 24049 |
|||
-1.50: ********************************** 27555 |
|||
-1.40: **************************************** 32153 |
|||
-1.30: ********************************************* 36707 |
|||
-1.20: *************************************************** 40921 |
|||
-1.10: ********************************************************* 45928 |
|||
-1.00: *************************************************************** 50707 |
|||
-0.90: ********************************************************************* 55697 |
|||
-0.80: *************************************************************************** 60377 |
|||
-0.70: ******************************************************************************** 64358 |
|||
-0.60: ************************************************************************************ 67928 |
|||
-0.50: ***************************************************************************************** 71911 |
|||
-0.40: ********************************************************************************************* 75054 |
|||
-0.30: ************************************************************************************************ 77073 |
|||
-0.20: ************************************************************************************************** 78768 |
|||
-0.10: *************************************************************************************************** 79732 |
|||
0.00: **************************************************************************************************** 79952 |
|||
0.10: *************************************************************************************************** 79412 |
|||
0.20: ************************************************************************************************ 77511 |
|||
0.30: ********************************************************************************************* 74487 |
|||
0.40: ****************************************************************************************** 72250 |
|||
0.50: ************************************************************************************** 68789 |
|||
0.60: ******************************************************************************** 64408 |
|||
0.70: *************************************************************************** 60122 |
|||
0.80: ********************************************************************* 55619 |
|||
0.90: *************************************************************** 50869 |
|||
1.00: ********************************************************* 45883 |
|||
1.10: **************************************************** 41586 |
|||
1.20: ********************************************** 37145 |
|||
1.30: *************************************** 31715 |
|||
1.40: ********************************** 27779 |
|||
1.50: ****************************** 24270 |
|||
1.60: ************************* 20516 |
|||
1.70: ********************* 17221 |
|||
1.80: ***************** 14344 |
|||
1.90: ************** 11789 |
|||
2.00: ************ 9796 |
|||
2.10: ********* 7922 |
|||
2.20: ******* 6331 |
|||
2.30: ****** 5138 |
|||
2.40: ***** 4044 |
|||
2.50: *** 3065 |
|||
2.60: ** 2397 |
|||
2.70: ** 1846 |
|||
2.80: * 1462 |
|||
2.90: * 1001 |
|||
3.00: 765 |
|||
3.10: 587 |
|||
3.20: 393 |
|||
3.30: 299 |
|||
3.40: 197 |
|||
3.50: 132 |
|||
3.60: 100 |
|||
3.70: 74 |
|||
3.80: 59 |
|||
3.90: 32 |
|||
4.00: 29 |
|||
4.10: 12 |
|||
4.20: 15 |
|||
4.30: 6 |
|||
4.40: 3 |
|||
4.50: 4 |
|||
4.60: 3 |
|||
4.70: 2 |
|||
4.80: 1 |
|||
</pre> |
</pre> |
||
Line 523: | Line 890: | ||
{{works with|Fortran|95 and later}} |
{{works with|Fortran|95 and later}} |
||
Using the Marsaglia polar method |
Using the Marsaglia polar method |
||
< |
<syntaxhighlight lang="fortran">program Normal_Distribution |
||
implicit none |
implicit none |
||
Line 571: | Line 938: | ||
end do |
end do |
||
end program</ |
end program</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 646: | Line 1,013: | ||
=={{header|FreeBASIC}}== |
=={{header|FreeBASIC}}== |
||
< |
<syntaxhighlight lang="freebasic">' FB 1.05.0 Win64 |
||
Const pi As Double = 3.141592653589793 |
Const pi As Double = 3.141592653589793 |
||
Line 732: | Line 1,099: | ||
Print |
Print |
||
Print "Press any key to quit" |
Print "Press any key to quit" |
||
Sleep</ |
Sleep</syntaxhighlight> |
||
Sample output: |
Sample output: |
||
{{out}} |
{{out}} |
||
Line 807: | Line 1,174: | ||
=={{header|Go}}== |
=={{header|Go}}== |
||
Box-Muller method shown here. Go has a normally distributed random function in the standard library, as shown in the Go [[Random numbers]] solution. It uses the ziggurat method. |
Box-Muller method shown here. Go has a normally distributed random function in the standard library, as shown in the Go [[Random numbers]] solution. It uses the ziggurat method. |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 850: | Line 1,217: | ||
fmt.Println(strings.Repeat("*", p/scale)) |
fmt.Println(strings.Repeat("*", p/scale)) |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
Output: |
Output: |
||
<pre> |
<pre> |
||
Line 869: | Line 1,236: | ||
=={{header|Haskell}}== |
=={{header|Haskell}}== |
||
< |
<syntaxhighlight lang="haskell">import Data.Map (Map, empty, insert, findWithDefault, toList) |
||
import Data.Maybe (fromMaybe) |
import Data.Maybe (fromMaybe) |
||
import Text.Printf (printf) |
import Text.Printf (printf) |
||
Line 929: | Line 1,296: | ||
main = do |
main = do |
||
runTest 1000 |
runTest 1000 |
||
runTest 2000000</ |
runTest 2000000</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,104: | Line 1,471: | ||
=={{header|J}}== |
=={{header|J}}== |
||
'''Solution''' |
'''Solution''' |
||
< |
<syntaxhighlight lang="j">runif01=: ?@$ 0: NB. random uniform number generator |
||
rnorm01=. (2 o. 2p1 * runif01) * [: %: _2 * ^.@runif01 NB. random normal number generator (Box-Muller) |
rnorm01=. (2 o. 2p1 * runif01) * [: %: _2 * ^.@runif01 NB. random normal number generator (Box-Muller) |
||
mean=: +/ % # NB. mean |
mean=: +/ % # NB. mean |
||
stddev=: (<:@# %~ +/)&.:*:@(- mean) NB. standard deviation |
stddev=: (<:@# %~ +/)&.:*:@(- mean) NB. standard deviation |
||
histogram=: <:@(#/.~)@(i.@#@[ , I.)</ |
histogram=: <:@(#/.~)@(i.@#@[ , I.)</syntaxhighlight> |
||
'''Example Usage''' |
'''Example Usage''' |
||
< |
<syntaxhighlight lang="j"> DataSet=: rnorm01 1e5 |
||
(mean , stddev) DataSet |
(mean , stddev) DataSet |
||
0.000781667 1.00154 |
0.000781667 1.00154 |
||
require 'plot' |
require 'plot' |
||
plot (5 %~ i: 25) ([;histogram) DataSet</ |
plot (5 %~ i: 25) ([;histogram) DataSet</syntaxhighlight> |
||
=={{header|Java}}== |
=={{header|Java}}== |
||
{{trans|D}} |
{{trans|D}} |
||
{{works with|Java|8}} |
{{works with|Java|8}} |
||
< |
<syntaxhighlight lang="java">import static java.lang.Math.*; |
||
import static java.util.Arrays.stream; |
import static java.util.Arrays.stream; |
||
import java.util.Locale; |
import java.util.Locale; |
||
Line 1,197: | Line 1,564: | ||
.toArray()); |
.toArray()); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
<pre>Mean: -0.001870, SD: 0.500539 |
<pre>Mean: -0.001870, SD: 0.500539 |
||
0.0: ** |
0.0: ** |
||
Line 1,209: | Line 1,576: | ||
0.8: ******* |
0.8: ******* |
||
0.9: **</pre> |
0.9: **</pre> |
||
=={{header|jq}}== |
|||
'''Adapted from [[#Wren|Wren]]''' |
|||
{{works with|jq}} |
|||
'''Works with gojq, the Go implementation of jq''' (*) |
|||
Since jq does not have a built-in PRNG, this entry uses an external source |
|||
for entropy. For the sake of illustration, we will use /dev/urandom |
|||
as follows: |
|||
cat /dev/urandom | tr -cd '0-9' | fold -w 10 | |
|||
jq -nRr -f normal-distribution.jq |
|||
To save space, the function that generates the sample does not retain the observations, and for |
|||
simplicity, computes the sum of squared observations on the |
|||
assumption that overflow will not be an an issue, which is |
|||
reasonable as jq arithmetic uses IEEE 754 64-bit numbers. |
|||
(*) gojq requires an enormous amount of memory to complete the task for N=100,000, |
|||
and takes about 20 times longer. |
|||
'''Preliminaries''' |
|||
<syntaxhighlight lang="jq"># Pretty print a number to facilitate alignment of the decimal point. |
|||
# Input: a number without an exponent |
|||
# Output: a string holding the reformatted number so that there are at least `left` characters |
|||
# to the left of the decimal point, and exactly `right` characters to its right. |
|||
# Spaces are used for padding on the left, and zeros for padding on the right. |
|||
# No left-truncation occurs, so `left` can be specified as 0 to prevent left-padding. |
|||
def pp(left; right): |
|||
def lpad: tostring | (left - length) as $l | (" " * $l)[:$l] + .; |
|||
def rpad: |
|||
if (right > length) then . + ((right - length) * "0") |
|||
else .[:right] |
|||
end; |
|||
tostring as $s |
|||
| $s |
|||
| index(".") as $ix |
|||
| ((if $ix then $s[0:$ix] else $s end) | lpad) + "." + |
|||
(if $ix then $s[$ix+1:] | .[:right] else "" end | rpad) ; |
|||
def sigma( stream ): reduce stream as $x (0; . + $x); |
|||
# Input: {n, sum, ss} |
|||
# Output: augmented object with {mean, variance} |
|||
def sample_mean_and_variance: |
|||
.mean = (.sum/.n) |
|||
| .variance = ((.ss / .n) - .mean*.mean);</syntaxhighlight> |
|||
'''The Task''' |
|||
<syntaxhighlight lang="jq"># Task parameters |
|||
def parameters: { |
|||
N: 100000, |
|||
NUM_BINS: 12, |
|||
HIST_CHAR: "■", |
|||
HIST_CHAR_ALT: "-", |
|||
HIST_CHAR_SIZE: null, # null means compute dynamically |
|||
binSize: 0.1, |
|||
mu: 0.5, |
|||
sigma: 0.25 } |
|||
| .bins = [range(0; .NUM_BINS) | 0] ; |
|||
# input: an array of two iid rvs on [0,1] |
|||
# output: [z0, z1] as per the Box-Muller method -- see |
|||
# https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform |
|||
def normal(mu; sigma): |
|||
def pi: (1|atan) * 4; |
|||
. as [$u1, $u2] |
|||
| pi as $pi |
|||
| (sigma * ((-2 * ($u1|log))|sqrt)) as $mag |
|||
| [ $mag * ((2 * $pi * $u2)|cos) + mu, |
|||
$mag * ((2 * $pi * $u2)|sin) + mu ] ; |
|||
# Generate a random sample as specified by ., the task object (see `parameters`). |
|||
# Output: updated task object with sample statistics and .bins for creating a histogram. |
|||
# Each call to `input` should yield a string of random decimal digits |
|||
# such that the ensemble of ("0." + input | tonumber) can be considered to be iid rv on [0,1]. |
|||
def generate: |
|||
# uniformly distributed random variable on [0,1]: |
|||
def udrv: "0." + input | tonumber; |
|||
# Maybe compute the bucket size: |
|||
(.HIST_CHAR_SIZE = (.HIST_CHAR_SIZE // (.N / (.NUM_BINS * 20) | ceil))) as $p |
|||
| reduce range(0; $p.N/2) as $i ($p; |
|||
([udrv, udrv] | normal($p.mu; $p.sigma)) as $rns |
|||
| reduce (0,1) as $j (.; |
|||
$rns[$j] as $rn |
|||
| .n += 1 |
|||
| .sum += $rn |
|||
| .ss += ($rn*$rn) |
|||
| (if $rn < 0 then 0 |
|||
elif $rn >= 1 then ($p.NUM_BINS - 1) |
|||
else ($rn/.binSize)|floor + 1 |
|||
end ) as $bn |
|||
| .bins[$bn] += 1 |
|||
# to retain the observations: .samples[$i*2 + $j] = $rn |
|||
)) ; |
|||
# Input: an object with |
|||
# {NUM_BINS, HIST_CHAR_SIZE, HIST_CHAR, HIST_CHAR_ALT, binSize, bins} |
|||
# Output: a stream of strings |
|||
def histogram: |
|||
def tidy: pp(2;1); |
|||
range(0; .NUM_BINS) as $i |
|||
| ((.bins[$i] / .HIST_CHAR_SIZE)|floor) as $bs |
|||
| (if $i == 0 or $i == .NUM_BINS -1 |
|||
then .HIST_CHAR_ALT else .HIST_CHAR end) as $char |
|||
| (if $bs == 0 then "" else $char * $bs end) as $hist |
|||
| if $i == 0 |
|||
then " -∞ ..< 0.0 \($hist)" # .bins[0] |
|||
elif ($i < .NUM_BINS - 1) |
|||
then "\(.binSize * ($i-1) | tidy) ..<\(.binSize * $i|tidy) \($hist)" # .bins[$i]] |
|||
else " 1.0 .. +∞ \($hist)" # .bins[.NUM_BINS - 1] |
|||
end; |
|||
def task: |
|||
parameters |
|||
| generate |
|||
| sample_mean_and_variance |
|||
| (if .HIST_CHAR_SIZE == 1 then "" else "s" end) as $plural |
|||
| "Summary statistics for \(.N) observations from N(\(.mu), \(.sigma)):", |
|||
" mean: \(.mean | pp(2;4))", |
|||
" variance: \(.variance | pp(2;4))", |
|||
" unadjusted stddev: \(.variance | sqrt | pp(2;4))", |
|||
" Range Number of observations (each \(.HIST_CHAR) represents \(.HIST_CHAR_SIZE) observation\($plural))", |
|||
histogram ; |
|||
task</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Summary statistics for 100000 observations from N(0.5, 0.25): |
|||
mean: 0.5001 |
|||
variance: 0.0622 |
|||
unadjusted stddev: 0.2495 |
|||
Range Number of observations (each ■ represents 417 observations) |
|||
-∞ ..< 0.0 ----- |
|||
0.0 ..< 0.1 ■■■■■■■■ |
|||
0.1 ..< 0.2 ■■■■■■■■■■■■■■ |
|||
0.2 ..< 0.3 ■■■■■■■■■■■■■■■■■■■■■■■ |
|||
0.3 ..< 0.4 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ |
|||
0.4 ..< 0.5 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ |
|||
0.5 ..< 0.6 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ |
|||
0.6 ..< 0.7 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ |
|||
0.7 ..< 0.8 ■■■■■■■■■■■■■■■■■■■■■■■ |
|||
0.8 ..< 0.9 ■■■■■■■■■■■■■■ |
|||
0.9 ..< 1.0 ■■■■■■■■ |
|||
1.0 .. +∞ ------ |
|||
</pre> |
|||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
Julia has the builtin package "Distributions" to generate random numbers from a standard distribution (Normal, Chisq etc.). |
Julia has the builtin package "Distributions" to generate random numbers from a standard distribution (Normal, Chisq etc.). |
||
< |
<syntaxhighlight lang="julia">using Printf, Distributions, Gadfly |
||
data = rand(Normal(0, 1), 1000) |
data = rand(Normal(0, 1), 1000) |
||
Line 1,219: | Line 1,731: | ||
@printf("range = (%2.2f, %2.2f\n)", minimum(data), maximum(data)) |
@printf("range = (%2.2f, %2.2f\n)", minimum(data), maximum(data)) |
||
h = plot(x=data, Geom.histogram) |
h = plot(x=data, Geom.histogram) |
||
draw(PNG("norm_hist.png", 10cm, 10cm), h)</ |
draw(PNG("norm_hist.png", 10cm, 10cm), h)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,228: | Line 1,740: | ||
=={{header|Kotlin}}== |
=={{header|Kotlin}}== |
||
{{trans|FreeBASIC}} |
{{trans|FreeBASIC}} |
||
< |
<syntaxhighlight lang="scala">// version 1.1.2 |
||
val rand = java.util.Random() |
val rand = java.util.Random() |
||
Line 1,284: | Line 1,796: | ||
val sampleSizes = intArrayOf(100, 1_000, 10_000, 100_000) |
val sampleSizes = intArrayOf(100, 1_000, 10_000, 100_000) |
||
for (sampleSize in sampleSizes) normalStats(sampleSize) |
for (sampleSize in sampleSizes) normalStats(sampleSize) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,358: | Line 1,870: | ||
=={{header|Lasso}}== |
=={{header|Lasso}}== |
||
< |
<syntaxhighlight lang="lasso">define stat1(a) => { |
||
if(#a->size) => { |
if(#a->size) => { |
||
local(mean = (with n in #a sum #n) / #a->size) |
local(mean = (with n in #a sum #n) / #a->size) |
||
Line 1,419: | Line 1,931: | ||
histogram(#n) |
histogram(#n) |
||
'\r\r' |
'\r\r' |
||
^}</ |
^}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,474: | Line 1,986: | ||
=={{header|Liberty BASIC}}== |
=={{header|Liberty BASIC}}== |
||
Uses LB Statistics/Basic |
Uses LB Statistics/Basic |
||
< |
<syntaxhighlight lang="lb">call sample 100000 |
||
end |
end |
||
Line 1,525: | Line 2,037: | ||
v =rnd( 1) |
v =rnd( 1) |
||
normalDist =( -2 *log( u))^0.5 *cos( 2 *3.14159265 *v) |
normalDist =( -2 *log( u))^0.5 *cos( 2 *3.14159265 *v) |
||
end function</ |
end function</syntaxhighlight> |
||
100000 data terms used. |
100000 data terms used. |
||
Largest term was 4.12950792 & smallest was -4.37934139 |
Largest term was 4.12950792 & smallest was -4.37934139 |
||
Line 1,571: | Line 2,083: | ||
=={{header|Lua}}== |
=={{header|Lua}}== |
||
Lua provides math.random() to generate uniformly distributed random numbers. The function gaussian() shown here uses math.random() to generate normally distributed random numbers with given mean and variance. |
Lua provides math.random() to generate uniformly distributed random numbers. The function gaussian() shown here uses math.random() to generate normally distributed random numbers with given mean and variance. |
||
< |
<syntaxhighlight lang="lua">function gaussian (mean, variance) |
||
return math.sqrt(-2 * variance * math.log(math.random())) * |
return math.sqrt(-2 * variance * math.log(math.random())) * |
||
math.cos(2 * math.pi * math.random()) + mean |
math.cos(2 * math.pi * math.random()) + mean |
||
Line 1,616: | Line 2,128: | ||
print("Mean:", mean(t) .. ", expected " .. average) |
print("Mean:", mean(t) .. ", expected " .. average) |
||
print("StdDev:", std(t) .. ", expected " .. math.sqrt(variance) .. "\n") |
print("StdDev:", std(t) .. ", expected " .. math.sqrt(variance) .. "\n") |
||
showHistogram(t)</ |
showHistogram(t)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Mean: 50.008328894275, expected 50 |
<pre>Mean: 50.008328894275, expected 50 |
||
Line 1,644: | Line 2,156: | ||
=={{header|Maple}}== |
=={{header|Maple}}== |
||
Maple has a built-in for sampling directly from [http://www.maplesoft.com/support/help/Maple/view.aspx?path=Statistics/Distributions/Normal Normal] distributions: |
Maple has a built-in for sampling directly from [http://www.maplesoft.com/support/help/Maple/view.aspx?path=Statistics/Distributions/Normal Normal] distributions: |
||
< |
<syntaxhighlight lang="maple">with(Statistics): |
||
n := 100000: |
n := 100000: |
||
X := Sample( Normal(0,1), n ); |
X := Sample( Normal(0,1), n ); |
||
Mean( X ); |
Mean( X ); |
||
StandardDeviation( X ); |
StandardDeviation( X ); |
||
Histogram( X );</ |
Histogram( X );</syntaxhighlight> |
||
=={{header|Mathematica}}== |
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
||
< |
<syntaxhighlight lang="mathematica">x:= RandomReal[1] |
||
SampleNormal[n_] := (Print[#//Length, " numbers, Mean : ", #//Mean, ", StandardDeviation : ", #//StandardDeviation]; |
SampleNormal[n_] := (Print[#//Length, " numbers, Mean : ", #//Mean, ", StandardDeviation : ", #//StandardDeviation]; |
||
Histogram[#, BarOrigin -> Left,Axes -> False])& [(Table[(-2*Log[x])^0.5*Cos[2*Pi*x], {n} ]] |
Histogram[#, BarOrigin -> Left,Axes -> False])& [(Table[(-2*Log[x])^0.5*Cos[2*Pi*x], {n} ]] |
||
Invocation: |
Invocation: |
||
SampleNormal[ 10000 ] |
SampleNormal[ 10000 ] |
||
->10000 numbers, Mean : -0.0122308, StandardDeviation : 1.00646 |
->10000 numbers, Mean : -0.0122308, StandardDeviation : 1.00646 |
||
</syntaxhighlight> |
|||
</lang> |
|||
[[File:Mma_NormalDistribution.png]] |
[[File:Mma_NormalDistribution.png]] |
||
=={{header|MATLAB}} / {{header|Octave}}== |
=={{header|MATLAB}} / {{header|Octave}}== |
||
< |
<syntaxhighlight lang="matlab"> N = 100000; |
||
x = randn(N,1); |
x = randn(N,1); |
||
mean(x) |
mean(x) |
||
std(x) |
std(x) |
||
[nn,xx] = hist(x,100); |
[nn,xx] = hist(x,100); |
||
bar(xx,nn);</ |
bar(xx,nn);</syntaxhighlight> |
||
=={{header|Nim}}== |
|||
In module “random”, Nim provides two procedures named <code>gauss</code> to generate random values following normal distribution and following Gauss distribution with given mean and standard deviation. |
|||
Here is a way to generate random values following normal distribution from random values following uniform distribution. It uses the Basic form of the Box-Muller transform. |
|||
<syntaxhighlight lang="nim">import math, random, sequtils, stats, strformat, strutils |
|||
proc drawHistogram(ns: seq[float]) = |
|||
# Distribute values in bins. |
|||
const NBins = 50 |
|||
var minval = min(ns) |
|||
var maxval = max(ns) |
|||
var h = newSeq[int](NBins + 1) |
|||
for n in ns: |
|||
let pos = ((n - minval) * NBins / (maxval - minval)).toInt |
|||
inc h[pos] |
|||
# Eliminate extremes values. |
|||
const MaxWidth = 50 |
|||
let mx = max(h) |
|||
var first = 0 |
|||
while (h[first] / mx * MaxWidth).toInt == 0: inc first |
|||
var last = h.high |
|||
while (h[last] / mx * MaxWidth).toInt == 0: dec last |
|||
# Draw the histogram. |
|||
echo "" |
|||
for n in first..last: |
|||
echo repeat('+', (h[n] / mx * MaxWidth).toInt) |
|||
echo "" |
|||
const N = 100_000 |
|||
randomize() |
|||
let u1, u2 = newSeqWith(N, rand(1.0)) |
|||
var z = newSeq[float](N) |
|||
for i in 0..<N: |
|||
z[i] = sqrt(-2 * ln(u1[i])) * cos(2 * PI * u2[i]) |
|||
echo &"μ = {z.mean:.12f} σ = {z.standardDeviation:.12f}" |
|||
z.drawHistogram()</syntaxhighlight> |
|||
{{out}} |
|||
<pre>μ = -0.001105836229 σ = 0.999906544722 |
|||
+ |
|||
+ |
|||
++ |
|||
+++ |
|||
+++++ |
|||
+++++++ |
|||
+++++++++ |
|||
++++++++++++ |
|||
++++++++++++++++ |
|||
+++++++++++++++++++++ |
|||
++++++++++++++++++++++++++ |
|||
+++++++++++++++++++++++++++++++ |
|||
+++++++++++++++++++++++++++++++++++++ |
|||
++++++++++++++++++++++++++++++++++++++++ |
|||
+++++++++++++++++++++++++++++++++++++++++++++ |
|||
++++++++++++++++++++++++++++++++++++++++++++++++ |
|||
++++++++++++++++++++++++++++++++++++++++++++++++++ |
|||
++++++++++++++++++++++++++++++++++++++++++++++++++ |
|||
+++++++++++++++++++++++++++++++++++++++++++++++ |
|||
+++++++++++++++++++++++++++++++++++++++++++ |
|||
+++++++++++++++++++++++++++++++++++++++ |
|||
++++++++++++++++++++++++++++++++++++ |
|||
+++++++++++++++++++++++++++++++ |
|||
+++++++++++++++++++++++++ |
|||
++++++++++++++++++++ |
|||
++++++++++++++++ |
|||
++++++++++++ |
|||
+++++++++ |
|||
++++++ |
|||
+++++ |
|||
+++ |
|||
++ |
|||
+ |
|||
+ |
|||
</pre> |
|||
=={{header|PARI/GP}}== |
=={{header|PARI/GP}}== |
||
{{works with|PARI/GP|2.4.3 and above}} |
{{works with|PARI/GP|2.4.3 and above}} |
||
< |
<syntaxhighlight lang="parigp">rnormal()={ |
||
my(u1=random(1.),u2=random(1.); |
my(u1=random(1.),u2=random(1.)); |
||
sqrt(-2*log(u1))*cos(2*Pi* |
sqrt(-2*log(u1))*cos(2*Pi*u2) |
||
\\ Could easily be extended with a second normal at very little cost. |
\\ Could easily be extended with a second normal at very little cost. |
||
}; |
}; |
||
Line 1,694: | Line 2,290: | ||
for(i=1,#h,for(j=1,h[i]\sz,print1("#"));print()); |
for(i=1,#h,for(j=1,h[i]\sz,print1("#"));print()); |
||
}; |
}; |
||
show(10^4)</ |
show(10^4)</syntaxhighlight> |
||
For versions before 2.4.3, define |
For versions before 2.4.3, define |
||
< |
<syntaxhighlight lang="parigp">rreal()={ |
||
my(pr=32*ceil(default(realprecision)*log(10)/log(4294967296))); \\ Current precision |
my(pr=32*ceil(default(realprecision)*log(10)/log(4294967296))); \\ Current precision |
||
random(2^pr)*1.>>pr |
random(2^pr)*1.>>pr |
||
};</ |
};</syntaxhighlight> |
||
and use <code>rreal()</code> in place of <code>random(1.)</code>. |
and use <code>rreal()</code> in place of <code>random(1.)</code>. |
||
A PARI implementation: |
A PARI implementation: |
||
< |
<syntaxhighlight lang="c">GEN |
||
rnormal(long prec) |
rnormal(long prec) |
||
{ |
{ |
||
Line 1,717: | Line 2,313: | ||
ret = gerepileupto(ltop, ret); |
ret = gerepileupto(ltop, ret); |
||
return ret; |
return ret; |
||
}</ |
}</syntaxhighlight> |
||
Use <code>mpsincos</code> and caching to generate two values at nearly the same cost. |
Use <code>mpsincos</code> and caching to generate two values at nearly the same cost. |
||
Line 1,727: | Line 2,323: | ||
From [http://www.freepascal.org/docs-html/rtl/math/randg.html Free Pascal Docs unit math] |
From [http://www.freepascal.org/docs-html/rtl/math/randg.html Free Pascal Docs unit math] |
||
< |
<syntaxhighlight lang="pascal">Program Example40; |
||
{$IFDEF FPC} |
{$IFDEF FPC} |
||
{$MOde objFPC} |
{$MOde objFPC} |
||
Line 1,858: | Line 2,454: | ||
mySol := getSol(@rnorm,cMean,cStdDiv,1000000); |
mySol := getSol(@rnorm,cMean,cStdDiv,1000000); |
||
Histo(mySol,cHistCnt,cColLen); |
Histo(mySol,cHistCnt,cColLen); |
||
end.</ |
end.</syntaxhighlight> |
||
{{out}}<pre> |
{{out}}<pre> |
||
function randg |
function randg |
||
Line 1,919: | Line 2,515: | ||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
{{trans|Raku}} |
{{trans|Raku}} |
||
< |
<syntaxhighlight lang="perl">use constant pi => 3.14159265; |
||
use List::Util qw(sum reduce min max); |
use List::Util qw(sum reduce min max); |
||
Line 1,950: | Line 2,546: | ||
print "$i\t$t1$t2\n"; |
print "$i\t$t1$t2\n"; |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre style="height:35ex">32 ⎸ |
<pre style="height:35ex">32 ⎸ |
||
Line 1,990: | Line 2,586: | ||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
{{trans|Liberty_BASIC}} |
{{trans|Liberty_BASIC}} |
||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<lang Phix>procedure sample(integer n) |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
-- show mean, standard deviation. Find max, min. |
|||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">sample</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
sequence dat = repeat(0,n) |
|||
<span style="color: #000080;font-style:italic;">-- show mean, standard deviation. Find max, min.</span> |
|||
for i=1 to n do |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">dat</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;">n</span><span style="color: #0000FF;">)</span> |
|||
dat[i] = sqrt(-2*log(rnd()))*cos(2*PI*rnd()) |
|||
<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: #000000;">n</span> <span style="color: #008080;">do</span> |
|||
end for |
|||
<span style="color: #000000;">dat</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: #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: #7060A8;">cos</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">*</span><span style="color: #7060A8;">rnd</span><span style="color: #0000FF;">())</span> |
|||
printf(1,"%d data terms used.\n",{n}) |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%d data terms used.\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">})</span> |
|||
atom mean = sum(dat)/n, |
|||
mx = max(dat), |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">mean</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dat</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> |
|||
mn = min(dat), |
|||
<span style="color: #000000;">mx</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">max</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dat</span><span style="color: #0000FF;">),</span> |
|||
range = mx-mn |
|||
<span style="color: #000000;">mn</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">min</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dat</span><span style="color: #0000FF;">),</span> |
|||
printf(1,"Largest term was %g & smallest was %g\n",{mx,mn}) |
|||
<span style="color: #000000;">range</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mx</span><span style="color: #0000FF;">-</span><span style="color: #000000;">mn</span> |
|||
printf(1,"Mean = %g\n",{mean}) |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Largest term was %g & smallest was %g\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">mx</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mn</span><span style="color: #0000FF;">})</span> |
|||
printf(1,"Stddev = %g\n",sqrt(sum(sq_mul(dat,dat))/n-mean*mean)) |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Mean = %g\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">mean</span><span style="color: #0000FF;">})</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Stddev = %g\n"</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dat</span><span style="color: #0000FF;">,</span><span style="color: #000000;">dat</span><span style="color: #0000FF;">))/</span><span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">mean</span><span style="color: #0000FF;">*</span><span style="color: #000000;">mean</span><span style="color: #0000FF;">))</span> |
|||
-- show histogram |
|||
integer nBins = 50 |
|||
<span style="color: #000080;font-style:italic;">-- show histogram</span> |
|||
sequence bins = repeat(0,nBins+1) |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">nBins</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">50</span> |
|||
for i=1 to n do |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">bins</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;">nBins</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
bins[floor((dat[i]-mn)/range*nBins)+1] += 1 |
|||
<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: #000000;">n</span> <span style="color: #008080;">do</span> |
|||
end for |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">bdx</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">((</span><span style="color: #000000;">dat</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]-</span><span style="color: #000000;">mn</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">range</span><span style="color: #0000FF;">*</span><span style="color: #000000;">nBins</span><span style="color: #0000FF;">)+</span><span style="color: #000000;">1</span> |
|||
for b=1 to nBins do |
|||
<span style="color: #000000;">bins</span><span style="color: #0000FF;">[</span><span style="color: #000000;">bdx</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span> |
|||
puts(1,repeat('#',floor(nBins*bins[b]/n*30))&"\n") |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
end for |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">nBins</span> <span style="color: #008080;">do</span> |
|||
end procedure |
|||
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #008000;">'#'</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nBins</span><span style="color: #0000FF;">*</span><span style="color: #000000;">bins</span><span style="color: #0000FF;">[</span><span style="color: #000000;">b</span><span style="color: #0000FF;">]/</span><span style="color: #000000;">n</span><span style="color: #0000FF;">*</span><span style="color: #000000;">30</span><span style="color: #0000FF;">))&</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
sample(100000)</lang> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
|||
<span style="color: #000000;">sample</span><span style="color: #0000FF;">(</span><span style="color: #000000;">100000</span><span style="color: #0000FF;">)</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{Out}} |
{{Out}} |
||
<pre> |
<pre> |
||
Line 2,062: | Line 2,662: | ||
</pre> |
</pre> |
||
{{trans|Lua}} |
{{trans|Lua}} |
||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<lang Phix>function gaussian(atom mean, atom variance) |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
return sqrt(-2 * variance * log(rnd())) * |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">gaussian</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">mean</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">variance</span><span style="color: #0000FF;">)</span> |
|||
cos(2 * variance * PI * rnd()) + mean |
|||
<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: #000000;">variance</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> |
|||
end function |
|||
<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;">variance</span> <span style="color: #0000FF;">*</span> <span style="color: #004600;">PI</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: #000000;">mean</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
function mean(sequence t) |
|||
return sum(t)/length(t) |
|||
end function |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">mean</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> |
|||
function std(sequence t) |
|||
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)/</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> |
|||
atom squares = 0, |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
avg = mean(t) |
|||
for i=1 to length(t) do |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">std</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> |
|||
squares += power(avg-t[i],2) |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">squares</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> |
|||
end for |
|||
<span style="color: #000000;">avg</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mean</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> |
|||
atom variance = squares/length(t) |
|||
<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;">t</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
return sqrt(variance) |
|||
<span style="color: #000000;">squares</span> <span style="color: #0000FF;">+=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">avg</span><span style="color: #0000FF;">-</span><span style="color: #000000;">t</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span> |
|||
end function |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">variance</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">squares</span><span style="color: #0000FF;">/</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> |
|||
procedure showHistogram(sequence t) |
|||
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">variance</span><span style="color: #0000FF;">)</span> |
|||
for i=ceil(min(t)) to floor(max(t)) do |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
integer n = 0 |
|||
for k=1 to length(t) do |
|||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">showHistogram</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> |
|||
n += ceil(t[k]-0.5)=i |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #7060A8;">ceil</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">min</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">))</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">max</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">))</span> <span style="color: #008080;">do</span> |
|||
end for |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span> |
|||
integer l = floor(n/length(t)*200) |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">k</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;">t</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
printf(1,"%d %s %d\n",{i,repeat('=',l),n}) |
|||
<span style="color: #000000;">n</span> <span style="color: #0000FF;">+=</span> <span style="color: #7060A8;">ceil</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]-</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">i</span> |
|||
end for |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
end procedure |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">l</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">/</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">200</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%d %s %d\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #008000;">'='</span><span style="color: #0000FF;">,</span><span style="color: #000000;">l</span><span style="color: #0000FF;">),</span><span style="color: #000000;">n</span><span style="color: #0000FF;">})</span> |
|||
sequence t = repeat(0,100000) |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
integer avg = 50, variance = 10 |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
|||
for i=1 to length(t) do |
|||
t[i] = gaussian(avg, variance) |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">t</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;">100000</span><span style="color: #0000FF;">)</span> |
|||
end for |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">avg</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">50</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">variance</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">10</span> |
|||
printf(1,"Mean: %g, expected %g\n",{mean(t),avg}) |
|||
<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;">t</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
printf(1,"StdDev: %g, expected %g\n",{std(t),sqrt(variance)}) |
|||
<span style="color: #000000;">t</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;">gaussian</span><span style="color: #0000FF;">(</span><span style="color: #000000;">avg</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">variance</span><span style="color: #0000FF;">)</span> |
|||
showHistogram(t)</lang> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Mean: %g, expected %g\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">mean</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">),</span><span style="color: #000000;">avg</span><span style="color: #0000FF;">})</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"StdDev: %g, expected %g\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">std</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">),</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">variance</span><span style="color: #0000FF;">)})</span> |
|||
<span style="color: #000000;">showHistogram</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{Out}} |
{{Out}} |
||
<pre> |
<pre> |
||
Line 2,134: | Line 2,737: | ||
=={{header|PureBasic}}== |
=={{header|PureBasic}}== |
||
< |
<syntaxhighlight lang="purebasic">Procedure.f randomf(resolution = 2147483647) |
||
ProcedureReturn Random(resolution) / resolution |
ProcedureReturn Random(resolution) / resolution |
||
EndProcedure |
EndProcedure |
||
Line 2,198: | Line 2,801: | ||
Print(#CRLF$ + #CRLF$ + "Press ENTER to exit"): Input() |
Print(#CRLF$ + #CRLF$ + "Press ENTER to exit"): Input() |
||
CloseConsole() |
CloseConsole() |
||
EndIf</ |
EndIf</syntaxhighlight> |
||
Sample output: |
Sample output: |
||
<pre>100000 data terms used. |
<pre>100000 data terms used. |
||
Line 2,237: | Line 2,840: | ||
=={{header|Python}}== |
=={{header|Python}}== |
||
This uses the external [http://matplotlib.org/ matplotlib] package as well as the built-in standardlib function [http://docs.python.org/2/library/random.html?highlight=gauss#random.gauss random.gauss]. |
This uses the external [http://matplotlib.org/ matplotlib] package as well as the built-in standardlib function [http://docs.python.org/2/library/random.html?highlight=gauss#random.gauss random.gauss]. |
||
< |
<syntaxhighlight lang="python">from __future__ import division |
||
import matplotlib.pyplot as plt |
import matplotlib.pyplot as plt |
||
import random |
import random |
||
Line 2,251: | Line 2,854: | ||
% (mn, sd, max(data), min(data), size)) |
% (mn, sd, max(data), min(data), size)) |
||
plt.hist(data,bins=50)</ |
plt.hist(data,bins=50)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 2,259: | Line 2,862: | ||
=={{header|R}}== |
=={{header|R}}== |
||
R can generate random normal distributed numbers using the [https://stat.ethz.ch/R-manual/R-devel/library/stats/html/Normal.html rnorm] command: |
|||
Generate normal random numbers from uniform random numbers, using the [https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform Box-Muller transform]. Both x and y are normally distributed, and independent. |
|||
<lang r>n <- 100000 |
|||
<syntaxhighlight lang="r">n <- 100000 |
|||
u <- sqrt(-2*log(runif(n))) |
|||
v <- 2*pi*runif(n) |
|||
x <- u*cos(v) |
|||
y <- v*sin(v) |
|||
hist(x)</syntaxhighlight> |
|||
R can generate random normal distributed numbers using the [https://stat.ethz.ch/R-manual/R-devel/library/stats/html/Normal.html rnorm] function: |
|||
<syntaxhighlight lang="r">n <- 100000 |
|||
x <- rnorm(n, mean=0, sd=1) |
x <- rnorm(n, mean=0, sd=1) |
||
mean(x) |
mean(x) |
||
sd(x) |
sd(x) |
||
hist(x)</ |
hist(x)</syntaxhighlight> |
||
=={{header|Racket}}== |
=={{header|Racket}}== |
||
Line 2,270: | Line 2,883: | ||
compute statistics and plot a histogram. |
compute statistics and plot a histogram. |
||
[[File:histogram-racket.png|thumb|right]] |
[[File:histogram-racket.png|thumb|right]] |
||
< |
<syntaxhighlight lang="racket"> |
||
#lang racket |
#lang racket |
||
(require math (planet williams/science/histogram-with-graphics)) |
(require math (planet williams/science/histogram-with-graphics)) |
||
Line 2,284: | Line 2,897: | ||
(for ([x data]) (histogram-increment! h x)) |
(for ([x data]) (histogram-increment! h x)) |
||
(histogram-plot h "Normal distribution μ=50 σ=4") |
(histogram-plot h "Normal distribution μ=50 σ=4") |
||
</syntaxhighlight> |
|||
</lang> |
|||
The other part of the task was to produce normal distributed numbers from a unit distribution. |
The other part of the task was to produce normal distributed numbers from a unit distribution. |
||
Line 2,290: | Line 2,903: | ||
version of [http://planet.plt-scheme.org/package-source/schematics/random.plt/1/0/random.ss code] |
version of [http://planet.plt-scheme.org/package-source/schematics/random.plt/1/0/random.ss code] |
||
originally written by Sebastian Egner. |
originally written by Sebastian Egner. |
||
< |
<syntaxhighlight lang="racket"> |
||
#lang racket |
#lang racket |
||
(require math) |
(require math) |
||
Line 2,310: | Line 2,923: | ||
(set! next (* scale v2)) |
(set! next (* scale v2)) |
||
(+ μ (* σ scale v1))]))))))) |
(+ μ (* σ scale v1))]))))))) |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Raku}}== |
=={{header|Raku}}== |
||
(formerly Perl 6) |
(formerly Perl 6) |
||
{{works with|Rakudo|2018.03}} |
{{works with|Rakudo|2018.03}} |
||
<lang |
<syntaxhighlight lang="raku" line>sub normdist ($m, $σ) { |
||
my $r = sqrt -2 * log rand; |
my $r = sqrt -2 * log rand; |
||
my $Θ = τ * rand; |
my $Θ = τ * rand; |
||
Line 2,339: | Line 2,952: | ||
say $i, "\t", '█' x $full, @subbar[$part]; |
say $i, "\t", '█' x $full, @subbar[$part]; |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>"m" => 50.006107405837142e0 |
<pre>"m" => 50.006107405837142e0 |
||
Line 2,382: | Line 2,995: | ||
The REXX language doesn't have any "higher math" BIF functions like SIN, COS, LN, LOG, SQRT, EXP, POW, etc, |
The REXX language doesn't have any "higher math" BIF functions like SIN, COS, LN, LOG, SQRT, EXP, POW, etc, |
||
<br>so we hoi polloi programmers have to roll our own. |
<br>so we hoi polloi programmers have to roll our own. |
||
< |
<syntaxhighlight lang="rexx">/*REXX program generates 10,000 normally distributed numbers (Gaussian distribution).*/ |
||
numeric digits 20 /*use |
numeric digits 20 /*use twenty decimal digs for accuracy.*/ |
||
parse arg n seed . /*obtain optional arguments from the CL*/ |
parse arg n seed . /*obtain optional arguments from the CL*/ |
||
if n=='' | n=="," then n= 10000 /*Not specified? Then use the default.*/ |
if n=='' | n=="," then n= 10000 /*Not specified? Then use the default.*/ |
||
Line 2,393: | Line 3,006: | ||
mn= #.1; mx= mn; noise= n * .0005 /*calculate the noise: 1/20th % of N.*/ |
mn= #.1; mx= mn; noise= n * .0005 /*calculate the noise: 1/20th % of N.*/ |
||
ss= 0 |
ss= 0 |
||
do j=1 for n; _=#.j |
do j=1 for n; _= #.j /*the sum, and the sum of squares. */ |
||
s= s + _; ss= ss + _ * _ /*the sum, and the sum of squares. */ |
|||
mn= min(mn, _); mx= max(mx, _) /*find the minimum and the maximum. */ |
mn= min(mn, _); mx= max(mx, _) /*find the minimum and the maximum. */ |
||
end /*j*/ |
end /*j*/ |
||
!.= 0 |
!.= 0 |
||
say 'number of data points = ' |
say 'number of data points = ' fmt(n ) |
||
say ' minimum = ' |
say ' minimum = ' fmt(mn ) |
||
say ' maximum = ' |
say ' maximum = ' fmt(mx ) |
||
say ' arithmetic mean = ' |
say ' arithmetic mean = ' fmt(s/n) |
||
say ' standard deviation = ' |
say ' standard deviation = ' fmt(sqrt( ss/n - (s/n) **2) ) |
||
?mn= !.1; ?mx= ?mn /*define minimum & maximum value so far*/ |
?mn= !.1; ?mx= ?mn /*define minimum & maximum value so far*/ |
||
parse value scrSize() with sd sw . /*obtain the (true) screen size of term*/ /*◄──not all REXXes have this BIF*/ |
parse value scrSize() with sd sw . /*obtain the (true) screen size of term*/ /*◄──not all REXXes have this BIF*/ |
||
sdE= sd - 4 /*the effective ( |
sdE= sd - 4 /*the effective (usable) screen depth. */ |
||
swE= sw - 1 /* " " " " width.*/ |
swE= sw - 1 /* " " " " width.*/ |
||
$= 1 / max(1, mx-mn) * sdE /*$ is used for scaling depth of histo*/ |
$= 1 / max(1, mx-mn) * sdE /*$ is used for scaling depth of histo*/ |
||
do i=1 for n; |
do i=1 for n; ?= trunc((#.i-mn) *$) /*calculate the relative line. */ |
||
!.?= !.? + 1 |
!.?= !.? + 1 /*bump the counter. */ |
||
?mn= min(?mn, !.?) |
?mn= min(?mn, !.?) /*find the minimum. */ |
||
?mx= max(?mx, !.?) /* " " maximum. */ |
|||
end /*i*/ |
end /*i*/ |
||
f=swE/?mx |
f= swE/?mx /*limit graph to 1 full screen*/ |
||
do h=0 for sdE; |
do h=0 for sdE; _= !.h /*obtain a data point. */ |
||
if _>noise then say copies('─', trunc(_*f) ) /*display a bar of histogram. */ |
if _>noise then say copies('─', trunc(_*f) ) /*display a bar of histogram. */ |
||
end /*h*/ /*[↑] use a hyphen for histo.*/ |
end /*h*/ /*[↑] use a hyphen for histo.*/ |
||
exit /*stick a fork in it, we're all done. */ |
exit /*stick a fork in it, we're all done. */ |
||
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ |
|||
/*──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ |
|||
fmt: parse arg @; return left('', (@>=0) + 2 * datatype(@, 'W'))@ /*prepend a blank if #>=0, add 2 blanks if whole.*/ |
|||
e: |
e: e = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535; return e |
||
pi: |
pi: pi= 3.1415926535897932384626433832795028841971693993751058209749445923078164062862; return pi |
||
r2r: |
r2r: return arg(1) // (pi() * 2) /*normalize the given angle (in radians) to ±2pi.*/ |
||
rand: |
rand: return random(1, 1e5) / 1e5 /*REXX generates uniform random postive integers.*/ |
||
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ |
|||
.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 |
|||
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ |
|||
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; z= 1; _= 1 |
|||
x= x*x; p= z; do k=2 by 2; _= -_ * x / (k*(k-1)); z= z + _; if z=p then leave; p= z; end; return z |
|||
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ |
|||
/*──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ |
|||
sqrt: |
sqrt: procedure; parse arg x; if x=0 then return 0; d= digits(); m.= 9; numeric digits; numeric form; h= d+6 |
||
parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g*.5'e'_%2; 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</syntaxhighlight> |
|||
This REXX program makes use of '''scrsize''' REXX program (or BIF) which is used to determine the screen size of the terminal (console); this is to aid in maximizing the width of the horizontal histogram. |
This REXX program makes use of '''scrsize''' REXX program (or BIF) which is used to determine the screen size of the terminal (console); this is to aid in maximizing the width of the horizontal histogram. |
||
Line 2,499: | Line 3,114: | ||
── |
── |
||
─ |
─ |
||
</pre> |
|||
=={{header|Ruby}}== |
|||
{{works with|Ruby|2.7}} |
|||
'''The Implementation''' |
|||
<syntaxhighlight lang="ruby"># Class to implement a Normal distribution, generated from a Uniform distribution. |
|||
# Uses the Marsaglia polar method. |
|||
class NormalFromUniform |
|||
# Initialize an instance. |
|||
def initialize() |
|||
@next = nil |
|||
end |
|||
# Generate and return the next Normal distribution value. |
|||
def rand() |
|||
if @next |
|||
retval, @next = @next, nil |
|||
return retval |
|||
else |
|||
u = v = s = nil |
|||
loop do |
|||
u = Random.rand(-1.0..1.0) |
|||
v = Random.rand(-1.0..1.0) |
|||
s = u**2 + v**2 |
|||
break if (s > 0.0) && (s <= 1.0) |
|||
end |
|||
f = Math.sqrt(-2.0 * Math.log(s) / s) |
|||
@next = v * f |
|||
return u * f |
|||
end |
|||
end |
|||
end</syntaxhighlight> |
|||
'''The Task''' |
|||
{{libheader|enumerable-statistics}} |
|||
<syntaxhighlight lang="ruby">require('enumerable/statistics') |
|||
def show_stats_and_histogram(data, bins) |
|||
puts("size = #{data.length} mean = #{data.mean()} stddev = #{data.stdev()}") |
|||
hist = data.histogram(bins) |
|||
scale = 100.0 / hist.weights.max |
|||
inx_beg = nil |
|||
inx_end = nil |
|||
hist.weights.length.times do |inx| |
|||
histstars = (0.5 + (scale * hist.weights[inx])).to_i |
|||
inx_beg = inx if !inx_beg && (histstars > 0) |
|||
inx_end = inx if (histstars > 0) |
|||
end |
|||
(inx_beg..inx_end).each do |inx| |
|||
bincenter = 0.5 * (hist.edges[inx] + hist.edges[inx + 1]) |
|||
histstars = (0.5 + (scale * hist.weights[inx])).to_i |
|||
puts('%6.2f: %s' % [bincenter, '*' * histstars]) |
|||
end |
|||
end |
|||
puts |
|||
puts('Uniform random number generator:') |
|||
show_stats_and_histogram(1000000.times.map { Random.rand(-1.0..1.0) }, 20) |
|||
puts |
|||
puts('Normal random numbers using the Marsaglia polar method:') |
|||
gen_normal = NormalFromUniform.new |
|||
show_stats_and_histogram(100.times.map { gen_normal.rand }, 40) |
|||
show_stats_and_histogram(10000.times.map { gen_normal.rand }, 60) |
|||
show_stats_and_histogram(1000000.times.map { gen_normal.rand }, 120)</syntaxhighlight> |
|||
{{out}} |
|||
<pre style="height: 200ex; overflow: scroll"> |
|||
Uniform random number generator: |
|||
size = 1000000 mean = 0.0001483724103528628 stddev = 0.5773085740222473 |
|||
-0.95: **************************************************************************************************** |
|||
-0.85: **************************************************************************************************** |
|||
-0.75: *************************************************************************************************** |
|||
-0.65: *************************************************************************************************** |
|||
-0.55: *************************************************************************************************** |
|||
-0.45: **************************************************************************************************** |
|||
-0.35: **************************************************************************************************** |
|||
-0.25: **************************************************************************************************** |
|||
-0.15: **************************************************************************************************** |
|||
-0.05: *************************************************************************************************** |
|||
0.05: **************************************************************************************************** |
|||
0.15: **************************************************************************************************** |
|||
0.25: **************************************************************************************************** |
|||
0.35: *************************************************************************************************** |
|||
0.45: *************************************************************************************************** |
|||
0.55: **************************************************************************************************** |
|||
0.65: **************************************************************************************************** |
|||
0.75: **************************************************************************************************** |
|||
0.85: **************************************************************************************************** |
|||
0.95: *************************************************************************************************** |
|||
Normal random numbers using the Marsaglia polar method: |
|||
size = 100 mean = 0.1611961166227389 stddev = 0.9827581078952005 |
|||
-2.30: ********** |
|||
-2.10: |
|||
-1.90: ********** |
|||
-1.70: ******************** |
|||
-1.50: |
|||
-1.30: ********** |
|||
-1.10: ********************************************************************** |
|||
-0.90: ************************************************************ |
|||
-0.70: ********************************************************************** |
|||
-0.50: ************************************************************ |
|||
-0.30: ******************************************************************************** |
|||
-0.10: ******************** |
|||
0.10: ******************************************************************************** |
|||
0.30: **************************************************************************************************** |
|||
0.50: ****************************************************************************************** |
|||
0.70: ******************************************************************************** |
|||
0.90: ************************************************************ |
|||
1.10: ****************************** |
|||
1.30: ************************************************** |
|||
1.50: |
|||
1.70: ******************** |
|||
1.90: ************************************************** |
|||
2.10: ********** |
|||
2.30: ******************** |
|||
size = 10000 mean = -0.004863294071004369 stddev = 0.9984395022628252 |
|||
-3.30: * |
|||
-3.10: * |
|||
-2.90: ** |
|||
-2.70: ** |
|||
-2.50: **** |
|||
-2.30: ********* |
|||
-2.10: *********** |
|||
-1.90: **************** |
|||
-1.70: *********************** |
|||
-1.50: ***************************** |
|||
-1.30: ***************************************** |
|||
-1.10: ************************************************* |
|||
-0.90: ****************************************************************** |
|||
-0.70: ****************************************************************************** |
|||
-0.50: *************************************************************************************** |
|||
-0.30: ********************************************************************************************* |
|||
-0.10: *********************************************************************************************** |
|||
0.10: **************************************************************************************************** |
|||
0.30: ************************************************************************************** |
|||
0.50: ************************************************************************************ |
|||
0.70: ******************************************************************************* |
|||
0.90: **************************************************************** |
|||
1.10: *************************************************** |
|||
1.30: ******************************************** |
|||
1.50: ******************************* |
|||
1.70: ********************** |
|||
1.90: **************** |
|||
2.10: ********** |
|||
2.30: ****** |
|||
2.50: ***** |
|||
2.70: ** |
|||
2.90: * |
|||
3.10: * |
|||
size = 1000000 mean = 0.0007049206911295587 stddev = 1.0000356747411552 |
|||
-3.15: * |
|||
-3.05: * |
|||
-2.95: * |
|||
-2.85: ** |
|||
-2.75: ** |
|||
-2.65: *** |
|||
-2.55: **** |
|||
-2.45: ***** |
|||
-2.35: ****** |
|||
-2.25: ******** |
|||
-2.15: ********** |
|||
-2.05: ************ |
|||
-1.95: *************** |
|||
-1.85: ****************** |
|||
-1.75: ********************* |
|||
-1.65: ************************* |
|||
-1.55: ****************************** |
|||
-1.45: *********************************** |
|||
-1.35: **************************************** |
|||
-1.25: ********************************************* |
|||
-1.15: *************************************************** |
|||
-1.05: ********************************************************* |
|||
-0.95: *************************************************************** |
|||
-0.85: ********************************************************************* |
|||
-0.75: ************************************************************************** |
|||
-0.65: ********************************************************************************* |
|||
-0.55: ************************************************************************************* |
|||
-0.45: ***************************************************************************************** |
|||
-0.35: ******************************************************************************************** |
|||
-0.25: ************************************************************************************************ |
|||
-0.15: ************************************************************************************************** |
|||
-0.05: **************************************************************************************************** |
|||
0.05: *************************************************************************************************** |
|||
0.15: ************************************************************************************************** |
|||
0.25: ************************************************************************************************ |
|||
0.35: ********************************************************************************************* |
|||
0.45: ****************************************************************************************** |
|||
0.55: ************************************************************************************* |
|||
0.65: ******************************************************************************** |
|||
0.75: ************************************************************************** |
|||
0.85: ********************************************************************* |
|||
0.95: *************************************************************** |
|||
1.05: ********************************************************* |
|||
1.15: **************************************************** |
|||
1.25: ********************************************** |
|||
1.35: **************************************** |
|||
1.45: ********************************** |
|||
1.55: ****************************** |
|||
1.65: ************************* |
|||
1.75: ********************** |
|||
1.85: ****************** |
|||
1.95: *************** |
|||
2.05: ************ |
|||
2.15: ********** |
|||
2.25: ******** |
|||
2.35: ****** |
|||
2.45: ***** |
|||
2.55: **** |
|||
2.65: *** |
|||
2.75: ** |
|||
2.85: ** |
|||
2.95: * |
|||
3.05: * |
|||
3.15: * |
|||
</pre> |
</pre> |
||
=={{header|Run BASIC}}== |
=={{header|Run BASIC}}== |
||
< |
<syntaxhighlight lang="runbasic"> |
||
s = 100000 |
s = 100000 |
||
h$ = "=============================================================" |
h$ = "=============================================================" |
||
Line 2,544: | Line 3,372: | ||
print using("##",b);" ";using("#####",bins(b));" ";left$(h$,(bins(b) / mb) * 90) |
print using("##",b);" ";using("#####",bins(b));" ";left$(h$,(bins(b) / mb) * 90) |
||
next b |
next b |
||
END</ |
END</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 2,589: | Line 3,417: | ||
= |
= |
||
</pre> |
|||
=={{header|Rust}}== |
|||
{{libheader|math}} |
|||
{{libheader|rand}} |
|||
{{libheader|rand_distr}} |
|||
<syntaxhighlight lang="rust">//! Rust rosetta example for normal distribution |
|||
use math::{histogram::Histogram, traits::ToIterator}; |
|||
use rand; |
|||
use rand_distr::{Distribution, Normal}; |
|||
/// Returns the mean of the provided samples |
|||
/// |
|||
/// ## Arguments |
|||
/// * data -- reference to float32 array |
|||
fn mean(data: &[f32]) -> Option<f32> { |
|||
let sum: f32 = data.iter().sum(); |
|||
Some(sum / data.len() as f32) |
|||
} |
|||
/// Returns standard deviation of the provided samples |
|||
/// |
|||
/// ## Arguments |
|||
/// * data -- reference to float32 array |
|||
fn standard_deviation(data: &[f32]) -> Option<f32> { |
|||
let mean = mean(data).expect("invalid mean"); |
|||
let sum = data.iter().fold(0.0, |acc, &x| acc + (x - mean).powi(2)); |
|||
Some((sum / data.len() as f32).sqrt()) |
|||
} |
|||
/// Prints a histogram in the shell |
|||
/// |
|||
/// ## Arguments |
|||
/// * data -- reference to float32 array |
|||
/// * maxwidth -- the maxwidth of the histogram in # of characters |
|||
/// * bincount -- number of bins in the histogram |
|||
/// * ch -- character used to plot the graph |
|||
fn print_histogram(data: &[f32], maxwidth: usize, bincount: usize, ch: char) { |
|||
let min_val = data.iter().cloned().fold(f32::NAN, f32::min); |
|||
let max_val = data.iter().cloned().fold(f32::NAN, f32::max); |
|||
let histogram = Histogram::new(Some(&data.to_vec()), bincount, min_val, max_val).unwrap(); |
|||
let max_bin_value = histogram.get_counters().iter().max().unwrap(); |
|||
println!(); |
|||
for x in histogram.to_iter() { |
|||
let (bin_min, bin_max, freq) = x; |
|||
let bar_width = (((freq as f64) / (*max_bin_value as f64)) * (maxwidth as f64)) as u32; |
|||
let bar_as_string = (1..bar_width).fold(String::new(), |b, _| b + &ch.to_string()); |
|||
println!( |
|||
"({:>6},{:>6}) |{} {:.2}%", |
|||
format!("{:.2}", bin_min), |
|||
format!("{:.2}", bin_max), |
|||
bar_as_string, |
|||
(freq as f64) * 100.0 / (data.len() as f64) |
|||
); |
|||
} |
|||
println!(); |
|||
} |
|||
/// Runs the demo to generate normal distribution of three different sample sizes |
|||
fn main() { |
|||
let expected_mean: f32 = 0.0; |
|||
let expected_std_deviation: f32 = 4.0; |
|||
let normal = Normal::new(expected_mean, expected_std_deviation).unwrap(); |
|||
let mut rng = rand::thread_rng(); |
|||
for &number_of_samples in &[1000, 10_000, 1_000_000] { |
|||
let data: Vec<f32> = normal |
|||
.sample_iter(&mut rng) |
|||
.take(number_of_samples) |
|||
.collect(); |
|||
println!("Statistics for sample size {}:", number_of_samples); |
|||
println!("\tMean: {:?}", mean(&data).expect("invalid mean")); |
|||
println!( |
|||
"\tStandard deviation: {:?}", |
|||
standard_deviation(&data).expect("invalid standard deviation") |
|||
); |
|||
print_histogram(&data, 80, 40, '-'); |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre style="height: 120ex; overflow: scroll"> |
|||
Statistics for sample size 1000: |
|||
Mean: -0.11356559 |
|||
Standard deviation: 4.0244555 |
|||
(-13.81,-13.12) | 0.10% |
|||
(-13.12,-12.44) | 0.00% |
|||
(-12.44,-11.75) | 0.10% |
|||
(-11.75,-11.06) | 0.20% |
|||
(-11.06,-10.38) |- 0.30% |
|||
(-10.38, -9.69) | 0.10% |
|||
( -9.69, -9.01) |--- 0.50% |
|||
( -9.01, -8.32) |---- 0.60% |
|||
( -8.32, -7.64) |------ 0.80% |
|||
( -7.64, -6.95) |-------------- 1.60% |
|||
( -6.95, -6.27) |----------------- 1.90% |
|||
( -6.27, -5.58) |------------------------ 2.60% |
|||
( -5.58, -4.90) |----------------------- 2.50% |
|||
( -4.90, -4.21) |---------------------------------------- 4.20% |
|||
( -4.21, -3.53) |------------------------------------- 3.90% |
|||
( -3.53, -2.84) |------------------------------------------------- 5.10% |
|||
( -2.84, -2.15) |---------------------------------------------- 4.80% |
|||
( -2.15, -1.47) |------------------------------------------------------------------------ 7.40% |
|||
( -1.47, -0.78) |---------------------------------------------------------- 6.00% |
|||
( -0.78, -0.10) |----------------------------------------------------------------------- 7.30% |
|||
( -0.10, 0.59) |------------------------------------------------------------------------------- 8.10% |
|||
( 0.59, 1.27) |----------------------------------------------------------------------- 7.30% |
|||
( 1.27, 1.96) |------------------------------------------------- 5.10% |
|||
( 1.96, 2.64) |------------------------------------------------------------ 6.20% |
|||
( 2.64, 3.33) |----------------------------------------- 4.30% |
|||
( 3.33, 4.01) |----------------------------- 3.10% |
|||
( 4.01, 4.70) |------------------------------------- 3.90% |
|||
( 4.70, 5.39) |-------------------------- 2.80% |
|||
( 5.39, 6.07) |---------------------- 2.40% |
|||
( 6.07, 6.76) |---------------- 1.80% |
|||
( 6.76, 7.44) |---------------- 1.80% |
|||
( 7.44, 8.13) |--------- 1.10% |
|||
( 8.13, 8.81) |---------- 1.20% |
|||
( 8.81, 9.50) | 0.20% |
|||
( 9.50, 10.18) | 0.00% |
|||
( 10.18, 10.87) | 0.10% |
|||
( 10.87, 11.55) |- 0.30% |
|||
( 11.55, 12.24) | 0.10% |
|||
( 12.24, 12.92) | 0.10% |
|||
( 12.92, 13.61) | 0.10% |
|||
Statistics for sample size 10000: |
|||
Mean: 0.02012564 |
|||
Standard deviation: 3.9697735 |
|||
(-15.80,-14.99) | 0.02% |
|||
(-14.99,-14.18) | 0.04% |
|||
(-14.18,-13.37) | 0.04% |
|||
(-13.37,-12.56) | 0.04% |
|||
(-12.56,-11.75) | 0.09% |
|||
(-11.75,-10.94) | 0.08% |
|||
(-10.94,-10.13) |- 0.25% |
|||
(-10.13, -9.32) |--- 0.42% |
|||
( -9.32, -8.51) |----- 0.67% |
|||
( -8.51, -7.70) |--------- 1.10% |
|||
( -7.70, -6.89) |------------- 1.48% |
|||
( -6.89, -6.08) |------------------ 1.98% |
|||
( -6.08, -5.27) |-------------------------- 2.82% |
|||
( -5.27, -4.45) |------------------------------------ 3.80% |
|||
( -4.45, -3.64) |--------------------------------------------- 4.66% |
|||
( -3.64, -2.83) |------------------------------------------------------- 5.72% |
|||
( -2.83, -2.02) |------------------------------------------------------------------ 6.85% |
|||
( -2.02, -1.21) |---------------------------------------------------------------------------- 7.80% |
|||
( -1.21, -0.40) |---------------------------------------------------------------------------- 7.79% |
|||
( -0.40, 0.41) |------------------------------------------------------------------------------- 8.09% |
|||
( 0.41, 1.22) |----------------------------------------------------------------------------- 7.89% |
|||
( 1.22, 2.03) |--------------------------------------------------------------------------- 7.76% |
|||
( 2.03, 2.84) |-------------------------------------------------------------------- 7.06% |
|||
( 2.84, 3.65) |------------------------------------------------------- 5.74% |
|||
( 3.65, 4.46) |-------------------------------------------- 4.64% |
|||
( 4.46, 5.27) |-------------------------------------- 4.00% |
|||
( 5.27, 6.08) |---------------------------- 3.03% |
|||
( 6.08, 6.89) |------------------- 2.07% |
|||
( 6.89, 7.71) |-------------- 1.54% |
|||
( 7.71, 8.52) |-------- 0.97% |
|||
( 8.52, 9.33) |----- 0.61% |
|||
( 9.33, 10.14) |-- 0.36% |
|||
( 10.14, 10.95) |- 0.25% |
|||
( 10.95, 11.76) | 0.19% |
|||
( 11.76, 12.57) | 0.08% |
|||
( 12.57, 13.38) | 0.02% |
|||
( 13.38, 14.19) | 0.01% |
|||
( 14.19, 15.00) | 0.03% |
|||
( 15.00, 15.81) | 0.00% |
|||
( 15.81, 16.62) | 0.01% |
|||
Statistics for sample size 1000000: |
|||
Mean: -0.004743685 |
|||
Standard deviation: 4.0006065 |
|||
(-20.00,-19.02) | 0.00% |
|||
(-19.02,-18.04) | 0.00% |
|||
(-18.04,-17.06) | 0.00% |
|||
(-17.06,-16.07) | 0.00% |
|||
(-16.07,-15.09) | 0.00% |
|||
(-15.09,-14.11) | 0.01% |
|||
(-14.11,-13.13) | 0.03% |
|||
(-13.13,-12.15) | 0.06% |
|||
(-12.15,-11.16) | 0.14% |
|||
(-11.16,-10.18) |- 0.28% |
|||
(-10.18, -9.20) |--- 0.53% |
|||
( -9.20, -8.22) |------ 0.95% |
|||
( -8.22, -7.24) |----------- 1.51% |
|||
( -7.24, -6.25) |------------------ 2.40% |
|||
( -6.25, -5.27) |--------------------------- 3.48% |
|||
( -5.27, -4.29) |-------------------------------------- 4.82% |
|||
( -4.29, -3.31) |-------------------------------------------------- 6.27% |
|||
( -3.31, -2.32) |------------------------------------------------------------- 7.62% |
|||
( -2.32, -1.34) |----------------------------------------------------------------------- 8.77% |
|||
( -1.34, -0.36) |----------------------------------------------------------------------------- 9.58% |
|||
( -0.36, 0.62) |------------------------------------------------------------------------------- 9.74% |
|||
( 0.62, 1.60) |---------------------------------------------------------------------------- 9.39% |
|||
( 1.60, 2.59) |-------------------------------------------------------------------- 8.49% |
|||
( 2.59, 3.57) |---------------------------------------------------------- 7.30% |
|||
( 3.57, 4.55) |----------------------------------------------- 5.86% |
|||
( 4.55, 5.53) |----------------------------------- 4.45% |
|||
( 5.53, 6.51) |------------------------ 3.16% |
|||
( 6.51, 7.50) |---------------- 2.09% |
|||
( 7.50, 8.48) |--------- 1.34% |
|||
( 8.48, 9.46) |----- 0.81% |
|||
( 9.46, 10.44) |-- 0.46% |
|||
( 10.44, 11.42) | 0.23% |
|||
( 11.42, 12.41) | 0.11% |
|||
( 12.41, 13.39) | 0.06% |
|||
( 13.39, 14.37) | 0.02% |
|||
( 14.37, 15.35) | 0.01% |
|||
( 15.35, 16.34) | 0.00% |
|||
( 16.34, 17.32) | 0.00% |
|||
( 17.32, 18.30) | 0.00% |
|||
( 18.30, 19.28) | 0.00% |
|||
</pre> |
</pre> |
||
=={{header|SAS}}== |
=={{header|SAS}}== |
||
< |
<syntaxhighlight lang="sas">data test; |
||
n=100000; |
n=100000; |
||
twopi=2*constant('pi'); |
twopi=2*constant('pi'); |
||
Line 2,619: | Line 3,661: | ||
y 0.000023995 1.0019996 |
y 0.000023995 1.0019996 |
||
z 0.0012857 1.0056536 |
z 0.0012857 1.0056536 |
||
*/</ |
*/</syntaxhighlight> |
||
=={{header|Sidef}}== |
=={{header|Sidef}}== |
||
{{trans| |
{{trans|Raku}} |
||
< |
<syntaxhighlight lang="ruby">define τ = Num.tau |
||
func normdist (m, σ) { |
func normdist (m, σ) { |
||
Line 2,653: | Line 3,695: | ||
var part = (8 * (x - full)) |
var part = (8 * (x - full)) |
||
say (i, "\t", '█' * full, subbar[part]) |
say (i, "\t", '█' * full, subbar[part]) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 2,697: | Line 3,739: | ||
Pairs of normal numbers are generated from pairs of uniform numbers using the [https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform Box-Muller method]. A normal density is added to the histogram for comparison. See '''[http://www.stata.com/help.cgi?histogram histogram]''' in Stata help. A [https://en.wikipedia.org/wiki/Q%E2%80%93Q_plot Q-Q plot] is also drawn. |
Pairs of normal numbers are generated from pairs of uniform numbers using the [https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform Box-Muller method]. A normal density is added to the histogram for comparison. See '''[http://www.stata.com/help.cgi?histogram histogram]''' in Stata help. A [https://en.wikipedia.org/wiki/Q%E2%80%93Q_plot Q-Q plot] is also drawn. |
||
< |
<syntaxhighlight lang="stata">clear all |
||
set obs 100000 |
set obs 100000 |
||
gen u=runiform() |
gen u=runiform() |
||
Line 2,715: | Line 3,757: | ||
hist y, normal |
hist y, normal |
||
hist z, normal |
hist z, normal |
||
qqplot x z, msize(tiny)</ |
qqplot x z, msize(tiny)</syntaxhighlight> |
||
=={{header|Tcl}}== |
=={{header|Tcl}}== |
||
< |
<syntaxhighlight lang="tcl">package require Tcl 8.5 |
||
# Uses the Box-Muller transform to compute a pair of normal random numbers |
# Uses the Box-Muller transform to compute a pair of normal random numbers |
||
proc tcl::mathfunc::nrand {mean stddev} { |
proc tcl::mathfunc::nrand {mean stddev} { |
||
Line 2,756: | Line 3,798: | ||
stats 10000 |
stats 10000 |
||
puts "" |
puts "" |
||
stats 100000 20</ |
stats 100000 20</syntaxhighlight> |
||
Sample output: |
Sample output: |
||
<pre> |
<pre> |
||
Line 2,851: | Line 3,893: | ||
=={{header|VBA}}== |
=={{header|VBA}}== |
||
< |
<syntaxhighlight lang="vb">Public Sub standard_normal() |
||
Dim s() As Variant, bins(71) As Single |
Dim s() As Variant, bins(71) As Single |
||
ReDim s(20000) |
ReDim s(20000) |
||
Line 2,868: | Line 3,910: | ||
Debug.Print |
Debug.Print |
||
Next i |
Next i |
||
End Sub</ |
End Sub</syntaxhighlight>{{out}} |
||
<pre>sample size 20000 mean-5,26306310478751E-03 standard deviation 1,00355037427319 |
<pre>sample size 20000 mean-5,26306310478751E-03 standard deviation 1,00355037427319 |
||
-3,60--3,50 |
-3,60--3,50 |
||
Line 2,941: | Line 3,983: | ||
3,30-3,40 |
3,30-3,40 |
||
3,40-3,50 </pre> |
3,40-3,50 </pre> |
||
=={{header|Wren}}== |
|||
{{libheader|Wren-fmt}} |
|||
{{libheader|Wren-math}} |
|||
<syntaxhighlight lang="wren">import "random" for Random |
|||
import "./fmt" for Fmt |
|||
import "./math" for Nums |
|||
var rgen = Random.new() |
|||
// Box-Muller method from Wikipedia |
|||
var normal = Fn.new { |mu, sigma| |
|||
var u1 = rgen.float() |
|||
var u2 = rgen.float() |
|||
var mag = sigma * (-2 * u1.log).sqrt |
|||
var z0 = mag * (2 * Num.pi * u2).cos + mu |
|||
var z1 = mag * (2 * Num.pi * u2).sin + mu |
|||
return [z0, z1] |
|||
} |
|||
var N = 100000 |
|||
var NUM_BINS = 12 |
|||
var HIST_CHAR = "■" |
|||
var HIST_CHAR_SIZE = 250 |
|||
var bins = List.filled(NUM_BINS, 0) |
|||
var binSize = 0.1 |
|||
var samples = List.filled(N, 0) |
|||
var mu = 0.5 |
|||
var sigma = 0.25 |
|||
for (i in 0...N/2) { |
|||
var rns = normal.call(mu, sigma) |
|||
for (j in 0..1) { |
|||
var rn = rns[j] |
|||
var bn |
|||
if (rn < 0) { |
|||
bn = 0 |
|||
} else if (rn >= 1) { |
|||
bn = 11 |
|||
} else { |
|||
bn = (rn/binSize).floor + 1 |
|||
} |
|||
bins[bn] = bins[bn] + 1 |
|||
samples[i*2 + j] = rn |
|||
} |
|||
} |
|||
Fmt.print("Normal distribution with mean $0.2f and S/D $0.2f for $,d samples:\n", mu, sigma, N) |
|||
System.print(" Range Number of samples within that range") |
|||
for (i in 0...NUM_BINS) { |
|||
var hist = HIST_CHAR * (bins[i] / HIST_CHAR_SIZE).round |
|||
if (i == 0) { |
|||
Fmt.print(" -∞ ..< 0.00 $s $,d", hist, bins[0]) |
|||
} else if (i < NUM_BINS - 1) { |
|||
Fmt.print("$4.2f ..< $4.2f $s $,d", binSize * (i-1), binSize * i, hist, bins[i]) |
|||
} else { |
|||
Fmt.print("1.00 ... +∞ $s $,d", hist, bins[NUM_BINS - 1]) |
|||
} |
|||
} |
|||
Fmt.print("\nActual mean for these samples : $0.5f", Nums.mean(samples)) |
|||
Fmt.print("Actual S/D for these samples : $0.5f", Nums.stdDev(samples))</syntaxhighlight> |
|||
{{out}} |
|||
Specimen run: |
|||
<pre> |
|||
Normal distribution with mean 0.50 and S/D 0.25 for 100,000 samples: |
|||
Range Number of samples within that range |
|||
-∞ ..< 0.00 ■■■■■■■■■ 2,243 |
|||
0.00 ..< 0.10 ■■■■■■■■■■■■■ 3,250 |
|||
0.10 ..< 0.20 ■■■■■■■■■■■■■■■■■■■■■■■■ 5,977 |
|||
0.20 ..< 0.30 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 9,723 |
|||
0.30 ..< 0.40 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 13,104 |
|||
0.40 ..< 0.50 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 15,601 |
|||
0.50 ..< 0.60 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 15,469 |
|||
0.60 ..< 0.70 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 13,334 |
|||
0.70 ..< 0.80 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 9,659 |
|||
0.80 ..< 0.90 ■■■■■■■■■■■■■■■■■■■■■■■■ 6,119 |
|||
0.90 ..< 1.00 ■■■■■■■■■■■■■ 3,173 |
|||
1.00 ... +∞ ■■■■■■■■■ 2,348 |
|||
Actual mean for these samples : 0.50099 |
|||
Actual S/D for these samples : 0.25051 |
|||
</pre> |
|||
=={{header|zkl}}== |
=={{header|zkl}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
< |
<syntaxhighlight lang="zkl">fcn norm2{ // Box-Muller |
||
const PI2=(0.0).pi*2;; |
const PI2=(0.0).pi*2;; |
||
rnd:=(0.0).random.fp(1); // random number in [0,1), using partial application |
rnd:=(0.0).random.fp(1); // random number in [0,1), using partial application |
||
Line 2,957: | Line 4,082: | ||
b:=(v + SIG)*BINS/SIG/2; |
b:=(v + SIG)*BINS/SIG/2; |
||
if(0<=b<BINS) h[b]+=1; |
if(0<=b<BINS) h[b]+=1; |
||
};</ |
};</syntaxhighlight> |
||
Partial application: rnd() --> (0.0).random(1). Basically, the fp method fixes the call parameters, which are then used when the partial thing is run. |
Partial application: rnd() --> (0.0).random(1). Basically, the fp method fixes the call parameters, which are then used when the partial thing is run. |
||
< |
<syntaxhighlight lang="zkl">foreach i in (N/2){ v1,v2:=norm2(); accum(v1); accum(v2); } |
||
println("Samples: %,d".fmt(N)); |
println("Samples: %,d".fmt(N)); |
||
println("Mean: ", m:=sum/N); |
println("Mean: ", m:=sum/N); |
||
println("Stddev: ", (sumSq/N - m*m).sqrt()); |
println("Stddev: ", (sumSq/N - m*m).sqrt()); |
||
foreach p in (h){ println("*"*(p/SCALE)) }</ |
foreach p in (h){ println("*"*(p/SCALE)) }</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
Latest revision as of 10:29, 10 February 2024
You are encouraged to solve this task according to the task description, using any language you may know.
The Normal (or Gaussian) distribution is a frequently used distribution in statistics. While most programming languages provide a uniformly distributed random number generator, one can derive normally distributed random numbers from a uniform generator.
- The task
- Take a uniform random number generator and create a large (you decide how large) set of numbers that follow a normal (Gaussian) distribution. Calculate the dataset's mean and standard deviation, and show a histogram of the data.
- Mention any native language support for the generation of normally distributed random numbers.
- Reference
- You may refer to code in Statistics/Basic if available.
C
/*
* RosettaCode example: Statistics/Normal distribution in C
*
* The random number generator rand() of the standard C library is obsolete
* and should not be used in more demanding applications. There are plenty
* libraries with advanced features (eg. GSL) with functions to calculate
* the mean, the standard deviation, generating random numbers etc.
* However, these features are not the core of the standard C library.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#define NMAX 10000000
double mean(double* values, int n)
{
int i;
double s = 0;
for ( i = 0; i < n; i++ )
s += values[i];
return s / n;
}
double stddev(double* values, int n)
{
int i;
double average = mean(values,n);
double s = 0;
for ( i = 0; i < n; i++ )
s += (values[i] - average) * (values[i] - average);
return sqrt(s / (n - 1));
}
/*
* Normal random numbers generator - Marsaglia algorithm.
*/
double* generate(int n)
{
int i;
int m = n + n % 2;
double* values = (double*)calloc(m,sizeof(double));
double average, deviation;
if ( values )
{
for ( i = 0; i < m; i += 2 )
{
double x,y,rsq,f;
do {
x = 2.0 * rand() / (double)RAND_MAX - 1.0;
y = 2.0 * rand() / (double)RAND_MAX - 1.0;
rsq = x * x + y * y;
}while( rsq >= 1. || rsq == 0. );
f = sqrt( -2.0 * log(rsq) / rsq );
values[i] = x * f;
values[i+1] = y * f;
}
}
return values;
}
void printHistogram(double* values, int n)
{
const int width = 50;
int max = 0;
const double low = -3.05;
const double high = 3.05;
const double delta = 0.1;
int i,j,k;
int nbins = (int)((high - low) / delta);
int* bins = (int*)calloc(nbins,sizeof(int));
if ( bins != NULL )
{
for ( i = 0; i < n; i++ )
{
int j = (int)( (values[i] - low) / delta );
if ( 0 <= j && j < nbins )
bins[j]++;
}
for ( j = 0; j < nbins; j++ )
if ( max < bins[j] )
max = bins[j];
for ( j = 0; j < nbins; j++ )
{
printf("(%5.2f, %5.2f) |", low + j * delta, low + (j + 1) * delta );
k = (int)( (double)width * (double)bins[j] / (double)max );
while(k-- > 0) putchar('*');
printf(" %-.1f%%", bins[j] * 100.0 / (double)n);
putchar('\n');
}
free(bins);
}
}
int main(void)
{
double* seq;
srand((unsigned int)time(NULL));
if ( (seq = generate(NMAX)) != NULL )
{
printf("mean = %g, stddev = %g\n\n", mean(seq,NMAX), stddev(seq,NMAX));
printHistogram(seq,NMAX);
free(seq);
printf("\n%s\n", "press enter");
getchar();
return EXIT_SUCCESS;
}
return EXIT_FAILURE;
}
- Output:
mean = 0.000477941, stddev = 0.999945 (-3.05, -2.95) | 0.1% (-2.95, -2.85) | 0.1% (-2.85, -2.75) |* 0.1% (-2.75, -2.65) |* 0.1% (-2.65, -2.55) |* 0.1% (-2.55, -2.45) |** 0.2% (-2.45, -2.35) |** 0.2% (-2.35, -2.25) |*** 0.3% (-2.25, -2.15) |**** 0.4% (-2.15, -2.05) |***** 0.4% (-2.05, -1.95) |****** 0.5% (-1.95, -1.85) |******** 0.7% (-1.85, -1.75) |********* 0.8% (-1.75, -1.65) |*********** 0.9% (-1.65, -1.55) |************* 1.1% (-1.55, -1.45) |**************** 1.3% (-1.45, -1.35) |****************** 1.5% (-1.35, -1.25) |********************* 1.7% (-1.25, -1.15) |************************ 1.9% (-1.15, -1.05) |*************************** 2.2% (-1.05, -0.95) |****************************** 2.4% (-0.95, -0.85) |********************************* 2.7% (-0.85, -0.75) |************************************ 2.9% (-0.75, -0.65) |*************************************** 3.1% (-0.65, -0.55) |***************************************** 3.3% (-0.55, -0.45) |******************************************** 3.5% (-0.45, -0.35) |********************************************** 3.7% (-0.35, -0.25) |*********************************************** 3.8% (-0.25, -0.15) |************************************************* 3.9% (-0.15, -0.05) |************************************************* 4.0% (-0.05, 0.05) |************************************************** 4.0% ( 0.05, 0.15) |************************************************* 4.0% ( 0.15, 0.25) |************************************************* 3.9% ( 0.25, 0.35) |*********************************************** 3.8% ( 0.35, 0.45) |********************************************** 3.7% ( 0.45, 0.55) |******************************************** 3.5% ( 0.55, 0.65) |***************************************** 3.3% ( 0.65, 0.75) |*************************************** 3.1% ( 0.75, 0.85) |************************************ 2.9% ( 0.85, 0.95) |********************************* 2.7% ( 0.95, 1.05) |****************************** 2.4% ( 1.05, 1.15) |*************************** 2.2% ( 1.15, 1.25) |************************ 1.9% ( 1.25, 1.35) |********************* 1.7% ( 1.35, 1.45) |****************** 1.5% ( 1.45, 1.55) |**************** 1.3% ( 1.55, 1.65) |************* 1.1% ( 1.65, 1.75) |*********** 0.9% ( 1.75, 1.85) |********* 0.8% ( 1.85, 1.95) |******** 0.7% ( 1.95, 2.05) |****** 0.5% ( 2.05, 2.15) |***** 0.4% ( 2.15, 2.25) |**** 0.4% ( 2.25, 2.35) |*** 0.3% ( 2.35, 2.45) |** 0.2% ( 2.45, 2.55) |** 0.2% ( 2.55, 2.65) |* 0.1% ( 2.65, 2.75) |* 0.1% ( 2.75, 2.85) |* 0.1% ( 2.85, 2.95) | 0.1% press enter
C#
using System;
using MathNet.Numerics.Distributions;
using MathNet.Numerics.Statistics;
class Program
{
static void RunNormal(int sampleSize)
{
double[] X = new double[sampleSize];
var norm = new Normal(new Random());
norm.Samples(X);
const int numBuckets = 10;
var histogram = new Histogram(X, numBuckets);
Console.WriteLine("Sample size: {0:N0}", sampleSize);
for (int i = 0; i < numBuckets; i++)
{
string bar = new String('#', (int)(histogram[i].Count * 360 / sampleSize));
Console.WriteLine(" {0:0.00} : {1}", histogram[i].LowerBound, bar);
}
var statistics = new DescriptiveStatistics(X);
Console.WriteLine(" Mean: " + statistics.Mean);
Console.WriteLine("StdDev: " + statistics.StandardDeviation);
Console.WriteLine();
}
static void Main(string[] args)
{
RunNormal(100);
RunNormal(1000);
RunNormal(10000);
}
}
- Output:
Sample size: 100 -2.12 : ####### -1.66 : ############################ -1.19 : ####################################### -0.72 : ############################################## -0.26 : ############################################################################### 0.21 : ###################################################################################### 0.68 : ################################ 1.14 : ######################### 1.61 : ### 2.07 : ########## Mean: 0.0394411345297757 StdDev: 0.925286665513647 Sample size: 1,000 -2.98 : ## -2.34 : ########## -1.69 : ############################## -1.05 : ################################################################ -0.40 : ########################################################################################### 0.24 : ######################################################################################## 0.88 : ############################################## 1.53 : ################## 2.17 : ##### 2.82 : ## Mean: 0.0868718238400114 StdDev: 0.989120264661867 Sample size: 10,000 -3.88 : -3.12 : ## -2.35 : ################# -1.59 : #################################################### -0.82 : ################################################################################################ -0.06 : #################################################################################################### 0.71 : ############################################################### 1.47 : ##################### 2.23 : #### 3.00 : Mean: 0.0208920122989818 StdDev: 1.00046328880424
C++
showing features of C++11 here
#include <random>
#include <map>
#include <string>
#include <iostream>
#include <cmath>
#include <iomanip>
int main( ) {
std::random_device myseed ;
std::mt19937 engine ( myseed( ) ) ;
std::normal_distribution<> normDistri ( 2 , 3 ) ;
std::map<int , int> normalFreq ;
int sum = 0 ; //holds the sum of the randomly created numbers
double mean = 0.0 ;
double stddev = 0.0 ;
for ( int i = 1 ; i < 10001 ; i++ )
++normalFreq[ normDistri ( engine ) ] ;
for ( auto MapIt : normalFreq ) {
sum += MapIt.first * MapIt.second ;
}
mean = sum / 10000 ;
stddev = sqrt( sum / 10000 ) ;
std::cout << "The mean of the distribution is " << mean << " , the " ;
std::cout << "standard deviation " << stddev << " !\n" ;
std::cout << "And now the histogram:\n" ;
for ( auto MapIt : normalFreq ) {
std::cout << std::left << std::setw( 4 ) << MapIt.first <<
std::string( MapIt.second / 100 , '*' ) << std::endl ;
}
return 0 ;
}
Output:
The mean of the distribution is 1 , the standard deviation 1 ! And now the histogram: -10 -9 -8 -7 -6 -5 -4 * -3 ** -2 **** -1 ****** 0 ********************* 1 ************ 2 ************ 3 *********** 4 ********* 5 ****** 6 **** 7 ** 8 * 9 10 11 12 13
D
This uses the Box-Muller method as in the Go entry, and the module from the Statistics/Basic. A ziggurat-based normal generator for the Phobos standard library is in the works.
import std.stdio, std.random, std.math, std.range, std.algorithm,
statistics_basic;
struct Normals {
double mu, sigma;
double[2] state;
size_t index = state.length;
enum empty = false;
void popFront() pure nothrow { index++; }
@property double front() {
if (index >= state.length) {
immutable r = sqrt(-2 * uniform!"]["(0., 1.0).log) * sigma;
immutable x = 2 * PI * uniform01;
state = [mu + r * x.sin, mu + r * x.cos];
index = 0;
}
return state[index];
}
}
void main() {
const data = Normals(0.0, 0.5).take(100_000).array;
writefln("Mean: %8.6f, SD: %8.6f\n", data.meanStdDev[]);
data.map!q{ max(0.0, min(0.9999, a / 3 + 0.5)) }.showHistogram01;
}
- Output:
Mean: 0.000528, SD: 0.502245 0.0: * 0.1: ****** 0.2: ***************** 0.3: *********************************** 0.4: ************************************************* 0.5: ************************************************** 0.6: ********************************** 0.7: ***************** 0.8: ****** 0.9: *
EDSAC order code
[Normal distribution for Rosetta Code
EDSAC program, Initial Orders 2]
[==================================================================================
Uses an accept-reject method, which requires only logarithms (no trig functions).
Let u, v be independent uniform variates in (0, 1). Let x = -ln(u), y = -ln(v).
Accept x iff (x - 1)^2 <= 2y. If x is accepted, negate it with probability 1/2.
Then x is normally distributed with mean 0 and standard deviation 1.
The algorithm is modified for this EDSAC version:
(1) Uses EDSAC library subroutine L1 to calculate (1/32)log_2() rather than ln()
(2) Since real numbers on EDSAC are restricted to the interval [-1, 1), scales so
that the standard deviation is 1/4, and reject values >= 4 s.d. from the mean.
In the histogram, counts are divided by 16, with rounding.
On the EDSAC PC simulator, takes 2.75 EDSAC hours to find 4096 normal variates.
[=================================================================================]
[Arrange the storage]
T46K P70F [N parameter: library subroutine P1 to print +ve fraction]
T47K P200F [M parameter: main routine and dependent subroutines]
T49K P92F [L parameter: library subroutine L1 to calculate log_2]
T51K P130F [G parameter: generator for pseudo-random numbers]
T52K P168F [A parameter: library subroutine S2 for square root]
T54K P190F [C parameter: constants read in by library subroutine R9]
[Library subroutine R9, reads non-negative integers at load time.
Fractions can be read by converting to integers (multiply by 2^34).
15 locations, must be loaded at location 56.]
T56KGKT20FVDL8FA40DUDTFI40FA40FS39FG@S2FG23FA5@T5@E4@
[Library subroutine M3, prints header at load time.
M3 and header are then overwritten.]
PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF
*MEAN#!K0BL!!*SD#!K0M25BL@&#..PZ
[========================== C parameter ==========================]
[Tell R9 where to store integers read from tape]
E69K T#C ['T m D' in WWG, but this also works]
[List of integers; separated by F; end list with #TZ]
123456789F5F1549082005F11908177887F858993#TZ
[0] [seed for random number generator]
[2] [minimum argument for library subroutine L1]
[4] [1/(16*ln(2)) for accept-reject]
[6] [ln(2) for scaling standard deviation]
[8] [0.00005 for rounding in print routine]
[========================== M parameter ==========================]
E25K TM GK
[0] P4096F [number of data, in address field; code below assumes 4096]
[1] PF [negative count of data]
[2] PF [worlspace, low word]
[3] PF [workspace, high word]
[4] PF PF [auxiliary for accept-reject algorithm]
[6] PF PF [sum of value/count, for mean]
[8] PF PF [sum of value^2/count, for variance]
[10] PFPFPFPFPFPFPFPFPFPFPFPFPFPFPFPF [16 bins for histogram]
[26] CF [11110...0 binary, to isolate bin index; also prints colon]
[27] A10@ [A order for bin{0}]
[28] AF [A order for bin{16}]
[29] P8F [8 in address field]
[30] MF [subtract from A order to make T order; also prints dot]
[31] WF [1/8]
[32] FF [-15/16, mid value of histogram bin{0}]
[33] PF [mid value of current bin]
[Teleprinter]
[34] #F [set figures mode]
[35] PF
[36] AF [minus (in figures mode)]
[37] !F [space]
[38] @F [carriage return]
[39] &F [line feed]
[Enter with acc = 0]
[40]
S@ T1@ [store negated data count in address field]
A#C T4D [copy seed to 4D for PRNG]
[44] A44@ GG [initialize PRNG]
T6#@ T8#@ [clear sum and sum of squares]
O34@ [set teleprinter to figures]
[Start of loop to generate normal variates]
[49] TF [clear acc]
[Calculate and store logarithms of uniform variates]
[50] A50@ G176@ SD T2#@ [corresponds to x = -ln(u)]
[54] A54@ G176@ SD T4#@ [corresponds to y = -ln(v)]
[Accept or reject]
A4#C RD [acc := (1/32*ln(2))]
S2#@ [subtract x]
TD HD ND [acc := negated square]
H4#C V4#@ [add 2*(auxiliary value)]
G49@ [reject if result < 0]
[First log is accepted; multiply by 8*ln(2) so that s.d. becomes 0.25]
[Reject values outside (-1,1) (that is >= 4 s.d. from mean, unlikely)]
T6F [clear acc]
H6#C V2#@ [times ln(2)]
S31@ E49@ [if product >= 1/8 (unlikely) reject and try again]
A31@ [restore acc after test]
L2F [shift 3 left to complete scaling]
YF T2#@ [round and store back]
[Finally change sign with probability 1/2]
T6F T4D A2F T4F [pass range = 2 to PRNG]
[80] A80@ G1G [call PRNG, 0 or 1 returned in 0D]
SD E87@ [don't change sign if 0D = 0]
T6F S2#@ T2#@ [change sign, so value < 0]
[87] [Here with acc = 0.]
[To print the values, replace the following jump with a no-op (X F)]
E94@
A2#@ TD [pass abs(value) to print routine]
[90] A90@ G191@
O38@ O39@ [print CR, LF]
[94] A2#@ R1024F YF [load value, divide by count, round]
A6#@ T6#@ [update sum]
H2#@ V2#@ R1024F YF
A8#@ T8#@ [update sum of squares]
[Inc count in appropriate bin]
H26@ [mult reg := mask to isolate bin index]
C3@ [acc := top 4 bits of value]
R1024F [12 right, get bin -8..7 in address field]
A29@ [add 8 to address, bin index now 0..15]
A27@ [make A order for bin]
U113@ [plant in code]
S30@ [convert to T order]
T115@ [plant in code]
[113] AF [(planted) acc := count in bin]
A2F [inc by 1 in address field]
[115] TF [(planted) store updated bin count]
A1@ A2F U1@ [inc negative count of variates]
G49@ [loop till got required number]
[Print mean and standard deviation]
A6#@ TD [pass mean to print subroutine]
[122] A122@ G191@
O37@ O37@ O37@ [print spaces]
A8#@ H6#@
N6#@ T4D [calc variance, pass to square root s/r]
[131] A131@ GA [4D := standard deviation]
A4D U8#@ TD [pass to print subroutine]
[136] A136@ G191@ [print s.d.]
O38@ O39@ O38@ O39@ [print CR, LF twice]
[Print histogram]
A29@ LD A27@ [make A order for exclusive end bin]
T28@ [store as test for end]
A32@ T33@ [initialize mid-value of bin]
A27@ [A order for bin{0}]
[149] T158@ [loop: plant A order for current bin]
TD A33@ U1F [0D = mid-value of bin, extended to 35 bits]
A31@ T33@ [update mid-value for next time]
[155] A155@ G191@ [print mid-value]
O26@ [print colon]
[158] AF [(planted) load number of hits in address field]
A29@ [add 8 for rounding]
R4F [divide by 16]
E163@ [jump to middle of loop]
[162] O175@ [print plus sign]
[163] S2F E162@ [loop till printed enough plus signs]
O38@ O39@ [print CR, LF]
TF [clear acc]
A158@ A2F [make A order for next bin]
S28@ [compare with end order]
E174@ [exit if no more bins]
A28@ [restore acc after test]
G149@ [loop back (A order is negative)]
[174] O34@ [dummy character to flush print buffer]
[175] ZF [halt program; also serves as plus sign]
[Subroutine of main routine. Sets 0D := (1/32)log_2 of uniform variate]
[176] A3F T188@ [plant return link as usual]
[178] T4D [pass range = 0 to PRNG]
[179] A179@ G1G [0D := uniform variate]
AD S2#C [test for too small (unlikely)]
G189@ [jump to try again if so]
A2#C T6D [OK, pass uniform variate to logarithm subroutine]
[186] A186@ GL [0D := logarithm]
[188] ZF [(planted) return to caller]
[189] TF E178@ [if failed, clear acc and try again]
[Subroutine of main routine; prints signed fraction in 0D to 4 decimals.]
[Wrapper for library subroutine P1; adds sign, rounding, layout.]
[191] A3F T207@ [plant return link as usual]
AD G197@ [acc := value, jump if < 0]
O37@ E200@ [print space, jump to common code]
[197] O36@ [value < 0, print minus]
TD SD [acc := abs(valus)]
[200] A8#C TD [add 0.00005 for rounding, pass to P1 in 0D]
O35@ O30@ [print '0.']
[204] A204@ GN P4F [call P1 to print 5 decimals]
[207] ZF [(planted) jump back to caller]
[========================== G parameter ==========================]
[Linear congruential generator, same algorithm as Delphi 7 LCG.]
[38 storage locations.]
[Initialize: Call 0G with 35-bit seed in 4D.]
[Input: Call 1G with 35-bit range in 4D.]
[Output: if range > 0, 0D := random 35-bit integer in 0..(range - 1).]
[if range = 0, 0D := random 35-bit real in [0, 1)]
E25K TG
GKG10@G15@T2#ZPFT2ZI514DP257FT4#ZPFT4ZPDPFT6#ZPFT6ZPFRFA6#@S4#@
T6#@E25FE8ZPFT8ZPFPFA3FT14@A4DT8#@ZFA3FT37@H2#@V8#@L512FL512FL1024F
A4#@T8#@H6#@C8#@T8#@S4DG32@TDA8#@E35@H4DTDV8#@L1FTDZF
[==================== Library subroutines ============================]
E25K TL
[L1: Logarithm to base 2.]
[0D := (1/32)*log_2(6D), provided 6D > 2^(-32)]
[38 storage locations; working positions 4D and 8F.]
GKA3FT33@E11@IFP1024FP512FA3@LDT6DADS4@TDS3@A6DG6@T8FS5@T4D
H6DV6DS3@E34@A3@LDYFT6DA4DADTDA4DRDYFG17@EFA3@YFT6DE29@
E25K TN
[P1: Print non-negative fraction in 0D, without layout or rounding]
[21 storage locations.]
GKA18@U17@S20@T5@H19@PFT5@VDUFOFFFSFL4FTDA5@A2FG6@EFU3FJFM1F
E25K TA
[S2: square root.]
[Closed: 22 storage locations, working position 0D.]
[Forms sqrt( C(4D)) where C(4D) > 0, and places result in 4D.]
GKA3FT20@A4DS9@A6@UDHDR1FS21@TDN4DRDA4DYFT4DVDTDVDYFG5@EFSF
[======================= M parameter again ======================]
E25K TM GK
E40Z [define entry point]
PF [acc = 0 on entry]
- Output:
MEAN (0?) SD (0.25?) -0.0015 0.2488 -0.9375: -0.8125:+ -0.6875:++ -0.5625:++++ -0.4375:+++++++++++ -0.3125:++++++++++++++++++++++ -0.1875:++++++++++++++++++++++++++++++++++++++++ -0.0625:++++++++++++++++++++++++++++++++++++++++++++++++ 0.0625:++++++++++++++++++++++++++++++++++++++++++++++++++ 0.1875:++++++++++++++++++++++++++++++++++++++++ 0.3125:+++++++++++++++++++++++ 0.4375:+++++++++++ 0.5625:++++ 0.6875:+ 0.8125: 0.9375:
Elixir
defmodule Statistics do
def normal_distribution(n, w\\5) do
{sum, sum2, hist} = generate(n, w)
mean = sum / n
stddev = :math.sqrt(sum2 / n - mean*mean)
IO.puts "size: #{n}"
IO.puts "mean: #{mean}"
IO.puts "stddev: #{stddev}"
{min, max} = Map.to_list(hist)
|> Enum.filter_map(fn {_k,v} -> v >= n/120/w end, fn {k,_v} -> k end)
|> Enum.min_max
Enum.each(min..max, fn i ->
bar = String.duplicate("=", trunc(120 * w * Map.get(hist, i, 0) / n))
:io.fwrite "~4.1f: ~s~n", [i/w, bar]
end)
IO.puts ""
end
defp generate(n, w) do
Enum.reduce(1..n, {0, 0, %{}}, fn _,{sum, sum2, hist} ->
z = :rand.normal
{sum+z, sum2+z*z, Map.update(hist, round(w*z), 1, &(&1+1))}
end)
end
end
Enum.each([100,1000,10000], fn n ->
Statistics.normal_distribution(n)
end)
- Output:
size: 100 mean: 0.027742416007234007 stddev: 1.0209597927405403 -2.6: ============ -2.4: -2.2: ============ -2.0: ====== -1.8: -1.6: -1.4: ============================== -1.2: ====== -1.0: ============================== -0.8: ========================================== -0.6: ========================================== -0.4: ================================================ -0.2: ================================================ 0.0: ============================== 0.2: ==================================== 0.4: ========================================== 0.6: ====================================================== 0.8: ========================================== 1.0: ================================================ 1.2: ============================== 1.4: ====== 1.6: ============ 1.8: ============ 2.0: 2.2: 2.4: ====== 2.6: ====== size: 1000 mean: -0.025562168667763084 stddev: 1.0338288521306742 -3.2: = -3.0: -2.8: = -2.6: === -2.4: == -2.2: ====== -2.0: == -1.8: ============= -1.6: =============== -1.4: ================= -1.2: ================= -1.0: ==================================== -0.8: =================================== -0.6: ============================================ -0.4: ============================================ -0.2: =============================================== 0.0: ========================================= 0.2: =========================================== 0.4: ============================================= 0.6: ======================================= 0.8: ================================ 1.0: ============================ 1.2: ======================== 1.4: ================== 1.6: ========== 1.8: ===== 2.0: ======== 2.2: ==== 2.4: ===== 2.6: = 2.8: = size: 10000 mean: -0.009132420943007152 stddev: 0.9979508347451509 -2.6: = -2.4: === -2.2: ==== -2.0: ===== -1.8: ========= -1.6: ============== -1.4: ================ -1.2: ======================= -1.0: ============================ -0.8: ================================= -0.6: ============================================ -0.4: =========================================== -0.2: ============================================== 0.0: ================================================== 0.2: ============================================ 0.4: =========================================== 0.6: ======================================= 0.8: ===================================== 1.0: ============================ 1.2: ======================= 1.4: ================ 1.6: ============== 1.8: ========= 2.0: ====== 2.2: === 2.4: == 2.6: =
Factor
USING: assocs formatting kernel math math.functions
math.statistics random sequences sorting ;
2,000,000 [ 0 1 normal-random-float ] replicate ! make data set
dup [ mean ] [ population-std ] bi ! calculate and show
"Mean: %f\nStdev: %f\n\n" printf ! mean and stddev
[ 10 * floor 10 / ] map ! map data to buckets
histogram >alist [ first ] sort-with ! create histogram sorted by bucket (key)
dup values supremum ! find maximum count
[
[ /f 100 * >integer ] keepd ! how big should this histogram bar be?
[ [ CHAR: * ] "" replicate-as ] dip ! make the bar
"% 5.2f: %s %d\n" printf ! print a line of the histogram
] curry assoc-each
- Output:
Mean: 0.000798 Stdev: 1.000549 -4.90: 2 -4.80: 1 -4.70: 1 -4.60: 3 -4.50: 3 -4.40: 6 -4.30: 15 -4.20: 13 -4.10: 16 -4.00: 42 -3.90: 62 -3.80: 68 -3.70: 98 -3.60: 145 -3.50: 205 -3.40: 311 -3.30: 379 -3.20: 580 -3.10: 739 -3.00: * 1002 -2.90: * 1349 -2.80: ** 1893 -2.70: *** 2499 -2.60: **** 3211 -2.50: ***** 4035 -2.40: ****** 5141 -2.30: ******* 6392 -2.20: ********* 7869 -2.10: ************ 9780 -2.00: ************** 11787 -1.90: ****************** 14483 -1.80: ********************* 17183 -1.70: ************************* 20387 -1.60: ****************************** 24049 -1.50: ********************************** 27555 -1.40: **************************************** 32153 -1.30: ********************************************* 36707 -1.20: *************************************************** 40921 -1.10: ********************************************************* 45928 -1.00: *************************************************************** 50707 -0.90: ********************************************************************* 55697 -0.80: *************************************************************************** 60377 -0.70: ******************************************************************************** 64358 -0.60: ************************************************************************************ 67928 -0.50: ***************************************************************************************** 71911 -0.40: ********************************************************************************************* 75054 -0.30: ************************************************************************************************ 77073 -0.20: ************************************************************************************************** 78768 -0.10: *************************************************************************************************** 79732 0.00: **************************************************************************************************** 79952 0.10: *************************************************************************************************** 79412 0.20: ************************************************************************************************ 77511 0.30: ********************************************************************************************* 74487 0.40: ****************************************************************************************** 72250 0.50: ************************************************************************************** 68789 0.60: ******************************************************************************** 64408 0.70: *************************************************************************** 60122 0.80: ********************************************************************* 55619 0.90: *************************************************************** 50869 1.00: ********************************************************* 45883 1.10: **************************************************** 41586 1.20: ********************************************** 37145 1.30: *************************************** 31715 1.40: ********************************** 27779 1.50: ****************************** 24270 1.60: ************************* 20516 1.70: ********************* 17221 1.80: ***************** 14344 1.90: ************** 11789 2.00: ************ 9796 2.10: ********* 7922 2.20: ******* 6331 2.30: ****** 5138 2.40: ***** 4044 2.50: *** 3065 2.60: ** 2397 2.70: ** 1846 2.80: * 1462 2.90: * 1001 3.00: 765 3.10: 587 3.20: 393 3.30: 299 3.40: 197 3.50: 132 3.60: 100 3.70: 74 3.80: 59 3.90: 32 4.00: 29 4.10: 12 4.20: 15 4.30: 6 4.40: 3 4.50: 4 4.60: 3 4.70: 2 4.80: 1
Fortran
Using the Marsaglia polar method
program Normal_Distribution
implicit none
integer, parameter :: i64 = selected_int_kind(18)
integer, parameter :: r64 = selected_real_kind(15)
integer(i64), parameter :: samples = 1000000_i64
real(r64) :: mean, stddev
real(r64) :: sumn = 0, sumnsq = 0
integer(i64) :: n = 0
integer(i64) :: bin(-50:50) = 0
integer :: i, ind
real(r64) :: ur1, ur2, nr1, nr2, s
n = 0
do while(n <= samples)
call random_number(ur1)
call random_number(ur2)
ur1 = ur1 * 2.0 - 1.0
ur2 = ur2 * 2.0 - 1.0
s = ur1*ur1 + ur2*ur2
if(s >= 1.0_r64) cycle
nr1 = ur1 * sqrt(-2.0*log(s)/s)
ind = floor(5.0*nr1)
bin(ind) = bin(ind) + 1_i64
sumn = sumn + nr1
sumnsq = sumnsq + nr1*nr1
nr2 = ur2 * sqrt(-2.0*log(s)/s)
ind = floor(5.0*nr2)
bin(ind) = bin(ind) + 1_i64
sumn = sumn + nr2
sumnsq = sumnsq + nr2*nr2
n = n + 2_i64
end do
mean = sumn / n
stddev = sqrt(sumnsq/n - mean*mean)
write(*, "(a, i0)") "sample size = ", samples
write(*, "(a, f17.15)") "Mean : ", mean,
write(*, "(a, f17.15)") "Stddev : ", stddev
do i = -15, 15
write(*, "(f4.1, a, a)") real(i)/5.0, ": ", repeat("=", int(bin(i)*500/samples))
end do
end program
- Output:
sample size = 1000 Mean : 0.043096320705032 Stddev : 0.981532585231540 -3.0: -2.8: -2.6: == -2.4: == -2.2: ==== -2.0: ====== -1.8: ======= -1.6: ============ -1.4: ================ -1.2: ===================== -1.0: =========================== -0.8: ======================= -0.6: ================================== -0.4: ===================================== -0.2: ========================================== 0.0: =============================================== 0.2: ==================================== 0.4: ================================= 0.6: ================================== 0.8: ============================= 1.0: ==================== 1.2: ========================== 1.4: =========== 1.6: ========= 1.8: ==== 2.0: ====== 2.2: === 2.4: 2.6: 2.8: = 3.0: sample size = 1000000 Mean : 0.000166653231289 Stddev : 1.000025612171690 -3.0: -2.8: = -2.6: = -2.4: == -2.2: ==== -2.0: ====== -1.8: ========= -1.6: ============ -1.4: ================= -1.2: ===================== -1.0: ========================== -0.8: =============================== -0.6: =================================== -0.4: ====================================== -0.2: ======================================= 0.0: ======================================= 0.2: ====================================== 0.4: ================================== 0.6: =============================== 0.8: ========================== 1.0: ===================== 1.2: ================= 1.4: ============ 1.6: ========= 1.8: ====== 2.0: ==== 2.2: == 2.4: = 2.6: = 2.8: 3.0:
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
Sub normalStats(sampleSize As Integer)
If sampleSize < 1 Then Return
Dim r(1 To sampleSize) As Double
Dim h(-1 To 10) As Integer '' all zero by default
Dim sum As Double = 0.0
Dim hSum As Integer = 0
' Generate 'sampleSize' normally distributed random numbers with mean 0.5 and standard deviation 0.25
' calculate their sum
' and in which box they will fall when drawing the histogram
For i As Integer = 1 To sampleSize
r(i) = 0.5 + randomNormal / 4.0
sum += r(i)
If r(i) < 0.0 Then
h(-1) += 1
ElseIf r(i) >= 1.0 Then
h(10) += 1
Else
h(Int(r(i) * 10)) += 1
End If
Next
For i As Integer = -1 To 10 : hSum += h(i) : Next
' adjust one of the h() values if necessary to ensure hSum = sampleSize
Dim adj As Integer = sampleSize - hSum
If adj <> 0 Then
For i As Integer = -1 To 10
h(i) += adj
If h(i) >= 0 Then Exit For
h(i) -= adj
Next
End If
Dim mean As Double = sum / sampleSize
Dim sd As Double
sum = 0.0
' Now calculate their standard deviation
For i As Integer = 1 To sampleSize
sum += (r(i) - mean) ^ 2.0
Next
sd = Sqr(sum/sampleSize)
' Draw a histogram of the data with interval 0.1
Dim numStars As Integer
' If sample size > 300 then normalize histogram to 300
Dim scale As Double = 1.0
If sampleSize > 300 Then scale = 300.0 / sampleSize
Print "Sample size "; sampleSize
Print
Print Using " Mean #.######"; mean;
Print Using " SD #.######"; sd
Print
For i As Integer = -1 To 10
If i = -1 Then
Print Using "< 0.00 : ";
ElseIf i = 10 Then
Print Using ">=1.00 : ";
Else
Print Using " #.## : "; i/10.0;
End If
Print Using "##### " ; h(i);
numStars = Int(h(i) * scale + 0.5)
Print String(numStars, "*")
Next
End Sub
normalStats 100
Print
normalStats 1000
Print
normalStats 10000
Print
normalStats 100000
Print
Print "Press any key to quit"
Sleep
Sample output:
- Output:
Sample size 100 Mean 0.486977 SD 0.244147 < 0.00 : 2 ** 0.00 : 5 ***** 0.10 : 4 **** 0.20 : 14 ************** 0.30 : 12 ************ 0.40 : 15 *************** 0.50 : 17 ***************** 0.60 : 11 *********** 0.70 : 9 ********* 0.80 : 7 ******* 0.90 : 1 * >=1.00 : 3 *** Sample size 1000 Mean 0.489234 SD 0.247606 < 0.00 : 18 ***** 0.00 : 32 ********** 0.10 : 73 ********************** 0.20 : 111 ********************************* 0.30 : 138 ***************************************** 0.40 : 151 ********************************************* 0.50 : 153 ********************************************** 0.60 : 114 ********************************** 0.70 : 101 ****************************** 0.80 : 51 *************** 0.90 : 38 *********** >=1.00 : 20 ****** Sample size 10000 Mean 0.498239 SD 0.249235 < 0.00 : 225 ******* 0.00 : 333 ********** 0.10 : 589 ****************** 0.20 : 999 ****************************** 0.30 : 1320 **************************************** 0.40 : 1542 ********************************************** 0.50 : 1581 *********************************************** 0.60 : 1323 **************************************** 0.70 : 963 ***************************** 0.80 : 591 ****************** 0.90 : 314 ********* >=1.00 : 220 ******* Sample size 100000 Mean 0.500925 SD 0.248910 < 0.00 : 2173 ******* 0.00 : 3192 ********** 0.10 : 5938 ****************** 0.20 : 9715 ***************************** 0.30 : 13351 **************************************** 0.40 : 15399 ********************************************** 0.50 : 15680 *********************************************** 0.60 : 13422 **************************************** 0.70 : 9633 ***************************** 0.80 : 5993 ****************** 0.90 : 3207 ********** >=1.00 : 2297 *******
Go
Box-Muller method shown here. Go has a normally distributed random function in the standard library, as shown in the Go Random numbers solution. It uses the ziggurat method.
package main
import (
"fmt"
"math"
"math/rand"
"strings"
)
// Box-Muller
func norm2() (s, c float64) {
r := math.Sqrt(-2 * math.Log(rand.Float64()))
s, c = math.Sincos(2 * math.Pi * rand.Float64())
return s * r, c * r
}
func main() {
const (
n = 10000
bins = 12
sig = 3
scale = 100
)
var sum, sumSq float64
h := make([]int, bins)
for i, accum := 0, func(v float64) {
sum += v
sumSq += v * v
b := int((v + sig) * bins / sig / 2)
if b >= 0 && b < bins {
h[b]++
}
}; i < n/2; i++ {
v1, v2 := norm2()
accum(v1)
accum(v2)
}
m := sum / n
fmt.Println("mean:", m)
fmt.Println("stddev:", math.Sqrt(sumSq/float64(n)-m*m))
for _, p := range h {
fmt.Println(strings.Repeat("*", p/scale))
}
}
Output:
mean: -0.0034970888831523488 stddev: 1.0040682925006286 * **** ********* *************** ******************* ****************** ************** ********* **** *
Haskell
import Data.Map (Map, empty, insert, findWithDefault, toList)
import Data.Maybe (fromMaybe)
import Text.Printf (printf)
import Data.Function (on)
import Data.List (sort, maximumBy, minimumBy)
import Control.Monad.Random (RandomGen, Rand, evalRandIO, getRandomR)
import Control.Monad (replicateM)
-- Box-Muller
getNorm :: RandomGen g => Rand g Double
getNorm = do
u0 <- getRandomR (0.0, 1.0)
u1 <- getRandomR (0.0, 1.0)
let r = sqrt $ (-2.0) * log u0
theta = 2.0 * pi * u1
return $ r * sin theta
putInBin :: Double -> Map Int Int -> Double -> Map Int Int
putInBin binWidth t v =
let bin = round (v / binWidth)
count = findWithDefault 0 bin t
in insert bin (count+1) t
runTest :: Int -> IO ()
runTest n = do
rs <- evalRandIO $ replicateM n getNorm
let binWidth = 0.1
tally v (sv, sv2, t) = (sv+v, sv2 + v*v, putInBin binWidth t v)
(sum, sum2, tallies) = foldr tally (0.0, 0.0, empty) rs
tallyList = sort $ toList tallies
printStars tallies binWidth maxCount selection =
let count = findWithDefault 0 selection tallies
bin = binWidth * fromIntegral selection
maxStars = 100
starCount = if maxCount <= maxStars
then count
else maxStars * count `div` maxCount
stars = replicate starCount '*'
in printf "%5.2f: %s %d\n" bin stars count
mean = sum / fromIntegral n
stddev = sqrt (sum2/fromIntegral n - mean*mean)
printf "\n"
printf "sample count: %d\n" n
printf "mean: %9.7f\n" mean
printf "stddev: %9.7f\n" stddev
let maxCount = snd $ maximumBy (compare `on` snd) tallyList
maxBin = fst $ maximumBy (compare `on` fst) tallyList
minBin = fst $ minimumBy (compare `on` fst) tallyList
mapM_ (printStars tallies binWidth maxCount) [minBin..maxBin]
main = do
runTest 1000
runTest 2000000
- Output:
sample count: 1000 mean: -0.0269949 stddev: 0.9795285 -3.10: ** 2 -3.00: 0 -2.90: 0 -2.80: ** 2 -2.70: * 1 -2.60: **** 4 -2.50: ** 2 -2.40: ** 2 -2.30: 0 -2.20: *** 3 -2.10: ***** 5 -2.00: ****** 6 -1.90: ****** 6 -1.80: *********** 11 -1.70: ************ 12 -1.60: ******* 7 -1.50: ************* 13 -1.40: ***************** 17 -1.30: ******************** 20 -1.20: **************** 16 -1.10: ***************** 17 -1.00: ********************** 22 -0.90: *************************** 27 -0.80: ********************** 22 -0.70: ******************************** 32 -0.60: ********************************* 33 -0.50: ****************************************** 42 -0.40: *********************************************** 47 -0.30: ************************************************ 48 -0.20: *************************** 27 -0.10: ***************************** 29 0.00: *********************************************** 47 0.10: *************************************************** 51 0.20: ****************************************** 42 0.30: ******************************** 32 0.40: ********************************* 33 0.50: ***************************************** 41 0.60: ****************************************** 42 0.70: **************************** 28 0.80: ********************** 22 0.90: *************************** 27 1.00: ******************* 19 1.10: ********************** 22 1.20: ************************ 24 1.30: ******************** 20 1.40: ***************** 17 1.50: ********** 10 1.60: ************* 13 1.70: **** 4 1.80: *** 3 1.90: ******* 7 2.00: ****** 6 2.10: * 1 2.20: * 1 2.30: ******* 7 2.40: *** 3 2.50: 0 2.60: * 1 2.70: 0 2.80: 0 2.90: 0 3.00: * 1 3.10: 0 3.20: 0 3.30: * 1 sample count: 2000000 mean: 0.0001017 stddev: 0.9994329 -4.60: 3 -4.50: 2 -4.40: 3 -4.30: 9 -4.20: 18 -4.10: 19 -4.00: 20 -3.90: 41 -3.80: 77 -3.70: 84 -3.60: 105 -3.50: 189 -3.40: 245 -3.30: 350 -3.20: 460 -3.10: 619 -3.00: * 838 -2.90: * 1234 -2.80: * 1586 -2.70: ** 2063 -2.60: *** 2716 -2.50: **** 3503 -2.40: ***** 4345 -2.30: ******* 5678 -2.20: ******** 7160 -2.10: *********** 8856 -2.00: ************* 10915 -1.90: **************** 13299 -1.80: ******************* 15599 -1.70: *********************** 19004 -1.60: *************************** 22321 -1.50: ******************************** 25940 -1.40: ************************************* 29622 -1.30: ****************************************** 34213 -1.20: ************************************************ 38922 -1.10: ****************************************************** 43415 -1.00: ************************************************************ 48250 -0.90: ****************************************************************** 53210 -0.80: ************************************************************************ 58127 -0.70: ****************************************************************************** 62438 -0.60: *********************************************************************************** 66650 -0.50: **************************************************************************************** 70298 -0.40: ******************************************************************************************** 73739 -0.30: *********************************************************************************************** 75831 -0.20: ************************************************************************************************** 78222 -0.10: *************************************************************************************************** 79412 0.00: **************************************************************************************************** 79801 0.10: *************************************************************************************************** 79255 0.20: ************************************************************************************************* 78163 0.30: ************************************************************************************************ 76667 0.40: ******************************************************************************************** 73554 0.50: **************************************************************************************** 70391 0.60: *********************************************************************************** 66566 0.70: ****************************************************************************** 62857 0.80: ************************************************************************ 57962 0.90: ****************************************************************** 53407 1.00: ************************************************************ 48565 1.10: ****************************************************** 43496 1.20: ************************************************ 38799 1.30: ****************************************** 34156 1.40: ************************************* 29713 1.50: ******************************** 25946 1.60: *************************** 22264 1.70: *********************** 18843 1.80: ******************* 15780 1.90: **************** 13151 2.00: ************* 10905 2.10: ********** 8690 2.20: ******** 7102 2.30: ******* 5693 2.40: ***** 4459 2.50: **** 3550 2.60: *** 2603 2.70: ** 2155 2.80: ** 1619 2.90: * 1121 3.00: * 914 3.10: 607 3.20: 478 3.30: 349 3.40: 216 3.50: 170 3.60: 113 3.70: 79 3.80: 58 3.90: 48 4.00: 33 4.10: 20 4.20: 9 4.30: 8 4.40: 7 4.50: 3 4.60: 3 4.70: 0 4.80: 1 4.90: 1
J
Solution
runif01=: ?@$ 0: NB. random uniform number generator
rnorm01=. (2 o. 2p1 * runif01) * [: %: _2 * ^.@runif01 NB. random normal number generator (Box-Muller)
mean=: +/ % # NB. mean
stddev=: (<:@# %~ +/)&.:*:@(- mean) NB. standard deviation
histogram=: <:@(#/.~)@(i.@#@[ , I.)
Example Usage
DataSet=: rnorm01 1e5
(mean , stddev) DataSet
0.000781667 1.00154
require 'plot'
plot (5 %~ i: 25) ([;histogram) DataSet
Java
import static java.lang.Math.*;
import static java.util.Arrays.stream;
import java.util.Locale;
import java.util.function.DoubleSupplier;
import static java.util.stream.Collectors.joining;
import java.util.stream.DoubleStream;
import static java.util.stream.IntStream.range;
public class Test implements DoubleSupplier {
private double mu, sigma;
private double[] state = new double[2];
private int index = state.length;
Test(double m, double s) {
mu = m;
sigma = s;
}
static double[] meanStdDev(double[] numbers) {
if (numbers.length == 0)
return new double[]{0.0, 0.0};
double sx = 0.0, sxx = 0.0;
long n = 0;
for (double x : numbers) {
sx += x;
sxx += pow(x, 2);
n++;
}
return new double[]{sx / n, pow((n * sxx - pow(sx, 2)), 0.5) / n};
}
static String replicate(int n, String s) {
return range(0, n + 1).mapToObj(i -> s).collect(joining());
}
static void showHistogram01(double[] numbers) {
final int maxWidth = 50;
long[] bins = new long[10];
for (double x : numbers)
bins[(int) (x * bins.length)]++;
double maxFreq = stream(bins).max().getAsLong();
for (int i = 0; i < bins.length; i++)
System.out.printf(" %3.1f: %s%n", i / (double) bins.length,
replicate((int) (bins[i] / maxFreq * maxWidth), "*"));
System.out.println();
}
@Override
public double getAsDouble() {
index++;
if (index >= state.length) {
double r = sqrt(-2 * log(random())) * sigma;
double x = 2 * PI * random();
state = new double[]{mu + r * sin(x), mu + r * cos(x)};
index = 0;
}
return state[index];
}
public static void main(String[] args) {
Locale.setDefault(Locale.US);
double[] data = DoubleStream.generate(new Test(0.0, 0.5)).limit(100_000)
.toArray();
double[] res = meanStdDev(data);
System.out.printf("Mean: %8.6f, SD: %8.6f%n", res[0], res[1]);
showHistogram01(stream(data).map(a -> max(0.0, min(0.9999, a / 3 + 0.5)))
.toArray());
}
}
Mean: -0.001870, SD: 0.500539 0.0: ** 0.1: ******* 0.2: ****************** 0.3: ************************************ 0.4: *************************************************** 0.5: ************************************************** 0.6: *********************************** 0.7: ****************** 0.8: ******* 0.9: **
jq
Adapted from Wren
Works with gojq, the Go implementation of jq (*)
Since jq does not have a built-in PRNG, this entry uses an external source for entropy. For the sake of illustration, we will use /dev/urandom as follows:
cat /dev/urandom | tr -cd '0-9' | fold -w 10 | jq -nRr -f normal-distribution.jq
To save space, the function that generates the sample does not retain the observations, and for simplicity, computes the sum of squared observations on the assumption that overflow will not be an an issue, which is reasonable as jq arithmetic uses IEEE 754 64-bit numbers.
(*) gojq requires an enormous amount of memory to complete the task for N=100,000, and takes about 20 times longer.
Preliminaries
# Pretty print a number to facilitate alignment of the decimal point.
# Input: a number without an exponent
# Output: a string holding the reformatted number so that there are at least `left` characters
# to the left of the decimal point, and exactly `right` characters to its right.
# Spaces are used for padding on the left, and zeros for padding on the right.
# No left-truncation occurs, so `left` can be specified as 0 to prevent left-padding.
def pp(left; right):
def lpad: tostring | (left - length) as $l | (" " * $l)[:$l] + .;
def rpad:
if (right > length) then . + ((right - length) * "0")
else .[:right]
end;
tostring as $s
| $s
| index(".") as $ix
| ((if $ix then $s[0:$ix] else $s end) | lpad) + "." +
(if $ix then $s[$ix+1:] | .[:right] else "" end | rpad) ;
def sigma( stream ): reduce stream as $x (0; . + $x);
# Input: {n, sum, ss}
# Output: augmented object with {mean, variance}
def sample_mean_and_variance:
.mean = (.sum/.n)
| .variance = ((.ss / .n) - .mean*.mean);
The Task
# Task parameters
def parameters: {
N: 100000,
NUM_BINS: 12,
HIST_CHAR: "■",
HIST_CHAR_ALT: "-",
HIST_CHAR_SIZE: null, # null means compute dynamically
binSize: 0.1,
mu: 0.5,
sigma: 0.25 }
| .bins = [range(0; .NUM_BINS) | 0] ;
# input: an array of two iid rvs on [0,1]
# output: [z0, z1] as per the Box-Muller method -- see
# https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
def normal(mu; sigma):
def pi: (1|atan) * 4;
. as [$u1, $u2]
| pi as $pi
| (sigma * ((-2 * ($u1|log))|sqrt)) as $mag
| [ $mag * ((2 * $pi * $u2)|cos) + mu,
$mag * ((2 * $pi * $u2)|sin) + mu ] ;
# Generate a random sample as specified by ., the task object (see `parameters`).
# Output: updated task object with sample statistics and .bins for creating a histogram.
# Each call to `input` should yield a string of random decimal digits
# such that the ensemble of ("0." + input | tonumber) can be considered to be iid rv on [0,1].
def generate:
# uniformly distributed random variable on [0,1]:
def udrv: "0." + input | tonumber;
# Maybe compute the bucket size:
(.HIST_CHAR_SIZE = (.HIST_CHAR_SIZE // (.N / (.NUM_BINS * 20) | ceil))) as $p
| reduce range(0; $p.N/2) as $i ($p;
([udrv, udrv] | normal($p.mu; $p.sigma)) as $rns
| reduce (0,1) as $j (.;
$rns[$j] as $rn
| .n += 1
| .sum += $rn
| .ss += ($rn*$rn)
| (if $rn < 0 then 0
elif $rn >= 1 then ($p.NUM_BINS - 1)
else ($rn/.binSize)|floor + 1
end ) as $bn
| .bins[$bn] += 1
# to retain the observations: .samples[$i*2 + $j] = $rn
)) ;
# Input: an object with
# {NUM_BINS, HIST_CHAR_SIZE, HIST_CHAR, HIST_CHAR_ALT, binSize, bins}
# Output: a stream of strings
def histogram:
def tidy: pp(2;1);
range(0; .NUM_BINS) as $i
| ((.bins[$i] / .HIST_CHAR_SIZE)|floor) as $bs
| (if $i == 0 or $i == .NUM_BINS -1
then .HIST_CHAR_ALT else .HIST_CHAR end) as $char
| (if $bs == 0 then "" else $char * $bs end) as $hist
| if $i == 0
then " -∞ ..< 0.0 \($hist)" # .bins[0]
elif ($i < .NUM_BINS - 1)
then "\(.binSize * ($i-1) | tidy) ..<\(.binSize * $i|tidy) \($hist)" # .bins[$i]]
else " 1.0 .. +∞ \($hist)" # .bins[.NUM_BINS - 1]
end;
def task:
parameters
| generate
| sample_mean_and_variance
| (if .HIST_CHAR_SIZE == 1 then "" else "s" end) as $plural
| "Summary statistics for \(.N) observations from N(\(.mu), \(.sigma)):",
" mean: \(.mean | pp(2;4))",
" variance: \(.variance | pp(2;4))",
" unadjusted stddev: \(.variance | sqrt | pp(2;4))",
" Range Number of observations (each \(.HIST_CHAR) represents \(.HIST_CHAR_SIZE) observation\($plural))",
histogram ;
task
- Output:
Summary statistics for 100000 observations from N(0.5, 0.25): mean: 0.5001 variance: 0.0622 unadjusted stddev: 0.2495 Range Number of observations (each ■ represents 417 observations) -∞ ..< 0.0 ----- 0.0 ..< 0.1 ■■■■■■■■ 0.1 ..< 0.2 ■■■■■■■■■■■■■■ 0.2 ..< 0.3 ■■■■■■■■■■■■■■■■■■■■■■■ 0.3 ..< 0.4 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.4 ..< 0.5 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.5 ..< 0.6 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.6 ..< 0.7 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.7 ..< 0.8 ■■■■■■■■■■■■■■■■■■■■■■■ 0.8 ..< 0.9 ■■■■■■■■■■■■■■ 0.9 ..< 1.0 ■■■■■■■■ 1.0 .. +∞ ------
Julia
Julia has the builtin package "Distributions" to generate random numbers from a standard distribution (Normal, Chisq etc.).
using Printf, Distributions, Gadfly
data = rand(Normal(0, 1), 1000)
@printf("N = %i\n", length(data))
@printf("μ = %2.2f\tσ = %2.2f\n", mean(data), std(data))
@printf("range = (%2.2f, %2.2f\n)", minimum(data), maximum(data))
h = plot(x=data, Geom.histogram)
draw(PNG("norm_hist.png", 10cm, 10cm), h)
- Output:
N = 1000 μ = 0.02 σ = 0.97 range = (-2.76, 3.42)
Kotlin
// version 1.1.2
val rand = java.util.Random()
fun normalStats(sampleSize: Int) {
if (sampleSize < 1) return
val r = DoubleArray(sampleSize)
val h = IntArray(12) // all zero by default
/*
Generate 'sampleSize' normally distributed random numbers with mean 0.5 and SD 0.25
and calculate in which box they will fall when drawing the histogram
*/
for (i in 0 until sampleSize) {
r[i] = 0.5 + rand.nextGaussian() / 4.0
when {
r[i] < 0.0 -> h[0]++
r[i] >= 1.0 -> h[11]++
else -> h[1 + (r[i] * 10).toInt()]++
}
}
// adjust one of the h[] values if necessary to ensure they sum to sampleSize
val adj = sampleSize - h.sum()
if (adj != 0) {
for (i in 0..11) {
h[i] += adj
if (h[i] >= 0) break
h[i] -= adj
}
}
val mean = r.average()
val sd = Math.sqrt(r.map { (it - mean) * (it - mean) }.average())
// Draw a histogram of the data with interval 0.1
var numStars: Int
// If sample size > 300 then normalize histogram to 300
val scale = if (sampleSize <= 300) 1.0 else 300.0 / sampleSize
println("Sample size $sampleSize\n")
println(" Mean ${"%1.6f".format(mean)} SD ${"%1.6f".format(sd)}\n")
for (i in 0..11) {
when (i) {
0 -> print("< 0.00 : ")
11 -> print(">=1.00 : ")
else -> print(" %1.2f : ".format(i / 10.0))
}
print("%5d ".format(h[i]))
numStars = (h[i] * scale + 0.5).toInt()
println("*".repeat(numStars))
}
println()
}
fun main(args: Array<String>) {
val sampleSizes = intArrayOf(100, 1_000, 10_000, 100_000)
for (sampleSize in sampleSizes) normalStats(sampleSize)
}
- Output:
Sample size 100 Mean 0.525211 SD 0.266316 < 0.00 : 3 *** 0.10 : 1 * 0.20 : 3 *** 0.30 : 11 *********** 0.40 : 14 ************** 0.50 : 13 ************* 0.60 : 15 *************** 0.70 : 13 ************* 0.80 : 10 ********** 0.90 : 11 *********** 1.00 : 4 **** >=1.00 : 2 ** Sample size 1000 Mean 0.500948 SD 0.255757 < 0.00 : 29 ********* 0.10 : 35 *********** 0.20 : 70 ********************* 0.30 : 71 ********************* 0.40 : 138 ***************************************** 0.50 : 139 ****************************************** 0.60 : 168 ************************************************** 0.70 : 123 ************************************* 0.80 : 110 ********************************* 0.90 : 62 ******************* 1.00 : 32 ********** >=1.00 : 23 ******* Sample size 10000 Mean 0.501376 SD 0.248317 < 0.00 : 240 ******* 0.10 : 305 ********* 0.20 : 617 ******************* 0.30 : 927 **************************** 0.40 : 1291 *************************************** 0.50 : 1554 *********************************************** 0.60 : 1609 ************************************************ 0.70 : 1319 **************************************** 0.80 : 983 ***************************** 0.90 : 609 ****************** 1.00 : 324 ********** >=1.00 : 222 ******* Sample size 100000 Mean 0.499427 SD 0.250533 < 0.00 : 2341 ******* 0.10 : 3246 ********** 0.20 : 6005 ****************** 0.30 : 9718 ***************************** 0.40 : 13247 **************************************** 0.50 : 15595 *********************************************** 0.60 : 15271 ********************************************** 0.70 : 13405 **************************************** 0.80 : 9653 ***************************** 0.90 : 5990 ****************** 1.00 : 3257 ********** >=1.00 : 2272 *******
Lasso
define stat1(a) => {
if(#a->size) => {
local(mean = (with n in #a sum #n) / #a->size)
local(sdev = math_pow(((with n in #a sum Math_Pow((#n - #mean),2)) / #a->size),0.5))
return (:#sdev, #mean)
else
return (:0,0)
}
}
define stat2(a) => {
if(#a->size) => {
local(sx = 0, sxx = 0)
with x in #a do => {
#sx += #x
#sxx += #x*#x
}
local(sdev = math_pow((#a->size * #sxx - #sx * #sx),0.5) / #a->size)
return (:#sdev, #sx / #a->size)
else
return (:0,0)
}
}
define histogram(a) => {
local(
out = '\r',
h = array(0,0,0,0,0,0,0,0,0,0,0),
maxwidth = 50,
sc = 0
)
with n in #a do => {
if((#n * 10) <= 0) => {
#h->get(1) += 1
else((#n * 10) >= 10)
#h->get(#h->size) += 1
else
#h->get(integer(decimal(#n)*10)+1) += 1
}
}
local(mx = decimal(with n in #h max #n))
with i in #h do => {
#out->append((#sc/10.0)->asString(-precision=1)+': '+('+' * integer(#i / #mx * #maxwidth))+'\r')
#sc++
}
return #out
}
define normalDist(mean,sdev) => {
// Uses Box-Muller transform
return ((-2 * decimal_random->log)->sqrt * (2 * pi * decimal_random)->cos) * #sdev + #mean
}
with scale in array(100,1000,10000) do => {^
local(n = array)
loop(#scale) => { #n->insert(normalDist(0.5, 0.2)) }
local(sdev1,mean1) = stat1(#n)
local(sdev2,mean2) = stat2(#n)
#scale' numbers:\r'
'Naive method: sd: '+#sdev1+', mean: '+#mean1+'\r'
'Second method: sd: '+#sdev2+', mean: '+#mean2+'\r'
histogram(#n)
'\r\r'
^}
- Output:
100 numbers: Naive method: sd: 0.199518, mean: 0.506059 Second method: sd: 0.199518, mean: 0.506059 0.0: ++ 0.1: ++++ 0.2: +++++++++++++++++ 0.3: ++++++++++++++++++++++ 0.4: ++++++++++++++++++++++++++++++++++++++++++++++++++ 0.5: +++++++++++++++++++++++++++++++++++++++ 0.6: +++++++++++++++++++++++++++++++++ 0.7: ++++++++++++++++++++++++ 0.8: ++++++++++++++++++++ 0.9: ++++ 1.0: ++ 1000 numbers: Naive method: sd: 0.199653, mean: 0.504046 Second method: sd: 0.199653, mean: 0.504046 0.0: +++ 0.1: ++++++ 0.2: ++++++++++++++++ 0.3: ++++++++++++++++++++++++++++++ 0.4: +++++++++++++++++++++++++++++++++++++++++++++++ 0.5: ++++++++++++++++++++++++++++++++++++++++++++++++++ 0.6: ++++++++++++++++++++++++++++++++++++++++++++++ 0.7: +++++++++++++++++++++++++ 0.8: +++++++++++++++++++ 0.9: +++++++ 1.0: ++++ 10000 numbers: Naive method: sd: 0.202354, mean: 0.502519 Second method: sd: 0.202354, mean: 0.502519 0.0: +++ 0.1: +++++++ 0.2: +++++++++++++++ 0.3: +++++++++++++++++++++++++++++ 0.4: ++++++++++++++++++++++++++++++++++++++++++ 0.5: ++++++++++++++++++++++++++++++++++++++++++++++++++ 0.6: +++++++++++++++++++++++++++++++++++++++++++ 0.7: ++++++++++++++++++++++++++++++ 0.8: ++++++++++++++++ 0.9: +++++++ 1.0: ++++
Liberty BASIC
Uses LB Statistics/Basic
call sample 100000
end
sub sample n
dim dat( n)
for i =1 to n
dat( i) =normalDist( 1, 0.2)
next i
'// show mean, standard deviation. Find max, min.
mx =-1000
mn = 1000
sum =0
sSq =0
for i =1 to n
d =dat( i)
mx =max( mx, d)
mn =min( mn, d)
sum =sum +d
sSq =sSq +d^2
next i
print n; " data terms used."
mean =sum / n
print "Largest term was "; mx; " & smallest was "; mn
range =mx -mn
print "Mean ="; mean
print "Stddev ="; ( sSq /n -mean^2)^0.5
'// show histogram
nBins =50
dim bins( nBins)
for i =1 to n
z =int( ( dat( i) -mn) /range *nBins)
bins( z) =bins( z) +1
next i
for b =0 to nBins -1
for j =1 to int( nBins *bins( b)) /n *30)
print "#";
next j
print
next b
print
end sub
function normalDist( m, s) ' Box Muller method
u =rnd( 1)
v =rnd( 1)
normalDist =( -2 *log( u))^0.5 *cos( 2 *3.14159265 *v)
end function
100000 data terms used. Largest term was 4.12950792 & smallest was -4.37934139 Mean =-0.26785425e-2 Stddev =1.00097669
# ## ### ##### ######## ############ ################ ######################## ############################## ###################################### ############################################## ######################################################## ################################################################### ############################################################################## ####################################################################################### ################################################################################################ #################################################################################################### ######################################################################################################## ##################################################################################################### ############################################################################################## ######################################################################################### ################################################################################## ######################################################################### ############################################################## #################################################### ########################################## ################################# ########################## ################## ############# ######### ###### #### ## # #
Lua
Lua provides math.random() to generate uniformly distributed random numbers. The function gaussian() shown here uses math.random() to generate normally distributed random numbers with given mean and variance.
function gaussian (mean, variance)
return math.sqrt(-2 * variance * math.log(math.random())) *
math.cos(2 * math.pi * math.random()) + mean
end
function mean (t)
local sum = 0
for k, v in pairs(t) do
sum = sum + v
end
return sum / #t
end
function std (t)
local squares, avg = 0, mean(t)
for k, v in pairs(t) do
squares = squares + ((avg - v) ^ 2)
end
local variance = squares / #t
return math.sqrt(variance)
end
function showHistogram (t)
local lo = math.ceil(math.min(unpack(t)))
local hi = math.floor(math.max(unpack(t)))
local hist, barScale = {}, 200
for i = lo, hi do
hist[i] = 0
for k, v in pairs(t) do
if math.ceil(v - 0.5) == i then
hist[i] = hist[i] + 1
end
end
io.write(i .. "\t" .. string.rep('=', hist[i] / #t * barScale))
print(" " .. hist[i])
end
end
math.randomseed(os.time())
local t, average, variance = {}, 50, 10
for i = 1, 1000 do
table.insert(t, gaussian(average, variance))
end
print("Mean:", mean(t) .. ", expected " .. average)
print("StdDev:", std(t) .. ", expected " .. math.sqrt(variance) .. "\n")
showHistogram(t)
- Output:
Mean: 50.008328894275, expected 50 StdDev: 3.2374717425824, expected 3.1622776601684 41 3 42 = 8 43 == 11 44 ==== 22 45 ======= 38 46 ============ 60 47 ============== 73 48 ================== 92 49 ======================= 118 50 =========================== 136 51 ========================= 128 52 ================= 89 53 ================= 89 54 =========== 56 55 ======= 37 56 === 18 57 = 7 58 = 5 59 = 6 60 2
Maple
Maple has a built-in for sampling directly from Normal distributions:
with(Statistics):
n := 100000:
X := Sample( Normal(0,1), n );
Mean( X );
StandardDeviation( X );
Histogram( X );
Mathematica/Wolfram Language
x:= RandomReal[1]
SampleNormal[n_] := (Print[#//Length, " numbers, Mean : ", #//Mean, ", StandardDeviation : ", #//StandardDeviation];
Histogram[#, BarOrigin -> Left,Axes -> False])& [(Table[(-2*Log[x])^0.5*Cos[2*Pi*x], {n} ]]
Invocation:
SampleNormal[ 10000 ]
->10000 numbers, Mean : -0.0122308, StandardDeviation : 1.00646
MATLAB / Octave
N = 100000;
x = randn(N,1);
mean(x)
std(x)
[nn,xx] = hist(x,100);
bar(xx,nn);
Nim
In module “random”, Nim provides two procedures named gauss
to generate random values following normal distribution and following Gauss distribution with given mean and standard deviation.
Here is a way to generate random values following normal distribution from random values following uniform distribution. It uses the Basic form of the Box-Muller transform.
import math, random, sequtils, stats, strformat, strutils
proc drawHistogram(ns: seq[float]) =
# Distribute values in bins.
const NBins = 50
var minval = min(ns)
var maxval = max(ns)
var h = newSeq[int](NBins + 1)
for n in ns:
let pos = ((n - minval) * NBins / (maxval - minval)).toInt
inc h[pos]
# Eliminate extremes values.
const MaxWidth = 50
let mx = max(h)
var first = 0
while (h[first] / mx * MaxWidth).toInt == 0: inc first
var last = h.high
while (h[last] / mx * MaxWidth).toInt == 0: dec last
# Draw the histogram.
echo ""
for n in first..last:
echo repeat('+', (h[n] / mx * MaxWidth).toInt)
echo ""
const N = 100_000
randomize()
let u1, u2 = newSeqWith(N, rand(1.0))
var z = newSeq[float](N)
for i in 0..<N:
z[i] = sqrt(-2 * ln(u1[i])) * cos(2 * PI * u2[i])
echo &"μ = {z.mean:.12f} σ = {z.standardDeviation:.12f}"
z.drawHistogram()
- Output:
μ = -0.001105836229 σ = 0.999906544722 + + ++ +++ +++++ +++++++ +++++++++ ++++++++++++ ++++++++++++++++ +++++++++++++++++++++ ++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++ +++++++++++++++++++++++++ ++++++++++++++++++++ ++++++++++++++++ ++++++++++++ +++++++++ ++++++ +++++ +++ ++ + +
PARI/GP
rnormal()={
my(u1=random(1.),u2=random(1.));
sqrt(-2*log(u1))*cos(2*Pi*u2)
\\ Could easily be extended with a second normal at very little cost.
};
mean(v)={
sum(i=1,#v,v[i])/#v
};
stdev(v,mu="")={
if(mu=="",mu=mean(v));
sqrt(sum(i=1,#v,(v[i]-mu)^2))/#v
};
histogram(v,bins=16,low=0,high=1)={
my(u=vector(bins),width=(high-low)/bins);
for(i=1,#v,u[(v[i]-low)\width+1]++);
u
};
show(n)={
my(v=vector(n,i,rnormal()),m=mean(v),s=stdev(v,m),h,sz=ceil(n/300));
h=histogram(v,,vecmin(v)-.1,vecmax(v)+.1);
for(i=1,#h,for(j=1,h[i]\sz,print1("#"));print());
};
show(10^4)
For versions before 2.4.3, define
rreal()={
my(pr=32*ceil(default(realprecision)*log(10)/log(4294967296))); \\ Current precision
random(2^pr)*1.>>pr
};
and use rreal()
in place of random(1.)
.
A PARI implementation:
GEN
rnormal(long prec)
{
pari_sp ltop = avma;
GEN u1, u2, left, right, ret;
u1 = randomr(prec);
u2 = randomr(prec);
left = sqrtr_abs(shiftr(mplog(u1), 1));
right = mpcos(mulrr(shiftr(mppi(prec), 1), u2));
ret = mulrr(left, right);
ret = gerepileupto(ltop, ret);
return ret;
}
Use mpsincos
and caching to generate two values at nearly the same cost.
Pascal
//not neccessary include unit math if using function rnorm
got some trouble with math.randg needs this call randg(cMean,cMean*cStdDiv), whereas randg(cMean,cStdDiv) to get the same results??
From Free Pascal Docs unit math
Program Example40;
{$IFDEF FPC}
{$MOde objFPC}
{$ENDIF}
{ Program to demonstrate the randg function. }
Uses Math;
type
tTestData = extended;//because of math.randg
ttstfunc = function (mean, sd: tTestData): tTestData;
tExArray = Array of tTestData;
tSolution = record
SolExArr : tExArray;
SollowVal,
SolHighVal,
SolMean,
SolStdDiv : tTestData;
SolSmpCnt : LongInt;
end;
function getSol(genFunc:ttstfunc;Mean,StdDiv: tTestData;smpCnt: LongInt): tSolution;
var
GenValue,
sumValue,
sumsqrVal : extended;
Begin
with result do
Begin
SolSmpCnt := smpCnt;
SolMean := 0;
SolStdDiv := 0;
SolLowVal := Mean+50* StdDiv;
SolHighVal := Mean-50* StdDiv;
setlength(SolExArr,smpCnt);
if smpCnt <= 0 then
EXIT;
sumValue := 0;
sumsqrVal := 0;
repeat
GenValue := genFunc(Mean,StdDiv);
sumValue := sumvalue+GenValue;
sumsqrVal := sumsqrVal+sqr(GenValue);
IF GenValue < SollowVal then
SollowVal:= GenValue
else
IF GenValue > SolHighVal then
SolHighVal := GenValue;
dec(smpCnt);
SolExArr[smpCnt] := GenValue;
until smpCnt<= 0;
SolMean := sumValue/SolSmpCnt;
SolStdDiv := sqrt(sumsqrVal/SolSmpCnt-sqr(SolMean));
end;
end;
//http://wiki.freepascal.org/Generating_Random_Numbers#Normal_.28Gaussian.29_Distribution
function rnorm (mean, sd: tTestData): tTestData;
{Calculates Gaussian random numbers according to the Box-Müller approach}
var
u1, u2: extended;
begin
u1 := random;
u2 := random;
rnorm := mean * abs(1 + sqrt(-2 * (ln(u1))) * cos(2 * pi * u2) * sd);
end;
procedure Histo(const sol:TSolution;Colcnt,ColLen :LongInt);
var
CntHisto : array of integer;
LoLmt,HiLmt,span : tTestData;
i, j,cnt,maxCnt: LongInt;
sCross : Ansistring;
Begin
setlength(CntHisto,Colcnt);
with Sol do
Begin
span := solHighVal-solLowVal;
LoLmt := solLowVal;
writeln('Count: ',SolSmpCnt:10,' Mean ',SolMean:10:6,' StdDiv ',SolStdDIv:10:6);
writeln('span : ',span:10:5,' Low ',solLowVal:10:6,' high ',solHighVal:10:6);
end;
maxCnt := 0;
For j := 0 to Colcnt-1 do
Begin
HiLmt:= LoLmt+span/Colcnt;
cnt := 0;
with sol do
For i := 0 to High(SolExArr) do
IF (HiLmt > SolExArr[i]) AND (SolExArr[i]>= LoLmt) then
inc(cnt);
CntHisto[j] := cnt;
IF maxCnt < cnt then
maxCnt := cnt;
LoLmt:= HiLmt;
end;
inc(CntHisto[Colcnt]); // for HiLmt itself
writeln;
LoLmt := sol.solLowVal;
For i := 0 to Colcnt-1 do
Begin
Writeln(LoLmt:8:4,': ');
cnt:= Round(CntHisto[i]*ColLen/maxCnt);
setlength(sCross,cnt+3);
fillChar(sCross[1],3,' ');
fillChar(sCross[4],cnt,'#');
writeln(CntHisto[i]:10,sCross);
LoLmt := LoLmt+span/Colcnt;
end;
Writeln(sol.solHighVal:8:4,': ');
end;
const
cHistCnt = 11;
cColLen = 65;
cStdDiv = 0.25;
cMean = 20*cStdDiv;
var
mySol : tSolution;
begin
Randomize;
// test of randg of unit math
Writeln('function randg');
mySol := getSol(@randg,cMean,cMean*cStdDiv,100000);
Histo(mySol,cHistCnt,cColLen);
writeln;
// test of rnorm from wiki
Writeln('function rnorm');
mySol := getSol(@rnorm,cMean,cStdDiv,1000000);
Histo(mySol,cHistCnt,cColLen);
end.
- Output:
function randg Count: 100000 Mean 5.000326 StdDiv 1.250027 span : 10.65123 Low -0.333310 high 10.317922
-0.3333: 25 0.6350: 287 # 1.6033: 2291 ##### 2.5716: 9531 ##################### 3.5399: 22608 ################################################# 4.5082: 29953 ################################################################# 5.4765: 22917 ################################################## 6.4447: 9716 ##################### 7.4130: 2352 ##### 8.3813: 295 # 9.3496: 24 10.3179:function rnorm Count: 1000000 Mean 4.998391 StdDiv 1.251103 span : 11.08994 Low 0.001521 high 11.091461
0.0015: 704 1.0097: 7797 ## 2.0179: 49235 ########### 3.0261: 162761 #################################### 4.0342: 293242 ################################################################# 5.0424: 285818 ############################################################### 6.0506: 150781 ################################# 7.0588: 42641 ######### 8.0669: 6467 # 9.0751: 528 10.0833: 25 11.0915:
Perl
use constant pi => 3.14159265;
use List::Util qw(sum reduce min max);
sub normdist {
my($m, $sigma) = @_;
my $r = sqrt -2 * log rand;
my $theta = 2 * pi * rand;
$r * cos($theta) * $sigma + $m;
}
$size = 100000; $mean = 50; $stddev = 4;
push @dataset, normdist($mean,$stddev) for 1..$size;
my $m = sum(@dataset) / $size;
print "m = $m\n";
my $sigma = sqrt( (reduce { $a + $b **2 } 0,@dataset) / $size - $m**2 );
print "sigma = $sigma\n";
$hash{int $_}++ for @dataset;
my $scale = 180 * $stddev / $size;
my @subbar = < ⎸ ▏ ▎ ▍ ▌ ▋ ▊ ▉ █ >;
for $i (min(@dataset)..max(@dataset)) {
my $x = ($hash{$i} // 0) * $scale;
my $full = int $x;
my $part = 8 * ($x - $full);
my $t1 = '█' x $full;
my $t2 = $subbar[$part];
print "$i\t$t1$t2\n";
}
- Output:
32 ⎸ 33 ⎸ 34 ⎸ 35 ⎸ 36 ▎ 37 ▋ 38 █▏ 39 ██▍ 40 ████▍ 41 ███████▌ 42 ████████████⎸ 43 ███████████████████▏ 44 ████████████████████████████⎸ 45 ██████████████████████████████████████▎ 46 █████████████████████████████████████████████████▎ 47 ██████████████████████████████████████████████████████████▋ 48 ██████████████████████████████████████████████████████████████████▋ 49 ███████████████████████████████████████████████████████████████████████▍ 50 ██████████████████████████████████████████████████████████████████████▋ 51 ██████████████████████████████████████████████████████████████████▌ 52 ████████████████████████████████████████████████████████████▎ 53 ████████████████████████████████████████████████▏ 54 █████████████████████████████████████▊ 55 ███████████████████████████▍ 56 ███████████████████▊ 57 ████████████▌ 58 ███████▌ 59 ████▍ 60 ██▏ 61 █⎸ 62 ▌ 63 ▏ 64 ⎸ 65 ⎸ 66 ⎸
Phix
with javascript_semantics procedure sample(integer n) -- show mean, standard deviation. Find max, min. sequence dat = repeat(0,n) for i=1 to n do dat[i] = sqrt(-2*log(rnd()))*cos(2*PI*rnd()) end for printf(1,"%d data terms used.\n",{n}) atom mean = sum(dat)/n, mx = max(dat), mn = min(dat), range = mx-mn printf(1,"Largest term was %g & smallest was %g\n",{mx,mn}) printf(1,"Mean = %g\n",{mean}) printf(1,"Stddev = %g\n",sqrt(sum(sq_mul(dat,dat))/n-mean*mean)) -- show histogram integer nBins = 50 sequence bins = repeat(0,nBins+1) for i=1 to n do integer bdx = floor((dat[i]-mn)/range*nBins)+1 bins[bdx] += 1 end for for b=1 to nBins do puts(1,repeat('#',floor(nBins*bins[b]/n*30))&"\n") end for end procedure sample(100000)
- Output:
100000 data terms used. Largest term was 4.30779 & smallest was -4.11902 Mean = -0.00252597 Stddev = 1.00067 # ## #### ###### ########## ############# ################## ######################## ################################# ######################################## #################################################### ############################################################# ###################################################################### ############################################################################### ####################################################################################### ############################################################################################### ################################################################################################# ##################################################################################################### ################################################################################################### ################################################################################################ ######################################################################################## ############################################################################### ####################################################################### ############################################################ ################################################# ####################################### ############################## ######################### ################ ############ ######### ###### #### ## #
with javascript_semantics function gaussian(atom mean, atom variance) return sqrt(-2 * variance * log(rnd())) * cos(2 * variance * PI * rnd()) + mean end function function mean(sequence t) return sum(t)/length(t) end function function std(sequence t) atom squares = 0, avg = mean(t) for i=1 to length(t) do squares += power(avg-t[i],2) end for atom variance = squares/length(t) return sqrt(variance) end function procedure showHistogram(sequence t) for i=ceil(min(t)) to floor(max(t)) do integer n = 0 for k=1 to length(t) do n += ceil(t[k]-0.5)=i end for integer l = floor(n/length(t)*200) printf(1,"%d %s %d\n",{i,repeat('=',l),n}) end for end procedure sequence t = repeat(0,100000) integer avg = 50, variance = 10 for i=1 to length(t) do t[i] = gaussian(avg, variance) end for printf(1,"Mean: %g, expected %g\n",{mean(t),avg}) printf(1,"StdDev: %g, expected %g\n",{std(t),sqrt(variance)}) showHistogram(t)
- Output:
Mean: 50.0041, expected 50 StdDev: 3.1673, expected 3.16228 37 2 38 7 39 30 40 92 41 220 42 = 523 43 == 1098 44 ==== 2140 45 ======= 3690 46 =========== 5753 47 =============== 7906 48 ==================== 10299 49 ======================= 11813 50 ========================= 12555 51 ======================= 11934 52 ==================== 10327 53 ================ 8099 54 =========== 5733 55 ======= 3684 56 ==== 2126 57 == 1098 58 487 59 226 60 106 61 36 62 9 63 7
PureBasic
Procedure.f randomf(resolution = 2147483647)
ProcedureReturn Random(resolution) / resolution
EndProcedure
Procedure.f normalDist() ;Box Muller method
ProcedureReturn Sqr(-2 * Log(randomf())) * Cos(2 * #PI * randomf())
EndProcedure
Procedure sample(n, nBins = 50)
Protected i, maxBinValue, binNumber
Protected.f d, mean, sum, sumSq, mx, mn, range
Dim dat.f(n)
For i = 1 To n
dat(i) = normalDist()
Next
;show mean, standard deviation, find max & min.
mx = -1000
mn = 1000
sum = 0
sumSq = 0
For i = 1 To n
d = dat(i)
If d > mx: mx = d: EndIf
If d < mn: mn = d: EndIf
sum + d
sumSq + d * d
Next
PrintN(Str(n) + " data terms used.")
PrintN("Largest term was " + StrF(mx) + " & smallest was " + StrF(mn))
mean = sum / n
PrintN("Mean = " + StrF(mean))
PrintN("Stddev = " + StrF((sumSq / n) - Sqr(mean * mean)))
;show histogram
range = mx - mn
Dim bins(nBins)
For i = 1 To n
binNumber = Int(nBins * (dat(i) - mn) / range)
bins(binNumber) + 1
Next
maxBinValue = 1
For i = 0 To nBins
If bins(i) > maxBinValue
maxBinValue = bins(i)
EndIf
Next
#normalizedMaxValue = 70
For binNumber = 0 To nBins
tickMarks = Round(bins(binNumber) * #normalizedMaxValue / maxBinValue, #PB_Round_Nearest)
PrintN(ReplaceString(Space(tickMarks), " ", "#"))
Next
PrintN("")
EndProcedure
If OpenConsole()
sample(100000)
Print(#CRLF$ + #CRLF$ + "Press ENTER to exit"): Input()
CloseConsole()
EndIf
Sample output:
100000 data terms used. Largest term was 4.5352029800 & smallest was -4.5405135155 Mean = 0.0012346541 Stddev = 0.9959455132 # ### ###### ########## ################## ############################ ######################################### ##################################################### ################################################################ ###################################################################### ###################################################################### ################################################################ ##################################################### ######################################### ############################# ################## ########## ###### ### #
Python
This uses the external matplotlib package as well as the built-in standardlib function random.gauss.
from __future__ import division
import matplotlib.pyplot as plt
import random
mean, stddev, size = 50, 4, 100000
data = [random.gauss(mean, stddev) for c in range(size)]
mn = sum(data) / size
sd = (sum(x*x for x in data) / size
- (sum(data) / size) ** 2) ** 0.5
print("Sample mean = %g; Stddev = %g; max = %g; min = %g for %i values"
% (mn, sd, max(data), min(data), size))
plt.hist(data,bins=50)
- Output:
Sample mean = 49.9822; Stddev = 4.00938; max = 66.8091; min = 33.5283 for 100000 values
R
Generate normal random numbers from uniform random numbers, using the Box-Muller transform. Both x and y are normally distributed, and independent.
n <- 100000
u <- sqrt(-2*log(runif(n)))
v <- 2*pi*runif(n)
x <- u*cos(v)
y <- v*sin(v)
hist(x)
R can generate random normal distributed numbers using the rnorm function:
n <- 100000
x <- rnorm(n, mean=0, sd=1)
mean(x)
sd(x)
hist(x)
Racket
This shows how one would generate samples from a normal distribution, compute statistics and plot a histogram.
#lang racket
(require math (planet williams/science/histogram-with-graphics))
(define data (sample (normal-dist 50 4) 100000))
(displayln (~a "Mean:\t" (mean data)))
(displayln (~a "Stddev:\t" (stddev data)))
(displayln (~a "Max:\t" (apply max data)))
(displayln (~a "Min:\t" (apply min data)))
(define h (make-histogram-with-ranges-uniform 40 30 70))
(for ([x data]) (histogram-increment! h x))
(histogram-plot h "Normal distribution μ=50 σ=4")
The other part of the task was to produce normal distributed numbers from a unit distribution. The following code is an implementation of the polar method. It is a slightly modified version of code originally written by Sebastian Egner.
#lang racket
(require math)
(define random-normal
(let ([unit (uniform-dist)]
[next #f])
(λ (μ σ)
(if next
(begin0
(+ μ (* σ next))
(set! next #f))
(let loop ()
(let* ([v1 (- (* 2.0 (sample unit)) 1.0)]
[v2 (- (* 2.0 (sample unit)) 1.0)]
[s (+ (sqr v1) (sqr v2))])
(cond [(>= s 1) (loop)]
[else (define scale (sqrt (/ (* -2.0 (log s)) s)))
(set! next (* scale v2))
(+ μ (* σ scale v1))])))))))
Raku
(formerly Perl 6)
sub normdist ($m, $σ) {
my $r = sqrt -2 * log rand;
my $Θ = τ * rand;
$r * cos($Θ) * $σ + $m;
}
sub MAIN ($size = 100000, $mean = 50, $stddev = 4) {
my @dataset = normdist($mean,$stddev) xx $size;
my $m = [+](@dataset) / $size;
say (:$m);
my $σ = sqrt [+](@dataset X** 2) / $size - $m**2;
say (:$σ);
(my %hash){.round}++ for @dataset;
my $scale = 180 * $stddev / $size;
constant @subbar = < ⎸ ▏ ▎ ▍ ▌ ▋ ▊ ▉ █ >;
for %hash.keys».Int.minmax(+*) -> $i {
my $x = (%hash{$i} // 0) * $scale;
my $full = floor $x;
my $part = 8 * ($x - $full);
say $i, "\t", '█' x $full, @subbar[$part];
}
}
- Output:
"m" => 50.006107405837142e0 "σ" => 4.0814435639885254e0 33 ⎸ 34 ⎸ 35 ⎸ 36 ▏ 37 ▎ 38 ▊ 39 █▋ 40 ███⎸ 41 █████▊ 42 ██████████⎸ 43 ███████████████▋ 44 ███████████████████████▏ 45 ████████████████████████████████▌ 46 ███████████████████████████████████████████▍ 47 ██████████████████████████████████████████████████████▏ 48 ███████████████████████████████████████████████████████████████▏ 49 █████████████████████████████████████████████████████████████████████▋ 50 ███████████████████████████████████████████████████████████████████████▊ 51 █████████████████████████████████████████████████████████████████████▌ 52 ███████████████████████████████████████████████████████████████⎸ 53 ██████████████████████████████████████████████████████▎ 54 ███████████████████████████████████████████⎸ 55 ████████████████████████████████▌ 56 ███████████████████████▍ 57 ███████████████▉ 58 █████████▉ 59 █████▍ 60 ███▍ 61 █▋ 62 ▊ 63 ▍ 64 ▏ 65 ⎸ 66 ⎸ 67 ⎸
REXX
The REXX language doesn't have any "higher math" BIF functions like SIN, COS, LN, LOG, SQRT, EXP, POW, etc,
so we hoi polloi programmers have to roll our own.
/*REXX program generates 10,000 normally distributed numbers (Gaussian distribution).*/
numeric digits 20 /*use twenty decimal digs for accuracy.*/
parse arg n seed . /*obtain optional arguments from the CL*/
if n=='' | n=="," then n= 10000 /*Not specified? Then use the default.*/
if datatype(seed, 'W') then call random ,,seed /*seed is for repeatable RANDOM numbers*/
call pi /*call subroutine to define pi constant*/
do g=1 for n; #.g= sqrt( -2 * ln( rand() ) ) * cos( 2 * pi * rand() )
end /*g*/ /* [↑] uniform random number ───► #.g */
s= 0
mn= #.1; mx= mn; noise= n * .0005 /*calculate the noise: 1/20th % of N.*/
ss= 0
do j=1 for n; _= #.j /*the sum, and the sum of squares. */
s= s + _; ss= ss + _ * _ /*the sum, and the sum of squares. */
mn= min(mn, _); mx= max(mx, _) /*find the minimum and the maximum. */
end /*j*/
!.= 0
say 'number of data points = ' fmt(n )
say ' minimum = ' fmt(mn )
say ' maximum = ' fmt(mx )
say ' arithmetic mean = ' fmt(s/n)
say ' standard deviation = ' fmt(sqrt( ss/n - (s/n) **2) )
?mn= !.1; ?mx= ?mn /*define minimum & maximum value so far*/
parse value scrSize() with sd sw . /*obtain the (true) screen size of term*/ /*◄──not all REXXes have this BIF*/
sdE= sd - 4 /*the effective (usable) screen depth. */
swE= sw - 1 /* " " " " width.*/
$= 1 / max(1, mx-mn) * sdE /*$ is used for scaling depth of histo*/
do i=1 for n; ?= trunc((#.i-mn) *$) /*calculate the relative line. */
!.?= !.? + 1 /*bump the counter. */
?mn= min(?mn, !.?) /*find the minimum. */
?mx= max(?mx, !.?) /* " " maximum. */
end /*i*/
f= swE/?mx /*limit graph to 1 full screen*/
do h=0 for sdE; _= !.h /*obtain a data point. */
if _>noise then say copies('─', trunc(_*f) ) /*display a bar of histogram. */
end /*h*/ /*[↑] use a hyphen for histo.*/
exit /*stick a fork in it, we're all done. */
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/
fmt: parse arg @; return left('', (@>=0) + 2 * datatype(@, 'W'))@ /*prepend a blank if #>=0, add 2 blanks if whole.*/
e: e = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535; return e
pi: pi= 3.1415926535897932384626433832795028841971693993751058209749445923078164062862; return pi
r2r: return arg(1) // (pi() * 2) /*normalize the given angle (in radians) to ±2pi.*/
rand: return random(1, 1e5) / 1e5 /*REXX generates uniform random postive integers.*/
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/
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; z= 1; _= 1
x= x*x; p= z; do k=2 by 2; _= -_ * x / (k*(k-1)); z= z + _; if z=p then leave; p= z; end; return z
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x; if x=0 then return 0; d= digits(); m.= 9; numeric digits; numeric form; h= d+6
parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g*.5'e'_%2; 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
This REXX program makes use of scrsize REXX program (or BIF) which is used to determine the screen size of the terminal (console); this is to aid in maximizing the width of the horizontal histogram.
The SCRSIZE.REX REXX program is included here ──► SCRSIZE.REX.
- output when using the default input:
(The output shown when the screen size is 62x140.)
The graph is shown at 3/4 size.
number of data points = 10000 minimum = -3.8181072371544448250 maximum = 3.5445917138265268562 arithmetic mean = -0.01406470979976873427 standard deviation = 0.99486092191249231518 ─ ─ ─── ──── ───── ───── ──────── ─────────── ────────────── ───────────────────── ────────────────────── ────────────────────────────────── ──────────────────────────────────────── ─────────────────────────────────────────────── ───────────────────────────────────────────────────── ───────────────────────────────────────────────────────────────────────── ───────────────────────────────────────────────────────────────── ───────────────────────────────────────────────────────────────────────────────────── ────────────────────────────────────────────────────────────────────────────────────────────────── ────────────────────────────────────────────────────────────────────────────────────────────── ──────────────────────────────────────────────────────────────────────────────────────────────────────────────── ──────────────────────────────────────────────────────────────────────────────────────────────────────────────── ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ────────────────────────────────────────────────────────────────────────────────────────────────────────── ────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ──────────────────────────────────────────────────────────────────────────────────────────────────────────────── ─────────────────────────────────────────────────────────────────────────────────────────────── ──────────────────────────────────────────────────────────────────────────── ─────────────────────────────────────────────────────────────────────── ─────────────────────────────────────────────────────────────────── ───────────────────────────────────────────────────── ───────────────────────────────────────────────── ───────────────────────────────── ────────────────────────────────── ─────────────────────── ────────────────────── ────────────────── ─────────── ────────── ────── ─── ──── ── ─
Ruby
The Implementation
# Class to implement a Normal distribution, generated from a Uniform distribution.
# Uses the Marsaglia polar method.
class NormalFromUniform
# Initialize an instance.
def initialize()
@next = nil
end
# Generate and return the next Normal distribution value.
def rand()
if @next
retval, @next = @next, nil
return retval
else
u = v = s = nil
loop do
u = Random.rand(-1.0..1.0)
v = Random.rand(-1.0..1.0)
s = u**2 + v**2
break if (s > 0.0) && (s <= 1.0)
end
f = Math.sqrt(-2.0 * Math.log(s) / s)
@next = v * f
return u * f
end
end
end
The Task
require('enumerable/statistics')
def show_stats_and_histogram(data, bins)
puts("size = #{data.length} mean = #{data.mean()} stddev = #{data.stdev()}")
hist = data.histogram(bins)
scale = 100.0 / hist.weights.max
inx_beg = nil
inx_end = nil
hist.weights.length.times do |inx|
histstars = (0.5 + (scale * hist.weights[inx])).to_i
inx_beg = inx if !inx_beg && (histstars > 0)
inx_end = inx if (histstars > 0)
end
(inx_beg..inx_end).each do |inx|
bincenter = 0.5 * (hist.edges[inx] + hist.edges[inx + 1])
histstars = (0.5 + (scale * hist.weights[inx])).to_i
puts('%6.2f: %s' % [bincenter, '*' * histstars])
end
end
puts
puts('Uniform random number generator:')
show_stats_and_histogram(1000000.times.map { Random.rand(-1.0..1.0) }, 20)
puts
puts('Normal random numbers using the Marsaglia polar method:')
gen_normal = NormalFromUniform.new
show_stats_and_histogram(100.times.map { gen_normal.rand }, 40)
show_stats_and_histogram(10000.times.map { gen_normal.rand }, 60)
show_stats_and_histogram(1000000.times.map { gen_normal.rand }, 120)
- Output:
Uniform random number generator: size = 1000000 mean = 0.0001483724103528628 stddev = 0.5773085740222473 -0.95: **************************************************************************************************** -0.85: **************************************************************************************************** -0.75: *************************************************************************************************** -0.65: *************************************************************************************************** -0.55: *************************************************************************************************** -0.45: **************************************************************************************************** -0.35: **************************************************************************************************** -0.25: **************************************************************************************************** -0.15: **************************************************************************************************** -0.05: *************************************************************************************************** 0.05: **************************************************************************************************** 0.15: **************************************************************************************************** 0.25: **************************************************************************************************** 0.35: *************************************************************************************************** 0.45: *************************************************************************************************** 0.55: **************************************************************************************************** 0.65: **************************************************************************************************** 0.75: **************************************************************************************************** 0.85: **************************************************************************************************** 0.95: *************************************************************************************************** Normal random numbers using the Marsaglia polar method: size = 100 mean = 0.1611961166227389 stddev = 0.9827581078952005 -2.30: ********** -2.10: -1.90: ********** -1.70: ******************** -1.50: -1.30: ********** -1.10: ********************************************************************** -0.90: ************************************************************ -0.70: ********************************************************************** -0.50: ************************************************************ -0.30: ******************************************************************************** -0.10: ******************** 0.10: ******************************************************************************** 0.30: **************************************************************************************************** 0.50: ****************************************************************************************** 0.70: ******************************************************************************** 0.90: ************************************************************ 1.10: ****************************** 1.30: ************************************************** 1.50: 1.70: ******************** 1.90: ************************************************** 2.10: ********** 2.30: ******************** size = 10000 mean = -0.004863294071004369 stddev = 0.9984395022628252 -3.30: * -3.10: * -2.90: ** -2.70: ** -2.50: **** -2.30: ********* -2.10: *********** -1.90: **************** -1.70: *********************** -1.50: ***************************** -1.30: ***************************************** -1.10: ************************************************* -0.90: ****************************************************************** -0.70: ****************************************************************************** -0.50: *************************************************************************************** -0.30: ********************************************************************************************* -0.10: *********************************************************************************************** 0.10: **************************************************************************************************** 0.30: ************************************************************************************** 0.50: ************************************************************************************ 0.70: ******************************************************************************* 0.90: **************************************************************** 1.10: *************************************************** 1.30: ******************************************** 1.50: ******************************* 1.70: ********************** 1.90: **************** 2.10: ********** 2.30: ****** 2.50: ***** 2.70: ** 2.90: * 3.10: * size = 1000000 mean = 0.0007049206911295587 stddev = 1.0000356747411552 -3.15: * -3.05: * -2.95: * -2.85: ** -2.75: ** -2.65: *** -2.55: **** -2.45: ***** -2.35: ****** -2.25: ******** -2.15: ********** -2.05: ************ -1.95: *************** -1.85: ****************** -1.75: ********************* -1.65: ************************* -1.55: ****************************** -1.45: *********************************** -1.35: **************************************** -1.25: ********************************************* -1.15: *************************************************** -1.05: ********************************************************* -0.95: *************************************************************** -0.85: ********************************************************************* -0.75: ************************************************************************** -0.65: ********************************************************************************* -0.55: ************************************************************************************* -0.45: ***************************************************************************************** -0.35: ******************************************************************************************** -0.25: ************************************************************************************************ -0.15: ************************************************************************************************** -0.05: **************************************************************************************************** 0.05: *************************************************************************************************** 0.15: ************************************************************************************************** 0.25: ************************************************************************************************ 0.35: ********************************************************************************************* 0.45: ****************************************************************************************** 0.55: ************************************************************************************* 0.65: ******************************************************************************** 0.75: ************************************************************************** 0.85: ********************************************************************* 0.95: *************************************************************** 1.05: ********************************************************* 1.15: **************************************************** 1.25: ********************************************** 1.35: **************************************** 1.45: ********************************** 1.55: ****************************** 1.65: ************************* 1.75: ********************** 1.85: ****************** 1.95: *************** 2.05: ************ 2.15: ********** 2.25: ******** 2.35: ****** 2.45: ***** 2.55: **** 2.65: *** 2.75: ** 2.85: ** 2.95: * 3.05: * 3.15: *
Run BASIC
s = 100000
h$ = "============================================================="
h$ = h$ + h$
dim ndis(s)
' mean and standard deviation.
mx = -9999
mn = 9999
sum = 0
sumSqr = 0
for i = 1 to s ' find minimum and maximum
ms = rnd(1)
ss = rnd(1)
nd = (-2 * log(ms))^0.5 * cos(2 *3.14159265 * ss) ' normal distribution
ndis(i) = nd
mx = max(mx, nd)
mn = min(mn, nd)
sum = sum + nd
sumSqr = sumSqr + nd ^ 2
next i
mean = sum / s
range = mx - mn
print "Samples :"; s
print "Largest :"; mx
print "Smallest :"; mn
print "Range :"; range
print "Mean :"; mean
print "Stand Dev :"; (sumSqr /s -mean^2)^0.5
'Show chart of histogram
nBins = 50
dim bins(nBins)
for i = 1 to s
z = int((ndis(i) -mn) /range *nBins)
bins(z) = bins(z) + 1
mb = max(bins(z),mb)
next i
for b = 0 to nBins -1
print using("##",b);" ";using("#####",bins(b));" ";left$(h$,(bins(b) / mb) * 90)
next b
END
- Output:
Samples :100000 Largest :4.61187177 Smallest :-4.21695424 Range :8.82882601 Mean :-9.25042513e-4 Stand Dev :1.00680067 = == === ===== ======== ============= ================= ======================= ============================== ======================================= =============================================== ========================================================= =================================================================== =========================================================================== =================================================================================== ======================================================================================= ========================================================================================== ======================================================================================== ====================================================================================== ================================================================================= ============================================================================ ================================================================== ======================================================== ============================================== ===================================== ============================ ===================== =============== ========== ======= ===== === = =
Rust
//! Rust rosetta example for normal distribution
use math::{histogram::Histogram, traits::ToIterator};
use rand;
use rand_distr::{Distribution, Normal};
/// Returns the mean of the provided samples
///
/// ## Arguments
/// * data -- reference to float32 array
fn mean(data: &[f32]) -> Option<f32> {
let sum: f32 = data.iter().sum();
Some(sum / data.len() as f32)
}
/// Returns standard deviation of the provided samples
///
/// ## Arguments
/// * data -- reference to float32 array
fn standard_deviation(data: &[f32]) -> Option<f32> {
let mean = mean(data).expect("invalid mean");
let sum = data.iter().fold(0.0, |acc, &x| acc + (x - mean).powi(2));
Some((sum / data.len() as f32).sqrt())
}
/// Prints a histogram in the shell
///
/// ## Arguments
/// * data -- reference to float32 array
/// * maxwidth -- the maxwidth of the histogram in # of characters
/// * bincount -- number of bins in the histogram
/// * ch -- character used to plot the graph
fn print_histogram(data: &[f32], maxwidth: usize, bincount: usize, ch: char) {
let min_val = data.iter().cloned().fold(f32::NAN, f32::min);
let max_val = data.iter().cloned().fold(f32::NAN, f32::max);
let histogram = Histogram::new(Some(&data.to_vec()), bincount, min_val, max_val).unwrap();
let max_bin_value = histogram.get_counters().iter().max().unwrap();
println!();
for x in histogram.to_iter() {
let (bin_min, bin_max, freq) = x;
let bar_width = (((freq as f64) / (*max_bin_value as f64)) * (maxwidth as f64)) as u32;
let bar_as_string = (1..bar_width).fold(String::new(), |b, _| b + &ch.to_string());
println!(
"({:>6},{:>6}) |{} {:.2}%",
format!("{:.2}", bin_min),
format!("{:.2}", bin_max),
bar_as_string,
(freq as f64) * 100.0 / (data.len() as f64)
);
}
println!();
}
/// Runs the demo to generate normal distribution of three different sample sizes
fn main() {
let expected_mean: f32 = 0.0;
let expected_std_deviation: f32 = 4.0;
let normal = Normal::new(expected_mean, expected_std_deviation).unwrap();
let mut rng = rand::thread_rng();
for &number_of_samples in &[1000, 10_000, 1_000_000] {
let data: Vec<f32> = normal
.sample_iter(&mut rng)
.take(number_of_samples)
.collect();
println!("Statistics for sample size {}:", number_of_samples);
println!("\tMean: {:?}", mean(&data).expect("invalid mean"));
println!(
"\tStandard deviation: {:?}",
standard_deviation(&data).expect("invalid standard deviation")
);
print_histogram(&data, 80, 40, '-');
}
}
- Output:
Statistics for sample size 1000: Mean: -0.11356559 Standard deviation: 4.0244555 (-13.81,-13.12) | 0.10% (-13.12,-12.44) | 0.00% (-12.44,-11.75) | 0.10% (-11.75,-11.06) | 0.20% (-11.06,-10.38) |- 0.30% (-10.38, -9.69) | 0.10% ( -9.69, -9.01) |--- 0.50% ( -9.01, -8.32) |---- 0.60% ( -8.32, -7.64) |------ 0.80% ( -7.64, -6.95) |-------------- 1.60% ( -6.95, -6.27) |----------------- 1.90% ( -6.27, -5.58) |------------------------ 2.60% ( -5.58, -4.90) |----------------------- 2.50% ( -4.90, -4.21) |---------------------------------------- 4.20% ( -4.21, -3.53) |------------------------------------- 3.90% ( -3.53, -2.84) |------------------------------------------------- 5.10% ( -2.84, -2.15) |---------------------------------------------- 4.80% ( -2.15, -1.47) |------------------------------------------------------------------------ 7.40% ( -1.47, -0.78) |---------------------------------------------------------- 6.00% ( -0.78, -0.10) |----------------------------------------------------------------------- 7.30% ( -0.10, 0.59) |------------------------------------------------------------------------------- 8.10% ( 0.59, 1.27) |----------------------------------------------------------------------- 7.30% ( 1.27, 1.96) |------------------------------------------------- 5.10% ( 1.96, 2.64) |------------------------------------------------------------ 6.20% ( 2.64, 3.33) |----------------------------------------- 4.30% ( 3.33, 4.01) |----------------------------- 3.10% ( 4.01, 4.70) |------------------------------------- 3.90% ( 4.70, 5.39) |-------------------------- 2.80% ( 5.39, 6.07) |---------------------- 2.40% ( 6.07, 6.76) |---------------- 1.80% ( 6.76, 7.44) |---------------- 1.80% ( 7.44, 8.13) |--------- 1.10% ( 8.13, 8.81) |---------- 1.20% ( 8.81, 9.50) | 0.20% ( 9.50, 10.18) | 0.00% ( 10.18, 10.87) | 0.10% ( 10.87, 11.55) |- 0.30% ( 11.55, 12.24) | 0.10% ( 12.24, 12.92) | 0.10% ( 12.92, 13.61) | 0.10% Statistics for sample size 10000: Mean: 0.02012564 Standard deviation: 3.9697735 (-15.80,-14.99) | 0.02% (-14.99,-14.18) | 0.04% (-14.18,-13.37) | 0.04% (-13.37,-12.56) | 0.04% (-12.56,-11.75) | 0.09% (-11.75,-10.94) | 0.08% (-10.94,-10.13) |- 0.25% (-10.13, -9.32) |--- 0.42% ( -9.32, -8.51) |----- 0.67% ( -8.51, -7.70) |--------- 1.10% ( -7.70, -6.89) |------------- 1.48% ( -6.89, -6.08) |------------------ 1.98% ( -6.08, -5.27) |-------------------------- 2.82% ( -5.27, -4.45) |------------------------------------ 3.80% ( -4.45, -3.64) |--------------------------------------------- 4.66% ( -3.64, -2.83) |------------------------------------------------------- 5.72% ( -2.83, -2.02) |------------------------------------------------------------------ 6.85% ( -2.02, -1.21) |---------------------------------------------------------------------------- 7.80% ( -1.21, -0.40) |---------------------------------------------------------------------------- 7.79% ( -0.40, 0.41) |------------------------------------------------------------------------------- 8.09% ( 0.41, 1.22) |----------------------------------------------------------------------------- 7.89% ( 1.22, 2.03) |--------------------------------------------------------------------------- 7.76% ( 2.03, 2.84) |-------------------------------------------------------------------- 7.06% ( 2.84, 3.65) |------------------------------------------------------- 5.74% ( 3.65, 4.46) |-------------------------------------------- 4.64% ( 4.46, 5.27) |-------------------------------------- 4.00% ( 5.27, 6.08) |---------------------------- 3.03% ( 6.08, 6.89) |------------------- 2.07% ( 6.89, 7.71) |-------------- 1.54% ( 7.71, 8.52) |-------- 0.97% ( 8.52, 9.33) |----- 0.61% ( 9.33, 10.14) |-- 0.36% ( 10.14, 10.95) |- 0.25% ( 10.95, 11.76) | 0.19% ( 11.76, 12.57) | 0.08% ( 12.57, 13.38) | 0.02% ( 13.38, 14.19) | 0.01% ( 14.19, 15.00) | 0.03% ( 15.00, 15.81) | 0.00% ( 15.81, 16.62) | 0.01% Statistics for sample size 1000000: Mean: -0.004743685 Standard deviation: 4.0006065 (-20.00,-19.02) | 0.00% (-19.02,-18.04) | 0.00% (-18.04,-17.06) | 0.00% (-17.06,-16.07) | 0.00% (-16.07,-15.09) | 0.00% (-15.09,-14.11) | 0.01% (-14.11,-13.13) | 0.03% (-13.13,-12.15) | 0.06% (-12.15,-11.16) | 0.14% (-11.16,-10.18) |- 0.28% (-10.18, -9.20) |--- 0.53% ( -9.20, -8.22) |------ 0.95% ( -8.22, -7.24) |----------- 1.51% ( -7.24, -6.25) |------------------ 2.40% ( -6.25, -5.27) |--------------------------- 3.48% ( -5.27, -4.29) |-------------------------------------- 4.82% ( -4.29, -3.31) |-------------------------------------------------- 6.27% ( -3.31, -2.32) |------------------------------------------------------------- 7.62% ( -2.32, -1.34) |----------------------------------------------------------------------- 8.77% ( -1.34, -0.36) |----------------------------------------------------------------------------- 9.58% ( -0.36, 0.62) |------------------------------------------------------------------------------- 9.74% ( 0.62, 1.60) |---------------------------------------------------------------------------- 9.39% ( 1.60, 2.59) |-------------------------------------------------------------------- 8.49% ( 2.59, 3.57) |---------------------------------------------------------- 7.30% ( 3.57, 4.55) |----------------------------------------------- 5.86% ( 4.55, 5.53) |----------------------------------- 4.45% ( 5.53, 6.51) |------------------------ 3.16% ( 6.51, 7.50) |---------------- 2.09% ( 7.50, 8.48) |--------- 1.34% ( 8.48, 9.46) |----- 0.81% ( 9.46, 10.44) |-- 0.46% ( 10.44, 11.42) | 0.23% ( 11.42, 12.41) | 0.11% ( 12.41, 13.39) | 0.06% ( 13.39, 14.37) | 0.02% ( 14.37, 15.35) | 0.01% ( 15.35, 16.34) | 0.00% ( 16.34, 17.32) | 0.00% ( 17.32, 18.30) | 0.00% ( 18.30, 19.28) | 0.00%
SAS
data test;
n=100000;
twopi=2*constant('pi');
do i=1 to n;
u=ranuni(0);
v=ranuni(0);
r=sqrt(-2*log(u));
x=r*cos(twopi*v);
y=r*sin(twopi*v);
z=rannor(0);
output;
end;
keep x y z;
proc means mean stddev;
proc univariate;
histogram /normal;
run;
/*
Variable Mean Std Dev
----------------------------------------
x -0.0052720 0.9988467
y 0.000023995 1.0019996
z 0.0012857 1.0056536
*/
Sidef
define τ = Num.tau
func normdist (m, σ) {
var r = sqrt(-2 * 1.rand.log)
var Θ = (τ * 1.rand)
r * Θ.cos * σ + m
}
var size = 100_000
var mean = 50
var stddev = 4
var dataset = size.of { normdist(mean, stddev) }
var m = (dataset.sum / size)
say ("m: #{m}")
var σ = sqrt(dataset »**» 2 -> sum / size - m**2)
say ("s: #{σ}")
var hash = Hash()
dataset.each { |n| hash{ n.round } := 0 ++ }
var scale = (180 * stddev / size)
const subbar = < ⎸ ▏ ▎ ▍ ▌ ▋ ▊ ▉ █ >
for i in (hash.keys.map{.to_i}.sort) {
var x = (hash{i} * scale)
var full = x.int
var part = (8 * (x - full))
say (i, "\t", '█' * full, subbar[part])
}
- Output:
m: 49.99538275618550306540055142077589 s: 4.00295544816687358837821680496471 33 ⎸ 34 ⎸ 35 ⎸ 36 ▏ 37 ▎ 38 ▊ 39 █▋ 40 ███▏ 41 ██████▏ 42 █████████▍ 43 ███████████████▌ 44 ███████████████████████▋ 45 ████████████████████████████████▍ 46 ████████████████████████████████████████████▎ 47 █████████████████████████████████████████████████████▍ 48 ███████████████████████████████████████████████████████████████▍ 49 █████████████████████████████████████████████████████████████████████▌ 50 ████████████████████████████████████████████████████████████████████████▋ 51 █████████████████████████████████████████████████████████████████████▊ 52 ██████████████████████████████████████████████████████████████▏ 53 ████████████████████████████████████████████████████▉ 54 ███████████████████████████████████████████▉ 55 █████████████████████████████████▎ 56 ███████████████████████⎸ 57 ███████████████▋ 58 █████████▋ 59 █████▍ 60 ███▍ 61 █▊ 62 ▋ 63 ▍ 64 ▏ 65 ⎸ 66 ⎸
Stata
Pairs of normal numbers are generated from pairs of uniform numbers using the Box-Muller method. A normal density is added to the histogram for comparison. See histogram in Stata help. A Q-Q plot is also drawn.
clear all
set obs 100000
gen u=runiform()
gen v=runiform()
gen r=sqrt(-2*log(u))
gen x=r*cos(2*_pi*v)
gen y=r*sin(2*_pi*v)
gen z=rnormal()
sum x y z
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------------
x | 100,000 .0025861 1.002346 -4.508192 4.164336
y | 100,000 .0017389 1.001586 -4.631144 4.460274
z | 100,000 .005054 .9998861 -5.134265 4.449522
hist x, normal
hist y, normal
hist z, normal
qqplot x z, msize(tiny)
Tcl
package require Tcl 8.5
# Uses the Box-Muller transform to compute a pair of normal random numbers
proc tcl::mathfunc::nrand {mean stddev} {
variable savednormalrandom
if {[info exists savednormalrandom]} {
return [expr {$savednormalrandom*$stddev + $mean}][unset savednormalrandom]
}
set r [expr {sqrt(-2*log(rand()))}]
set theta [expr {2*3.1415927*rand()}]
set savednormalrandom [expr {$r*sin($theta)}]
expr {$r*cos($theta)*$stddev + $mean}
}
proc stats {size {slotfactor 10}} {
set sum 0.0
set sum2 0.0
for {set i 0} {$i < $size} {incr i} {
set r [expr { nrand(0.5, 0.2) }]
incr histo([expr {int(floor($r*$slotfactor))}])
set sum [expr {$sum + $r}]
set sum2 [expr {$sum2 + $r**2}]
}
set mean [expr {$sum / $size}]
set stddev [expr {sqrt($sum2/$size - $mean**2)}]
puts "$size numbers"
puts "Mean: $mean"
puts "StdDev: $stddev"
foreach i [lsort -integer [array names histo]] {
puts [string repeat "*" [expr {$histo($i)*350/int($size)}]]
}
}
stats 100
puts ""
stats 1000
puts ""
stats 10000
puts ""
stats 100000 20
Sample output:
100 numbers Mean: 0.49355955990390254 StdDev: 0.19651396178121985 *** ******* ************** *********************************** ******************************************************** ****************************************************************** ************************************************************************* ****************************************** ************************************** ************** 1000 numbers Mean: 0.5066940614105869 StdDev: 0.2016794788065389 * ***** ************** **************************** ********************************************************** **************************************************************** ************************************************************* ****************************************************** *********************************** ************ ********* * 10000 numbers Mean: 0.49980964730768285 StdDev: 0.1968441612522318 * ***** *************** ******************************* ***************************************************** ****************************************************************** ******************************************************************* **************************************************** ********************************* *************** ***** * 100000 numbers Mean: 0.49960438950922254 StdDev: 0.20060211160998606 * ** *** ****** ********* ************** ****************** *********************** ***************************** ******************************** ********************************** ********************************** ******************************** **************************** *********************** ****************** ************* ********* ****** *** ** *
The blank lines in the output are where the number of samples is too small to even merit a single unit on the histogram.
VBA
Public Sub standard_normal()
Dim s() As Variant, bins(71) As Single
ReDim s(20000)
For i = 1 To 20000
s(i) = WorksheetFunction.Norm_S_Inv(Rnd())
Next i
For i = -35 To 35
bins(i + 36) = i / 10
Next i
Debug.Print "sample size"; UBound(s), "mean"; mean(s), "standard deviation"; standard_deviation(s)
t = WorksheetFunction.Frequency(s, bins)
For i = -35 To 35
Debug.Print Format((i - 1) / 10, "0.00");
Debug.Print "-"; Format(i / 10, "0.00"),
Debug.Print String$(t(i + 36, 1) / 10, "X");
Debug.Print
Next i
End Sub
- Output:
sample size 20000 mean-5,26306310478751E-03 standard deviation 1,00355037427319 -3,60--3,50 -3,50--3,40 -3,40--3,30 -3,30--3,20 -3,20--3,10 -3,10--3,00 -3,00--2,90 XX -2,90--2,80 X -2,80--2,70 XX -2,70--2,60 XX -2,60--2,50 XXX -2,50--2,40 XXXX -2,40--2,30 XXXXX -2,30--2,20 XXXXXXXX -2,20--2,10 XXXXXXXX -2,10--2,00 XXXXXXXXXXX -2,00--1,90 XXXXXXXXXXXXX -1,90--1,80 XXXXXXXXXXXXXXX -1,80--1,70 XXXXXXXXXXXXXXXX -1,70--1,60 XXXXXXXXXXXXXXXXXXXX -1,60--1,50 XXXXXXXXXXXXXXXXXXXXXXXX -1,50--1,40 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -1,40--1,30 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -1,30--1,20 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -1,20--1,10 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -1,10--1,00 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -1,00--0,90 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,90--0,80 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,80--0,70 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,70--0,60 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,60--0,50 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,50--0,40 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,40--0,30 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,30--0,20 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,20--0,10 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX -0,10-0,00 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,00-0,10 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,10-0,20 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,20-0,30 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,30-0,40 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,40-0,50 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,50-0,60 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,60-0,70 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,70-0,80 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,80-0,90 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0,90-1,00 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 1,00-1,10 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 1,10-1,20 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 1,20-1,30 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 1,30-1,40 XXXXXXXXXXXXXXXXXXXXXXXXXXXXX 1,40-1,50 XXXXXXXXXXXXXXXXXXXXXXXXXX 1,50-1,60 XXXXXXXXXXXXXXXXXXXXXXXXX 1,60-1,70 XXXXXXXXXXXXXXXXXXXXXX 1,70-1,80 XXXXXXXXXXXXXXXXXX 1,80-1,90 XXXXXXXXXXXXXXX 1,90-2,00 XXXXXXXXXXX 2,00-2,10 XXXXXXXXXXXX 2,10-2,20 XXXXXXX 2,20-2,30 XXXXXX 2,30-2,40 XXXXX 2,40-2,50 XXX 2,50-2,60 XXXX 2,60-2,70 XX 2,70-2,80 XX 2,80-2,90 X 2,90-3,00 X 3,00-3,10 X 3,10-3,20 X 3,20-3,30 3,30-3,40 3,40-3,50
Wren
import "random" for Random
import "./fmt" for Fmt
import "./math" for Nums
var rgen = Random.new()
// Box-Muller method from Wikipedia
var normal = Fn.new { |mu, sigma|
var u1 = rgen.float()
var u2 = rgen.float()
var mag = sigma * (-2 * u1.log).sqrt
var z0 = mag * (2 * Num.pi * u2).cos + mu
var z1 = mag * (2 * Num.pi * u2).sin + mu
return [z0, z1]
}
var N = 100000
var NUM_BINS = 12
var HIST_CHAR = "■"
var HIST_CHAR_SIZE = 250
var bins = List.filled(NUM_BINS, 0)
var binSize = 0.1
var samples = List.filled(N, 0)
var mu = 0.5
var sigma = 0.25
for (i in 0...N/2) {
var rns = normal.call(mu, sigma)
for (j in 0..1) {
var rn = rns[j]
var bn
if (rn < 0) {
bn = 0
} else if (rn >= 1) {
bn = 11
} else {
bn = (rn/binSize).floor + 1
}
bins[bn] = bins[bn] + 1
samples[i*2 + j] = rn
}
}
Fmt.print("Normal distribution with mean $0.2f and S/D $0.2f for $,d samples:\n", mu, sigma, N)
System.print(" Range Number of samples within that range")
for (i in 0...NUM_BINS) {
var hist = HIST_CHAR * (bins[i] / HIST_CHAR_SIZE).round
if (i == 0) {
Fmt.print(" -∞ ..< 0.00 $s $,d", hist, bins[0])
} else if (i < NUM_BINS - 1) {
Fmt.print("$4.2f ..< $4.2f $s $,d", binSize * (i-1), binSize * i, hist, bins[i])
} else {
Fmt.print("1.00 ... +∞ $s $,d", hist, bins[NUM_BINS - 1])
}
}
Fmt.print("\nActual mean for these samples : $0.5f", Nums.mean(samples))
Fmt.print("Actual S/D for these samples : $0.5f", Nums.stdDev(samples))
- Output:
Specimen run:
Normal distribution with mean 0.50 and S/D 0.25 for 100,000 samples: Range Number of samples within that range -∞ ..< 0.00 ■■■■■■■■■ 2,243 0.00 ..< 0.10 ■■■■■■■■■■■■■ 3,250 0.10 ..< 0.20 ■■■■■■■■■■■■■■■■■■■■■■■■ 5,977 0.20 ..< 0.30 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 9,723 0.30 ..< 0.40 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 13,104 0.40 ..< 0.50 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 15,601 0.50 ..< 0.60 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 15,469 0.60 ..< 0.70 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 13,334 0.70 ..< 0.80 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 9,659 0.80 ..< 0.90 ■■■■■■■■■■■■■■■■■■■■■■■■ 6,119 0.90 ..< 1.00 ■■■■■■■■■■■■■ 3,173 1.00 ... +∞ ■■■■■■■■■ 2,348 Actual mean for these samples : 0.50099 Actual S/D for these samples : 0.25051
zkl
fcn norm2{ // Box-Muller
const PI2=(0.0).pi*2;;
rnd:=(0.0).random.fp(1); // random number in [0,1), using partial application
r,a:=(-2.0*rnd().log()).sqrt(), PI2*rnd();
return(r*a.cos(), r*a.sin()); // z0,z1
}
const N=100000, BINS=12, SIG=3, SCALE=500;
var sum=0.0,sumSq=0.0, h=BINS.pump(List(),0); // (0,0,0,...)
fcn accum(v){
sum+=v;
sumSq+=v*v;
b:=(v + SIG)*BINS/SIG/2;
if(0<=b<BINS) h[b]+=1;
};
Partial application: rnd() --> (0.0).random(1). Basically, the fp method fixes the call parameters, which are then used when the partial thing is run.
foreach i in (N/2){ v1,v2:=norm2(); accum(v1); accum(v2); }
println("Samples: %,d".fmt(N));
println("Mean: ", m:=sum/N);
println("Stddev: ", (sumSq/N - m*m).sqrt());
foreach p in (h){ println("*"*(p/SCALE)) }
- Output:
Samples: 100,000 Mean: 0.0005999 Stddev: 1.003 * *** ******** ****************** ***************************** ************************************** ************************************** ***************************** ****************** ******** *** *
- Programming Tasks
- Solutions by Programming Task
- C
- C sharp
- Math.Net
- C++
- D
- EDSAC order code
- Elixir
- Factor
- Fortran
- FreeBASIC
- Go
- Haskell
- J
- Java
- Jq
- Julia
- Kotlin
- Lasso
- Liberty BASIC
- Lua
- Maple
- Mathematica
- Wolfram Language
- MATLAB
- Octave
- Nim
- PARI/GP
- Pascal
- Perl
- Phix
- PureBasic
- Python
- R
- Racket
- Raku
- REXX
- Ruby
- Enumerable-statistics
- Run BASIC
- Rust
- Math
- Rand
- Rand distr
- SAS
- Sidef
- Stata
- Tcl
- VBA
- Wren
- Wren-fmt
- Wren-math
- Zkl