1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
|
import Prando from 'prando';
import { none } from 'fp-ts/lib/Option.js';
import * as G from './gen.js';
import { test, expect } from 'vitest';
test('can generate and shrink bytes', async () => {
const prng = new Prando(Math.random() * 1000);
const gen = new G.BytesGen(prng);
const bytes: number[][] = [];
for (let i = 1; i < 10; i++) {
bytes.push(gen.generate(i));
}
console.log('Generated bytes:', bytes);
const replay = gen.replay();
for (let i = 1; i < 10; i++) {
const replayed = replay.generate(i);
expect(replayed, `Replaying at ${i}`).toEqual(bytes[i - 1]);
}
});
test('can generate an arbitrary integer', () => {
const prng = new Prando(Math.random() * 1000);
const gen = G.genInt();
const value = gen.generate();
expect(typeof value).toBe('number');
});
test('can generate N arbitrary integers', () => {
const prng = new Prando(Math.random() * 1000);
const gen = G.genInt();
const value = gen.generateN(5);
expect(value.length).toBe(5);
});
test('can generate arbitrary integer within bounds', () => {
const seed = Math.random() * 1000;
const gen = G.genIntWithin(1000, 200000, seed);
for (let i = 0; i < 100; i++) {
const value = gen.generate();
expect(value).toBeGreaterThanOrEqual(1000);
expect(value).toBeLessThanOrEqual(200000);
}
});
test('generated integers are uniformly distributed', () => {
const seed = Math.random() * 1000;
const lb = 0;
const ub = 65536;
const gen = G.genIntWithin(lb, ub, seed);
const range = ub - lb;
const buckets = 128;
const numSamples = 100000;
// indexed by lower bound
const counts: Record<number, number> = {};
for (let i = 1; i <= buckets; i++) {
counts[Math.ceil(range / buckets) * i + lb] = 0;
}
const values = [];
for (let i = 0; i < numSamples; i++) {
const value = gen.generate();
values.push(value);
const bucket = Object.keys(counts).find((b) => value < parseInt(b));
if (!bucket) {
console.log('ignoring bucket for value', value);
} else {
const b = parseInt(bucket);
counts[b] = (counts[b] || 0) + 1;
}
}
// compute χ² statistic
const expectedCount = numSamples / buckets;
const chiSquared = Object.values(counts).reduce((acc, count) => {
const diff = count - expectedCount;
return acc + (diff * diff) / expectedCount;
});
// compute p-value and Cramér's V.
// In our case the p-value will
// systematically be negligible and therefore will reject the null
// hypothesis, eg. that the frequency distribution of numbers we
// generate is truly uniform. However, Cramér's V is a better
// measure of effect size and will be small if the distribution is
// uniform.
const cramerv = computeCramersV(chiSquared, numSamples, buckets);
const pValue = computeChiSquaredPValue(chiSquared, buckets);
console.log(`Chi-squared: ${chiSquared}, p-value: ${pValue}, Cramer V: ${cramerv}`);
expect(cramerv).toBeLessThan(0.01);
});
/**
* Computes the right-tailed p-value for a Chi-Squared distribution
* using the Wilson-Hilferty approximation.
*
* @param chiSquareValue The calculated empirical chi-square statistic (e.g., 1190.372)
* @param numBuckets The total number of bins/categories in your test (e.g., 100)
* @returns The p-value between 0.0 and 1.0
*/
function computeChiSquaredPValue(chiSquareValue: number, numBuckets: number): number {
// Degrees of Freedom = Number of buckets - 1
const df = numBuckets - 1;
if (df <= 0) {
throw new Error("Number of buckets must be greater than 1.");
}
if (chiSquareValue < 0) {
return 1.0;
}
// 1. Compute the Wilson-Hilferty Normal Approximation Z-score
const term1 = Math.pow(chiSquareValue / df, 1 / 3);
const term2 = 1 - 2 / (9 * df);
const denominator = Math.sqrt(2 / (9 * df));
const zScore = (term1 - term2) / denominator;
// 2. Compute right-tailed probability from Z-score using standard Normal CDF
return standardNormalRightTail(zScore);
}
/**
* Standard Normal Right-Tail Probability Q(z) = P(Z > z)
* Uses a highly precise numerical approximation for the error function (erf).
*/
function standardNormalRightTail(z: number): number {
// Standard Normal Cumulative Distribution Function (CDF)
const cdf = 0.5 * (1 + errorFunction(z / Math.sqrt(2)));
// Right tail is 1 - CDF
return 1 - cdf;
}
/**
* Abramowitz and Stegun approximation for the Error Function (erf),
* accurate to 1.5 x 10^-7
*
* See https://en.wikipedia.org/wiki/Error_function#Bounds_and_numerical_approximations
*/
function errorFunction(x: number): number {
const sign = x < 0 ? -1 : 1;
const absX = Math.abs(x);
// Constants for the numerical approximation
const p = 0.3275911;
const a1 = 0.254829592;
const a2 = -0.284496736;
const a3 = 1.421413741;
const a4 = -1.453152027;
const a5 = 1.061405429;
const t = 1.0 / (1.0 + p * absX);
const y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.exp(-absX * absX);
return sign * y;
}
/**
* Calculates Cramér's V to determine the practical effect size (strength of bias).
* See https://en.wikipedia.org/wiki/Cramér%27s_V
*
* @param chiSquare The empirical chi-square value (e.g., 1190.372)
* @param totalSamples The total number of data points (e.g., 100000)
* @param numBuckets The total number of bins/buckets (e.g., 100)
* @returns A value between 0.0 (no bias) and 1.0 (extreme bias)
*/
function computeCramersV(chiSquare: number, totalSamples: number, numBuckets: number): number {
const df = numBuckets - 1;
if (df <= 0 || totalSamples <= 0) return 0;
// Core Formula: sqrt( chi2 / (N * df) )
return Math.sqrt(chiSquare / (totalSamples * df));
}
|
