Тест на однородность
Первый тест самый простой - проверка на однородность. О нем мы уже говорили. Фактически случайные числа будут проверяться на равномерность распределения по диапазону от 0.0 до 1.0. Разобьем весь диапазон на 100 поддиапазонов, сформируем набор из 1000000 случайных чисел и вычислим количество значений, попавших в каждый поддиапазон. В поддиапазоне 0 будут находиться значения от 0.00 до 0.01, в поддиапазоне 1 - значения от 0.01 до 0.02 и т.д. Вероятность попадания случайного числа в любой поддиапазон составляет 0.01. Для полученного распределения вычислим значение параметра хи-квадрат и сравним его с данными для стандартного распределения хи-квадрат, находящимися в строке, для 99 степеней свободы.
Листинг 6.5. Тест на однородность
procedure UnifomityTest(RandGen : TtdBasePRNG;
var ChiSquare : double; var DegsFreedo : integer);
var
BucketNumber, i : integer;
Expected, ChiSqVal : double;
Bucket : array [0..pred(Uniformitylntervals) ] of integer;
begin
{вычислить количество чисел в каждом поддиапазоне}
FillChar(Bucket, sizeof(Bucket), 0);
for i := 0 to pred(UniformityCount) do
begin
BucketNumber := trunc(RandGen.AsDouble * Uniformitylntervals);
inc (Bucket [BucketNumber]);
end;
{вычислить значение параметра xu-квадрат}
Expected := UniformityCount / Uniformitylntervals;
ChiSqVal := 0.0;
for i := 0 to pred(Uniformitylntervals) do
ChiSqVal := ChiSqVal + (Sqr (Expected - Bucket [i]) / Expected);
{вернуть значения}
ChiSquare := ChiSqVal;
DegsFreedom := pred(Uniformitylntervals);
end;
Второй тест, который мы проведем, - тест на пропуски - несколько сложнее первого. Тест на пропуски гарантирует, что последовательность случайных чисел не будет попадать сначала в один поддиапазон, а затем в другой, третий и т.д., несмотря на то, что в целом значения будут распределены равномерно по всему диапазону. Определим в диапазоне поддиапазон, скажем, первую половину - от 0.0 до 0.5. Сформируем набор случайных чисел. Для каждого генерируемого числа будем проверять, попадает ли оно в выбранный поддиапазон (попадание) или нет (промах). В результате проверок будет получена последовательность попаданий и промахов. Найдите последовательности из одного и большего количества промахов (такие последовательности называются пропусками, отсюда и название теста - тест на пропуски). Вы получите последовательности из одного, двух и даже большего количества промахов. Разбейте длины пропусков на категории. Если известно, что вероятность попадания равна p (в нашем случае она будет равна длине выбранного поддиапазона), то вероятность промаха будет (1 -p). На основе этих данных можно определить вероятность возникновения пропуска из одного промаха — (1 -p)p, двух промахов — (1 -p)(^2^)p, n промахов - (1 -p)(^n^)p, а, следовательно, вычислить ожидаемое количество пропусков любой длины. После этого применим тест по критерию хи-квадрат. Будем использовать 10 категорий пропусков (поскольку вероятность возникновения пропусков длиной 11 и более промахов очень мала, все пропуски длиной 10 и более будут учитываться в последней категории;
при этом, конечно, следует учитывать реальную вероятность попадания длины пропуска в эту последнюю категорию), следовательно, мы получим девять степеней свободы. Как правило, тест на пропуски проводится пять раз: для первой и второй половины диапазона, а также для первой, второй и третьей третей диапазона.
Листинг 6.6. Тест на пропуски
procedure GapTest(RandGen : TtdBasePRNG;
Lower, Upper : double;
var ChiSquare : double;
var DegsFreedom : integer);
var
NumGaps : integer;
GapLen : integer;
i : integer;
p : double;
Expected : double;
ChiSqVal : double;
R : double;
Bucket : array [0..pred(GapBucketCount) ] of integer;
begin
{вычислить длины пропусков и определить количество пропусков в каждой категории}
FillChar(Bucket, sizeof(Bucket), 0);
GapLen := 0;
NumGaps := 0;
while (NumGaps < GapsCount) do
begin
R := RandGen.AsDouble;
if (Lower <= R) and (R < Upper) then begin
if (GapLen >= GapBucketCount) then
GapLen := pred(GapBucketCount);
inc(Bucket[GapLen]);
inc(NumGaps);
GapLen := 0;
end else
if (GapLen < GapBucketCount) then
inc(GapLen);
end;
p := Upper - Lower;
ChiSqVal := 0.0;
{обработать все категории, кроме последней}
for i := 0 to GapBucketCount-2 do
begin
Expected := p * IntPower(1-p, i) * NumGaps;
ChiSqVal := ChiSqVal + (Sqr (Expected - Bucket [i]) / Expected);
end;
{обработать последнюю категорию}
i := pred(GapBucketCount);
Expected IntPower (1-p, i) * NumGaps;
ChiSqVal := ChiSqVal + (Sqr (Expected - Bucket [i]) / Expected);
{вернуть значения}
ChiSquare := ChiSqVal;
DegsFreedom := pred(GapBucketCount);
end;
Третий тест известен под названием "покер" (poker test). Случайные числа группируются в наборы по пять, а затем преобразуются в "карты", которые представляют собой цифры от 0 до 9. После этого определяется количество разных карт в каждом наборе (оно будет равно от одного до пяти) и полученные результаты разбиваются на категории. Поскольку вероятность пятикратного повторения одной и той же карты достаточно низка, случай выпадения только одной карты, как правило, включается в категорию "две разные цифры". К полученным четырем категориям применятся тест по критерию хи-квадрат (три степени свободы). Вероятность возникновения события для каждой категории вычислить не так уж легко (к тому же математические выкладки основаны на использовании комбинаторных значений, называемых числами Стерлинга), поэтому вычисления в этой книге не приводятся. Если вам интересно, то подробное описание можно найти в [11].
Листинг 6.7. Тест "покер"
procedure PokerTest(RandGen : TtdBasePRNG;
var ChiSquare : double;
var DegsFreedom : integer);
var
i, j, jlBucketNumber, NumFives : integer;
Accum, Divisor, Expected, ChiSqVal : double;
Bucket : array [0..4] of integer;
Flag : array [0..9] of boolean;
p : array [0..4] of double;
begin
{подготовительные операции}
FillChar(Bucket, sizeof(Bucket), 0);
NumFives PokerCount div 5;
{вычислить вероятности для каждой категории событий, алгоритм Кнута}
Accum := 1.0;
Divisor := IntPower(10.0, 5);
for i := 0 to 4 do
begin
Accum := Accum * (10.0 - i);
p[i] := Accum * Stirling(5, succ(i)) / Divisor;
end;
{для каждой группы из пяти случайных чисел преобразовать все значения и числа от 1 до 10, определить количество разных цифр}
for i := 1 to NumFives do
begin
FillChar(Flag, sizeof(Flag), 0);
for j := 1 to 5 do begin
Flag [trunc(RandGen.AsDouble * 10.0)] :=true;
end;
BucketNumber := -1;
for j := 0 to 9 do
if Flag[j] then
inc(BucketNumber);
inc(Bucket[BucketNumber]);
end;
{объединить две первые категории - это будет сумма категорий "все цифры одинаковы" и "две разные цифры"}
inc(Bucket[1], Bucket[0]);
Expected := (p[0]+p[1]) * NumFives;
ChiSqVal := Sqr(Expected - Bucket[1]) / Expected;
{обработать другие категории}
for i := 2 to 4 do
begin
Expected :=p[i] * NumFives;
ChiSqVal := ChiSqVal + (Sqr (Expected - Bucket [i]) / Expected);
end;
{вернуть значения}
ChiSquare := ChiSqVal;
DegsFreedom := 3;
end;
Четвертый тест, который мы будем проводить, называется "сбор купонов" (coupon collector's test). Случайные числа считываются по одному и преобразуются в "купоны" - числа от 0 до 4. Фиксируется длина последовательности до получения полного комплекта купонов (т.е. цифр от 0 до 4). Очевидно, что получаемые длины будут от пяти и выше. После набора полного комплекта сбор купонов начинается снова. Длины последовательностей разбиваются на категории, к которым затем применяется тест по критерию хи-квадрат. Как правило, используются категории для длин от 5 до 19 и еще одна дополнительная категория для больших длин. Таким образом, мы получаем 16 категорий, а, следовательно, 15 степеней свободы. Как и в тесте "покер", вычисление вероятностей для каждой из категорий включает сложные математические выкладки, которые в этой книге не приводятся. Соответствующие подробности можно найти в [11].
Листинг 6.8. Тест "сбор купонов"
procedure CouponCollectorsTest(RandGen : TtdBasePRNG;
var ChiSquare : double;
var DegsFreedom : integer);
var
NumSeqs, LenSeq, NumVals, NewVal, i : integer;
Expected, ChiSqVal : double;
Bucket : array [5..20] of integer;
Occurs : array [0..4] of boolean;
p : array [5..20] of double;
begin
{вычислить вероятности для каждой категории событий, алгоритм Кнута}
p[20] := 1.0;
for i := 5 to 19 do
begin
p[i] := (120.0 * Stirling(i-1, 4)) / IntPower(5.0, i);
p[20] := p[20] - p[i];
end;
NumSeqs := 0;
FillChar(Bucket, sizeof(Bucket), 0);
while (NumSeqs < CouponCount) do
begin
{продолжать сбор купонов (т.е. случайных чисел) до получения полного набора из пяти купонов}
LenSeq := 0;
NumVals := 0;
FillChar (Occurs, sizeof(Occurs), 0);
repeat
inc(LenSeq);
NewVal := trune(RandGen.AsDouble * 5);
if not Occurs [NewVal] then begin
Occurs[NewVal] := true;
inc(NumVals);
end;
until (NumVals = 5);
{обновить значение для соответствующей категории в зависимости от количества собранных купонов}
if (LenSeq > 20) then
LenSeq := 20;
inc(Bucket[LenSeq]);
inc (NumSeqs);
end;
{вычислить значение xu-квадрат}
ChiSqVal := 0.0;
for i := 5 to 20 do
begin
Expected := p [ i ] * NumSeqs;
ChiSqVal := ChiSqVal + (Sqr(Expected - Bucket [i]) / Expected);
end;
{вернуть значения}
ChiSquare := ChiSqVal;
DegsFreedom := 15;
end;
Результаты выполнения тестов
В разделе сопровождающих эту книгу материалов, который расположен на Web-сайте издательства, можно найти тестовую программу, которая применяет все рассмотренные нами тесты к стандартному генератору случайных чисел Delphi и минимальному стандартному генератору случайных чисел. На рис. 6.1 приведены результаты проведения одного из тестов для генератора случайных чисел Delphi.
