/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* 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 */