Verify distribution uniformity/Naive: Difference between revisions

Content added Content deleted
No edit summary
(→‎{{header|REXX}}: changed comments and whitespace, changed some variable names to be more descriptive, added commas to output for easier reading the results.)
Line 1,459: Line 1,459:
=={{header|REXX}}==
=={{header|REXX}}==
<lang rexx>/*REXX program simulates a number of trials of a random digit and show it's skew %. */
<lang rexx>/*REXX program simulates a number of trials of a random digit and show it's skew %. */
parse arg f t d s . /*obtain arguments (options) from C.L. */
parse arg func times delta seed . /*obtain arguments (options) from C.L. */
if f=='' | f=="," then f='RANDOM' /*function not specified? Use default.*/
if func=='' | func=="," then func= 'RANDOM' /*function not specified? Use default.*/
if t=='' | t=="," then t=1000000 /*times " " " " */
if times=='' | times=="," then times= 1000000 /*times " " " " */
if d=='' | d=="," then d=1/2 /*delta% " " " " */
if delta=='' | delta=="," then delta= 1/2 /*delta% " " " " */
if s\=='' then call random ,,s /*use some RAND seed for repeatability.*/
if datatype(seed, 'W') then call random ,,seed /*use some RAND seed for repeatability.*/
highDig=9 /*use this var for the highest digit. */
highDig=9 /*use this var for the highest digit. */
!.=0 /*initialize all possible random trials*/
!.=0 /*initialize all possible random trials*/
do t /* [↓] perform a bunch of trials. */
do times /* [↓] perform a bunch of trials. */
if f=='RANDOM' then ?=random(0,highDig) /*use the RANDOM BIF function.*/
if func=='RANDOM' then ?=random(highDig) /*use RANDOM function.*/
else interpret '?='f"(0,"highDig')' /*use the (specified) " */
else interpret '?='func"(0,"highDig')' /* " specified " */
!.?=!.?+1 /*bump the invocation counter.*/
!.?=!.? + 1 /*bump the invocation counter.*/
end /*t*/ /* [↑] store trials ───► pigeonholes. */
end /*t*/ /* [↑] store trials ───► pigeonholes. */
/* [↓] compute the digit's skewness. */
/* [↓] compute the digit's skewness. */
g=t / (1+highDig) /*calculate number of each digit throw.*/
g=times / (1 + highDig) /*calculate number of each digit throw.*/
OK?= 'OK skewed' /*words to show "skewed" or if "OK".*/
w=max(9, length( commas(times) ) ) /*maximum length of number of trials.*/
w=max(8, length(t)) /*maximum length of number of trials.*/
pad=left('', 9) /*this is used for output indentation. */
pad=left('', 9) /*this is used for output indentation. */
say pad 'digit' center("hits",w) ' skew ' "skew%" 'result' /*header. */
say pad 'digit' center(" hits", w) ' skew ' "skew %" 'result' /*header. */
say pad '─────' center('',w,'─') '──────' "─────" '──────' /*separator.*/
say pad '─────' center('', w, '─') '──────' "──────" '──────' /*separator*/
/** [↑] show header and the separator.*/
/** [↑] show header and the separator.*/
do k=0 to highDig /*process each of the possible digits. */
do k=0 to highDig /*process each of the possible digits. */
skew=g-!.k /*calculate the skew for the digit. */
skew=g - !.k /*calculate the skew for the digit. */
skewPC=(1- (g-abs(skew)) / g) * 100 /* " " " percentage for dig*/
skewPC=(1- (g-abs(skew)) / g) * 100 /* " " " percentage for dig*/
ok=center(word(ok?, 1 + (skewPC>d)), 6) /*it's gotta be one of skewed or xx%*/
say pad center(k, 5) right( commas(!.k), w) right(skew, 6) ,
say pad center(k,5) right(!.k,w) right(skew,6) format(skewPC,,3) ok
right( format(skewPC, , 3), 6) center( word('ok skewed', 1+(skewPC>delta)), 6)
end /*k*/
end /*k*/


say pad '─────' center('',w,'─') '──────' "─────" '──────' /*separator. */
say pad '─────' center('', w, '─') '──────' "─────" '──────' /*separator. */
y=5+1+w+1+6+1+6+1+6 /*the width. */
y=5+1+w+1+6+1+6+1+6 /*the width. */
say pad center(" (with " t ' trials)',y) /*# trials. */
say pad center(" (with " commas(times) ' trials)' , y) /*# trials. */
say pad center(" (skewed when exceeds " d'%)',y) /*skewed note.*/
say pad center(" (skewed when exceeds " delta'%)' , y) /*skewed note.*/
/*stick a fork in it, we're all done. */</lang>
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
Execution note: &nbsp; quite a few runs were needed and the skew'''%''' lowered before a skewed result was obtained.
commas: procedure; parse arg _; n=_'.9'; #=123456789
e=verify(n, #'0', , verify(n, #"0.",'M') ) - 4
do j=e to verify(n, #, "M") by -3; _=insert(',', _, j); end; return _</lang>
Execution note: &nbsp; quite a few runs were needed before a skewed result was obtained.
<br><br>
<br><br>
'''output''' when using the default inputs:
'''output''' when using the default inputs:
<pre>
<pre>
digit hits skew skew% result
digit hits skew skew % result
───── ──────── ────── ───── ──────
───── ───────── ────── ────── ──────
0 99757 243 0.243 OK
0 99,739 261 0.261 ok
1 100226 -226 0.226 OK
1 99,819 181 0.181 ok
2 100605 -605 0.605 skewed
2 100,463 -463 0.463 ok
3 100005 -5 0.005 OK
3 99,787 213 0.213 ok
4 99670 330 0.330 OK
4 99,632 368 0.368 ok
5 100011 -11 0.011 OK
5 100,787 -787 0.787 skewed
6 100100 -100 0.100 OK
6 99,704 296 0.296 ok
7 99513 487 0.487 OK
7 99,605 395 0.395 ok
8 99884 116 0.116 OK
8 100,488 -488 0.488 ok
9 100229 -229 0.229 OK
9 99,976 24 0.024 ok
───── ──────── ────── ───── ──────
───── ───────── ────── ───── ──────
(with 1000000 trials)
(with 1,000,000 trials)
(skewed when exceeds 0.5%)
(skewed when exceeds 0.5%)
</pre>
</pre>