Testing Random Number Generators

Listing 5: The serial correlation test

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  */
/* Repeat Step 1 through Step 4 N times, then proceed to Step 5:    */
/*                                                                  */
/*   Step 1: Generate T variates, Xk, 0 <= Xk <= 32767, 1 <= k <= T.*/
/*                                                                  */
/*   Step 2: Calculate NumCorr correlation coefficients as follows: */
/*                                                                  */
/*           err = bjcorr(X, T, Corr, NumCorr);                     */
/*                                                                  */
/*           where                                                  */
/*                                                                  */
/*               Corr      receiving area for correlation           */
/*             coefficients                                         */
/*               NumCorr   User-specified number <= T - 3           */
/*                                                                  */
/*   Step 3: Calculate a test statistic, Ut:                        */
/*                                                                  */
/*           Ut = T * (T + 2) * Sum{Corr[j]^2 / (T-j), j = 1, J}    */
/*                                                                  */
/*           where                                                  */
/*               J = NumCorr                                        */
/*                                                                  */
/*                                                                  */
/*   Step 4: Calculate a chi-square probability (Pt) associated     */
/*      with Ut, with J degrees of freedom ; save it in an          */
/*      array of size N,                                            */
/*                                                                  */
/*           When N probabilities, Pt, have been computed,          */
/*           go  to  Step  5.                                       */
/*                                                                  */
/*   Step 5: Run a KS analysis on Pt(j), 1 <= j <= N.               */
/*                                                                  */
/*  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mconf.h>
#include <srcrdefs.h>
#include <miscdefs.h>

#define LOW_PROB        1.0e-8

int     Variates[32768u];
double  CorCoefs[32768u];
/* ============================================================== */
/* CalcSerCorChiSq - Gets Chi-square p-value for NumCorr          */
/* Coefficients                                                   */
/* ============================================================== */
void
CalcSerCorChiSq(SERCOR_DATA_STRU * SerCorData)
{
    UINT    k;
    UINT    NumCoefs = SerCorData->NumCoefs,
            NumVar   = SerCorData->NumVarPerPass;

    double  ChiSqProb, ChiSqStat;
    double  CoefSos = 0, CoefSum = 0, RmsCoefs, Sum = 0;
    double  LoLimit, HiLimit;

    /* ----------------------------------------- */
    /* Step 1: Generate Array of Random Deviates */
    /* ----------------------------------------- */
    for (k = 0; k < NumVar; ++k)
    {
        Variates[k] = SerCorData->RandFun();
    }

    /* ------------------------------------------ */
    /* Step 2: Calculate Correlation Coefficients */
    /* ------------------------------------------ */
    SerCorData->CallStatus =
        BJCoef(Variates, NumVar, CorCoefs, NumCoefs);

    if (SerCorData->CallStatus == OK)
    {
        /* ---------------------------------------------------- */
        /* Calculate Limits on Acceptable Chi-Square Statistics */
        /* ---------------------------------------------------- */
        ChiSqDist(LOW_PROB, NumCoefs, &LoLimit, &HiLimit);

        /* -------------------------------- */
        /* Step 3: Calculate Test Statistic */
        /* -------------------------------- */
        for (k = 0; k < NumCoefs; ++k)
        {
            double  SqrCoef = SQR(CorCoefs[k]);

            /* ------------------- */
            /* Sum of Coefficients */
            /* ------------------- */
            CoefSum += CorCoefs[k];

            /* ------------------------------ */
            /* Sum of Squares of Coefficients */
            /* ------------------------------ */
            CoefSos += SqrCoef;

            /* ------------------------- */
            /* Accumulate Test Statistic */
            /* ------------------------- */
            Sum += SqrCoef / (double) (NumVar - (k+1));
        }

        /* ----------------------------- */
        /* Finish Step 3: Test Statistic */
        /* ----------------------------- */
        ChiSqStat = Sum * (double)NumVar * (double)(NumVar + 2);

        /* ---------------------------------------- */
        /* Step 4: Calculate Chi-Square Probability */
        /* ---------------------------------------- */
        if (ChiSqStat < LoLimit)
        {
              ChiSqProb = LOW_PROB;
        }       else if (ChiSqStat > HiLimit)
        {
              ChiSqProb = 1.0 - LOW_PROB;
        }
        else
        {
             ChiSqProb = chdtr(NumCoefs, ChiSqStat);
        }
        SerCorData->ChiSqProb = ChiSqProb;

        /* --------------------------------------------------- */
        /* Maintain Counts of Numbers of Variates and Outliers */
        /* --------------------------------------------------- */
        SerCorData->TotCoefs += NumCoefs;
        SerCorData->TotSos   += CoefSos;
        SerCorData->TotSum   += CoefSum;

        RmsCoefs = sqrt(CoefSos / (double)NumCoefs);

        for (k = 0; k < NumCoefs; ++k)
        {
            if (fabs(CorCoefs[k]) > 2*RmsCoefs)
            {
                ++SerCorData->NumOutliers;
            }
        }
    }
                                  /* End of if */
}
/* End of File */