/* nag_nonpar_gofstat_anddar (g08chc) Example Program.
*
* Copyright 2024 Numerical Algorithms Group.
*
* Mark 30.1, 2024.
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
#include <string.h>
int main(void) {
/* Scalars */
double a2, aa2, beta, nupper, p, sa2, sbeta;
const Integer lseed = 1, subid = -1;
Integer exit_status = 0, i, j, k, lstate = 17, n, nsim, n_pseudo;
/* Arrays */
double *x = 0, *xsim = 0, *y = 0;
Integer seed[1], state[17];
/* NAG types */
Nag_Boolean issort;
NagError fail;
printf("%s\n\n",
"nag_nonpar_gofstat_anddar (g08chc) Example Program Results");
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Read number of observations */
scanf("%" NAG_IFMT "", &n);
scanf("%*[^\n] ");
/* Memory allocation */
if (!(x = NAG_ALLOC(n, double)) || !(y = NAG_ALLOC(n, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Read observations */
for (i = 0; i < n; i++) {
scanf("%lf", x + i);
}
scanf("%*[^\n] ");
/* Maximum likelihood estimate of mean */
for (i = 0, beta = 0.0; i < n; i++) {
beta += x[i];
}
beta /= (double)n;
/* PIT, using exponential CDF with mean beta */
for (i = 0; i < n; i++) {
y[i] = 1.0 - exp(-x[i] / beta);
}
/* Let nag_nonpar_gofstat_anddar (g08chc) sort the (approximately)
uniform variates */
issort = Nag_FALSE;
/* Calculate the Anderson-Darling goodness-of-fit test statistic */
INIT_FAIL(fail);
/* nag_nonpar_gofstat_anddar (g08chc) */
a2 = nag_nonpar_gofstat_anddar(n, issort, y, &fail);
/* Correction due to estimated mean */
aa2 = (1.0 + 0.6 / (double)n) * a2;
/* Number of simulations; a suitably high number */
nsim = 888;
/* Generate exponential variates using a repeatable seed */
n_pseudo = n * nsim;
if (!(xsim = NAG_ALLOC(n_pseudo, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
INIT_FAIL(fail);
/* Initialize NAG's Basic pseudorandom number generator to give a
repeatable sequence */
seed[0] = 206033;
/* nag_rand_init_repeat (g05kfc) */
nag_rand_init_repeat(Nag_Basic, subid, (const Integer *)seed, lseed, state,
&lstate, &fail);
/* Generate a vector of pseudorandom numbers from an exponential
distribution */
/* nag_rand_dist_exp (g05sfc) */
nag_rand_dist_exp(n_pseudo, beta, state, xsim, &fail);
/* Simulations loop */
for (j = 0, nupper = 0.0; j < nsim; j++) {
/* Index in the simulated data */
k = j * n;
/* Maximum likelihood estimate of mean */
for (i = 0, sbeta = 0.0; i < n; i++) {
sbeta += xsim[k + i];
}
sbeta /= (double)n;
/* PIT */
for (i = 0; i < n; i++) {
y[i] = 1.0 - exp(-xsim[k + i] / sbeta);
}
/* Calculate A-squared */
/* nag_nonpar_gofstat_anddar (g08chc) */
sa2 = nag_nonpar_gofstat_anddar(n, issort, y, &fail);
if (sa2 > aa2) {
nupper++;
}
}
/* Simulated upper tail probability value */
p = nupper / (double)(nsim + 1);
/* Results */
printf("%s", " H0: data from exponential distribution with mean ");
printf("%g\n", beta);
printf("%s", " Test statistic, A-squared: ");
printf("%6g\n", a2);
printf("%s", " Upper tail probability: ");
printf("%6g\n", p);
END:
NAG_FREE(x);
NAG_FREE(xsim);
NAG_FREE(y);
return exit_status;
}