forked from ocaml/ocaml
-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Simple chi-square test for functions from module Random
This test appeared during review of PR ocaml#9489, but seems useful independently of this PR.
- Loading branch information
1 parent
be8ead9
commit dae1460
Showing
1 changed file
with
68 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,68 @@ | ||
(* TEST | ||
*) | ||
|
||
(* A basic chi-square test to detect simple errors in the Random module. *) | ||
|
||
(* Accumulate [n] samples from function [f] and check the chi-square. | ||
Only the low 8 bits of the result of [f] are sampled. *) | ||
|
||
let chisquare n f = | ||
let r = 256 in | ||
let freq = Array.make r 0 in | ||
for i = 0 to n - 1 do | ||
let t = f () land 0xFF in | ||
freq.(t) <- freq.(t) + 1 | ||
done; | ||
let expected = float n /. float r in | ||
let t = | ||
Array.fold_left | ||
(fun s x -> let d = float x -. expected in d *. d +. s) | ||
0.0 freq in | ||
let chi2 = t /. expected in | ||
let degfree = float r -. 1.0 in | ||
(* The degree of freedom is high, so we approximate as a normal | ||
distribution with mean equal to degfree and variance 2 * degfree. | ||
Four sigmas correspond to a 99.9968% confidence interval. | ||
(Without the approximation, the confidence interval seems to be 99.986%.) | ||
*) | ||
chi2 <= degfree +. 4.0 *. sqrt (2.0 *. degfree) | ||
|
||
let test name f = | ||
if not (chisquare 100_000 f) | ||
then Printf.printf "%s: suspicious result\n%!" name | ||
|
||
let _ = | ||
test "Random.bits (bits 0-7)" | ||
Random.bits; | ||
test "Random.bits (bits 12-19)" | ||
(fun () -> Random.bits() lsr 12); | ||
test "Random.bits (bits 22-29)" | ||
(fun () -> Random.bits() lsr 22); | ||
test "Random.int 2^26 (bits 0-7)" | ||
(fun () -> Random.int (1 lsl 26)); | ||
test "Random.int 2^26 (bits 18-25)" | ||
(fun () -> Random.int (1 lsl 26) lsr 18); | ||
test "Random.int (256 * p) / p" | ||
(fun () -> Random.int (256 * 853187) / 853187); | ||
test "Random.float 1.0 (first 8 bits)" | ||
(fun () -> int_of_float (Random.float 1.0 *. 256.0)); | ||
test "Random.float 1.0 (next 8 bits)" | ||
(fun () -> int_of_float (Random.float 1.0 *. 65536.0)); | ||
test "Random.int32 2^30 (bits 0-7)" | ||
(fun () -> Int32.to_int (Random.int32 0x40000000l)); | ||
test "Random.int32 2^30 (bits 20-27)" | ||
(fun () -> Int32.(to_int (shift_right (Random.int32 0x40000000l) 20))); | ||
test "Random.int32 (256 * p) / p" | ||
(let p = 7048673l in | ||
fun () -> Int32.(to_int (div (Random.int32 (mul 256l p)) p))); | ||
test "Random.int64 2^60 (bits 0-7)" | ||
(fun () -> Int64.to_int (Random.int64 0x1000000000000000L)); | ||
test "Random.int64 2^60 (bits 30-37)" | ||
(fun () -> Int64.(to_int (shift_right (Random.int64 0x1000000000000000L) | ||
30))); | ||
test "Random.int64 2^60 (bits 52-59)" | ||
(fun () -> Int64.(to_int (shift_right (Random.int64 0x1000000000000000L) | ||
52))); | ||
test "Random.int64 (256 * p) / p" | ||
(let p = 16430454264262693L in | ||
fun () -> Int64.(to_int (div (Random.int64 (mul 256L p)) p))) |