1 /*************************************************************************\
2 *
3 * Package: TestU01
4 * File: svaria.c
5 * Environment: ANSI C
6 *
7 * Copyright (c) 2002 Pierre L'Ecuyer, DIRO, Université de Montréal.
8 * e-mail: lecuyer@iro.umontreal.ca
9 * All rights reserved.
10 *
11 * Redistribution and use in source and binary forms, with or without
12 * modification, are permitted without a fee for private, research,
13 * academic, or other non-commercial purposes.
14 * Any use of this software in a commercial environment requires a
15 * written licence from the copyright owner.
16 *
17 * Any changes made to this package must be clearly identified as such.
18 *
19 * In scientific publications which used this software, a reference to it
20 * would be appreciated.
21 *
22 * Redistributions of source code must retain this copyright notice
23 * and the following disclaimer.
24 *
25 * THIS PACKAGE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR
26 * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
27 * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
28 *
29 \*************************************************************************/
30
31 #include "gdef.h"
32 #include "util.h"
33 #include "tables.h"
34 #include "chrono.h"
35 #include "num.h"
36 #include "num2.h"
37
38 #include "svaria.h"
39 #include "unif01.h"
40 #include "sres.h"
41 #include "wdist.h"
42 #include "swrite.h"
43 #include "smultin.h"
44
45 #include "fmass.h"
46 #include "gofs.h"
47 #include "gofw.h"
48
49 #include <math.h>
50 #include <float.h>
51 #include <limits.h>
52 #include <stdio.h>
53 #include <string.h>
54
55
56
57
58 lebool svaria_Timer = FALSE;
59
60
61
62
63 /*-------------------------------- Constants ------------------------------*/
64
65 /* Max string lengths */
66 #define LEN1 100
67 #define LEN2 200
68
69 /* Sample limit for normal approximation in svaria_SampleMean */
70 #define SAM_LIM 60
71
72 /* Arrays dimension in svaria_AppearanceSpacings */
73 #define AS_DIM 32
74
75 /* The Meschach package for matrix computations */
76 #undef MESCHACH
77
78
79
80
81
82 /*-------------------------------- Functions ------------------------------*/
83
84
InitFDistMeans(int n,double Coef[])85 static void InitFDistMeans (int n, double Coef[])
86 /*
87 * Initializes the distribution for svaria_SampleMean by computing the
88 * coefficients. We shall keep the value of n in element Coef[SAM_LIM],
89 * since we shall need it to get the value of the distribution at x.
90 */
91 {
92 int s;
93 double z;
94 fmass_INFO Q;
95
96 z = num2_Factorial (n);
97 /* This uses the binomial formulae, but is not the binomial probability
98 distribution since p + q != 1 */
99 Q = fmass_CreateBinomial (n, -1.0, 1.0);
100 for (s = 0; s <= n; s++)
101 Coef[s] = fmass_BinomialTerm2 (Q, s) / z;
102 fmass_DeleteBinomial (Q);
103 Coef[SAM_LIM] = n;
104
105 if (swrite_Classes) {
106 printf ("---------------------------------------\n");
107 for (s = 0; s <= n; s++) {
108 printf (" Coeff[%2d] = %14.6g\n", s, Coef[s]);
109 }
110 printf ("\n");
111 }
112 }
113
114 /*-------------------------------------------------------------------------*/
115
FDistMeans(double C[],double x)116 static double FDistMeans (
117 double C[], /* Coefficients and sample size n */
118 double x /* Argument */
119 )
120 /*
121 * Distribution function of sample mean as in Stephens (1966), p.235.
122 * This function is not very precise: the normal approximation is poor
123 * for small x, and the computation of the exact function is numerically
124 * unstable for large n. This will be used only for n < SAM_LIM. The
125 * value of n is in C[SAM_LIM].
126 */
127 {
128 double Sum;
129 int M;
130 int i;
131 double nLR = C[SAM_LIM];
132 int n = nLR;
133
134 if (x <= 0.0)
135 return 0.0;
136 if (x >= n)
137 return 1.0;
138
139 M = x;
140 Sum = 0.0;
141 if (x < n / 2.0) {
142 for (i = 0; i <= M; i++) {
143 Sum += C[i] * pow (x, nLR);
144 x -= 1.0;
145 }
146 } else {
147 x = -x + nLR;
148 for (i = n; i >= M + 1; i--) {
149 Sum += C[i] * pow (x, nLR);
150 x -= 1.0;
151 }
152 if (!(n & 1))
153 Sum = -Sum;
154 Sum += 1.0;
155 }
156 return Sum;
157 }
158
159 /*-------------------------------------------------------------------------*/
160
svaria_SampleMean(unif01_Gen * gen,sres_Basic * res,long N,long n,int r)161 void svaria_SampleMean (unif01_Gen * gen, sres_Basic * res,
162 long N, long n, int r)
163 {
164 long i;
165 long Seq;
166 double Sum;
167 double Coef[SAM_LIM + 1];
168 lebool localRes = FALSE;
169 chrono_Chrono *Timer;
170 char *TestName = "svaria_SampleMean test";
171
172 Timer = chrono_Create ();
173 if (swrite_Basic) {
174 swrite_Head (gen, TestName, N, n, r);
175 printf ("\n\n");
176 }
177 util_Assert (n > 1, "svaria_SampleMean: n < 2");
178
179 if (res == NULL) {
180 localRes = TRUE;
181 res = sres_CreateBasic ();
182 }
183 sres_InitBasic (res, N, "svaria_SampleMean");
184 if (n < SAM_LIM)
185 InitFDistMeans (n, Coef);
186
187 if (n < SAM_LIM)
188 statcoll_SetDesc (res->sVal1, "SampleMean sVal1: n*<U>");
189 else
190 statcoll_SetDesc (res->sVal1, "SampleMean sVal1: standard normal");
191
192 for (Seq = 1; Seq <= N; Seq++) {
193 Sum = 0.0;
194 for (i = 1; i <= n; i++)
195 Sum += unif01_StripD (gen, r);
196
197 if (n < SAM_LIM)
198 statcoll_AddObs (res->sVal1, Sum);
199 else
200 statcoll_AddObs (res->sVal1, sqrt (12.0 / n) * (Sum - 0.5 * n));
201 }
202
203 if (n < SAM_LIM) {
204 gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, FDistMeans, Coef,
205 res->sVal2, res->pVal2);
206 } else {
207 /* Normal approximation */
208 gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_Normal,
209 (double *) NULL, res->sVal2, res->pVal2);
210 }
211 res->pVal1->NObs = N;
212
213 if (swrite_Collectors)
214 statcoll_Write (res->sVal1, 5, 14, 4, 3);
215
216 if (swrite_Basic) {
217 gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2,
218 "Statistic value :");
219 swrite_Final (gen, Timer);
220 }
221 if (localRes)
222 sres_DeleteBasic (res);
223 chrono_Delete (Timer);
224 }
225
226
227 /*=========================================================================*/
228
svaria_SampleCorr(unif01_Gen * gen,sres_Basic * res,long N,long n,int r,int k)229 void svaria_SampleCorr (unif01_Gen * gen, sres_Basic * res,
230 long N, long n, int r, int k)
231 {
232 long i;
233 long Seq;
234 double U;
235 double Sum;
236 double *Pre; /* Previous k generated numbers */
237 int pos; /* Circular index to element at lag k */
238 lebool localRes = FALSE;
239 chrono_Chrono *Timer;
240 char *TestName = "svaria_SampleCorr test";
241
242 Timer = chrono_Create ();
243 if (swrite_Basic) {
244 swrite_Head (gen, TestName, N, n, r);
245 printf (", k = %d\n\n", k);
246 }
247 util_Assert (n > 2, "svaria_SampleCorr: n <= 2");
248
249 if (res == NULL) {
250 localRes = TRUE;
251 res = sres_CreateBasic ();
252 }
253 sres_InitBasic (res, N, "svaria_SampleCorr");
254 statcoll_SetDesc (res->sVal1,
255 "SampleCorr sVal1: asymptotic standard normal");
256
257 Pre = util_Calloc ((size_t) (k + 1), sizeof (double));
258
259 for (Seq = 1; Seq <= N; Seq++) {
260 /* Generate first k numbers U and keep them in Pre */
261 for (i = 0; i < k; i++)
262 Pre[i] = unif01_StripD (gen, r);
263
264 Sum = 0.0;
265 pos = 0;
266 /* Element Pre[pos] is at lag k from U */
267 for (i = k; i < n; i++) {
268 U = unif01_StripD (gen, r);
269 Sum += Pre[pos] * U - 0.25;
270 Pre[pos] = U;
271 pos++;
272 pos %= k;
273 }
274 /* Save standardized correlation */
275 statcoll_AddObs (res->sVal1, Sum * sqrt (12.0 / (n - k)));
276 }
277
278 gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_Normal,
279 (double *) NULL, res->sVal2, res->pVal2);
280 res->pVal1->NObs = N;
281 sres_GetNormalSumStat (res);
282
283 if (swrite_Collectors)
284 statcoll_Write (res->sVal1, 5, 14, 4, 3);
285
286 if (swrite_Basic) {
287 gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2,
288 "Normal statistic :");
289 swrite_NormalSumTest (N, res);
290 swrite_Final (gen, Timer);
291 }
292 util_Free (Pre);
293 if (localRes)
294 sres_DeleteBasic (res);
295 chrono_Delete (Timer);
296 }
297
298
299 /*=========================================================================*/
300
FDistProd(double Par[],double x)301 static double FDistProd (
302 double Par[], /* The parameter t = Par[0] */
303 double x /* The argument */
304 )
305 /*
306 * Compute the distribution function F for the product of t random variables
307 * U[0, 1], where
308 * t - 1
309 * __ j
310 * F[u1*u2*...ut <= x] = x \ (-ln(x))
311 * /__ -----------
312 * j = 0 j!
313 *
314 */
315 {
316 double vlog, jterm, vterm, Sum;
317 int j, t;
318
319 if (x >= 1.0)
320 return 1.0;
321 if (x <= 0.0)
322 return 0.0;
323
324 vlog = log (x);
325 t = Par[0];
326 Sum = 1.0;
327 vterm = 1.0;
328 jterm = 1.0;
329 for (j = 1; j < t; j++) {
330 vterm *= vlog;
331 jterm *= -j;
332 Sum += vterm / jterm;
333 if (vterm / jterm < DBL_EPSILON)
334 break;
335 }
336 return x * Sum;
337 }
338
339
340 /*-------------------------------------------------------------------------*/
341
svaria_SampleProd(unif01_Gen * gen,sres_Basic * res,long N,long n,int r,int t)342 void svaria_SampleProd (unif01_Gen * gen, sres_Basic * res,
343 long N, long n, int r, int t)
344 {
345 long i;
346 int j;
347 long Seq;
348 double *P;
349 double temp;
350 double Par[1];
351 lebool localRes = FALSE;
352 chrono_Chrono *Timer;
353 char *TestName = "svaria_SampleProd test";
354
355 Timer = chrono_Create ();
356 if (swrite_Basic) {
357 swrite_Head (gen, TestName, N, n, r);
358 printf (", t = %d\n\n", t);
359 }
360
361 if (res == NULL) {
362 localRes = TRUE;
363 res = sres_CreateBasic ();
364 }
365 sres_InitBasic (res, N, "svaria_SampleProd");
366
367 P = util_Calloc ((size_t) n + 1, sizeof (double));
368 statcoll_SetDesc (res->sVal1, "SampleProd sVal1: Uniform [0, 1]");
369 Par[0] = t;
370
371 for (Seq = 1; Seq <= N; Seq++) {
372 for (i = 1; i <= n; i++) {
373 temp = unif01_StripD (gen, r);
374 for (j = 2; j <= t; j++)
375 temp *= unif01_StripD (gen, r);
376 P[i] = temp;
377 }
378 gofw_ActiveTests1 (P, n, FDistProd, Par, res->sVal2, res->pVal2);
379 statcoll_AddObs (res->sVal1, res->pVal2[gofw_AD]);
380 }
381
382 gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_Unif,
383 (double *) NULL, res->sVal2, res->pVal2);
384 res->pVal1->NObs = N;
385
386 if (swrite_Collectors)
387 statcoll_Write (res->sVal1, 5, 14, 4, 3);
388
389 if (swrite_Basic) {
390 gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2,
391 "Anderson-Darling statistic :");
392 swrite_Final (gen, Timer);
393 }
394 util_Free (P);
395 if (localRes)
396 sres_DeleteBasic (res);
397 chrono_Delete (Timer);
398 }
399
400
401 /*=========================================================================*/
402
svaria_SumLogs(unif01_Gen * gen,sres_Chi2 * res,long N,long n,int r)403 void svaria_SumLogs (unif01_Gen * gen, sres_Chi2 * res,
404 long N, long n, int r)
405 {
406 const double Eps = DBL_EPSILON / 2.0; /* To avoid log(0) */
407 const double Epsilon = 1.E-100; /* To avoid underflow */
408 long i;
409 long Seq;
410 double u;
411 double Prod;
412 double Sum;
413 double V[1];
414 lebool localRes = FALSE;
415 chrono_Chrono *Timer;
416 char *TestName = "svaria_SumLogs test";
417 char chaine[LEN1 + 1] = "";
418 char str[LEN2 + 1];
419
420 Timer = chrono_Create ();
421 if (swrite_Basic) {
422 swrite_Head (gen, TestName, N, n, r);
423 printf ("\n\n");
424 }
425 util_Assert (n < LONG_MAX/2, "2n > LONG_MAX");
426 if (res == NULL) {
427 localRes = TRUE;
428 res = sres_CreateChi2 ();
429 }
430 sres_InitChi2 (res, N, -1, "svaria_SumLogs");
431
432 strncpy (chaine, "SumLogs sVal1: chi2 with ", (size_t) LEN1);
433 sprintf (str, "%ld", 2 * n);
434 strncat (chaine, str, (size_t) LEN2);
435 strncat (chaine, " degrees of freedom", (size_t) LEN1);
436 statcoll_SetDesc (res->sVal1, chaine);
437 res->degFree = 2 * n;
438 if (res->degFree < 1) {
439 util_Warning (TRUE, "Chi-square with 0 degree of freedom.");
440 if (localRes)
441 sres_DeleteChi2 (res);
442 return;
443 }
444
445 for (Seq = 1; Seq <= N; Seq++) {
446 Prod = 1.0;
447 Sum = 0.0;
448 for (i = 1; i <= n; i++) {
449 u = unif01_StripD (gen, r);
450 if (u < Eps)
451 u = Eps;
452 Prod *= u;
453 if (Prod < Epsilon) {
454 Sum += log (Prod);
455 Prod = 1.0;
456 }
457 }
458 statcoll_AddObs (res->sVal1, -2.0 * (Sum + log (Prod)));
459
460 }
461 V[0] = 2 * n;
462 gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_ChiSquare, V,
463 res->sVal2, res->pVal2);
464 res->pVal1->NObs = N;
465 sres_GetChi2SumStat (res);
466
467 if (swrite_Collectors)
468 statcoll_Write (res->sVal1, 5, 14, 4, 3);
469
470 if (swrite_Basic) {
471 swrite_AddStrChi (str, LEN2, res->degFree);
472 gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2, str);
473 swrite_Chi2SumTest (N, res);
474 swrite_Final (gen, Timer);
475 }
476 if (localRes)
477 sres_DeleteChi2 (res);
478 chrono_Delete (Timer);
479 }
480
481
482 /*=========================================================================*/
483
WriteDataWeight(unif01_Gen * gen,char * TestName,long N,long n,int r,long k,double Alpha,double Beta)484 static void WriteDataWeight (unif01_Gen * gen, char *TestName,
485 long N, long n, int r, long k, double Alpha, double Beta)
486 {
487 swrite_Head (gen, TestName, N, n, r);
488 printf (", k = %1ld, Alpha = %6.4g, Beta = %6.4g\n\n",
489 k, Alpha, Beta);
490 }
491
492
493 /*-------------------------------------------------------------------------*/
494
svaria_WeightDistrib(unif01_Gen * gen,sres_Chi2 * res,long N,long n,int r,long k,double Alpha,double Beta)495 void svaria_WeightDistrib (unif01_Gen * gen, sres_Chi2 * res,
496 long N, long n, int r, long k, double Alpha, double Beta)
497 {
498 long W;
499 long j;
500 long i;
501 long Seq;
502 double X;
503 double U;
504 double p;
505 double nLR = n;
506 double V[1];
507 long NbClasses;
508 long *Loc;
509 fmass_INFO Q;
510 lebool localRes = FALSE;
511 chrono_Chrono *Timer;
512 char *TestName = "svaria_WeightDistrib test";
513 char chaine[LEN1 + 1] = "";
514 char str[LEN2 + 1];
515
516 Timer = chrono_Create ();
517 if (swrite_Basic)
518 WriteDataWeight (gen, TestName, N, n, r, k, Alpha, Beta);
519
520 /* util_Assert (n >= 3.0 * gofs_MinExpected,
521 "svaria_WeightDistrib: n is too small"); */
522 util_Assert (Alpha <= 1.0 && Alpha >= 0.0,
523 "svaria_WeightDistrib: Alpha must be in [0, 1]");
524 util_Assert (Beta <= 1.0 && Beta >= 0.0,
525 "svaria_WeightDistrib: Beta must be in [0, 1]");
526 p = Beta - Alpha;
527
528 if (res == NULL) {
529 localRes = TRUE;
530 res = sres_CreateChi2 ();
531 }
532 sres_InitChi2 (res, N, k, "svaria_WeightDistrib");
533 Loc = res->Loc;
534
535 /* Compute binomial probabilities and multiply by n */
536 Q = fmass_CreateBinomial (k, p, 1.0 - p);
537 for (i = 0; i <= k; i++)
538 res->NbExp[i] = nLR * fmass_BinomialTerm2 (Q, i);
539 fmass_DeleteBinomial (Q);
540
541 res->jmin = 0;
542 res->jmax = k;
543 if (swrite_Classes)
544 gofs_WriteClasses (res->NbExp, Loc, res->jmin, res->jmax, 0);
545
546 /* Merge classes for the chi-square */
547 gofs_MergeClasses (res->NbExp, Loc, &res->jmin, &res->jmax, &NbClasses);
548
549 if (swrite_Classes)
550 gofs_WriteClasses (res->NbExp, Loc, res->jmin, res->jmax, NbClasses);
551
552 strncpy (chaine, "WeightDistrib sVal1: chi2 with ", (size_t) LEN1);
553 sprintf (str, "%ld", NbClasses - 1);
554 strncat (chaine, str, (size_t) LEN2);
555 strncat (chaine, " degrees of freedom", (size_t) LEN1);
556 statcoll_SetDesc (res->sVal1, chaine);
557 res->degFree = NbClasses - 1;
558 if (res->degFree < 1) {
559 if (localRes)
560 sres_DeleteChi2 (res);
561 return;
562 }
563
564 for (Seq = 1; Seq <= N; Seq++) {
565 for (i = 0; i <= k; i++)
566 res->Count[i] = 0;
567 for (i = 1; i <= n; i++) {
568 W = 0;
569 for (j = 1; j <= k; j++) {
570 U = unif01_StripD (gen, r);
571 if (U >= Alpha && U < Beta)
572 ++W;
573 }
574 if (W > res->jmax)
575 ++res->Count[res->jmax];
576 else
577 ++res->Count[Loc[W]];
578 }
579 if (swrite_Counters)
580 tables_WriteTabL (res->Count, res->jmin, res->jmax, 5, 10,
581 "Observed numbers:");
582
583 X = gofs_Chi2 (res->NbExp, res->Count, res->jmin, res->jmax);
584 statcoll_AddObs (res->sVal1, X);
585 }
586
587 V[0] = NbClasses - 1;
588 gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_ChiSquare, V,
589 res->sVal2, res->pVal2);
590 res->pVal1->NObs = N;
591 sres_GetChi2SumStat (res);
592
593 if (swrite_Collectors)
594 statcoll_Write (res->sVal1, 5, 14, 4, 3);
595
596 if (swrite_Basic) {
597 swrite_AddStrChi (str, LEN2, res->degFree);
598 gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2, str);
599 swrite_Chi2SumTest (N, res);
600 swrite_Final (gen, Timer);
601 }
602 if (localRes)
603 sres_DeleteChi2 (res);
604 chrono_Delete (Timer);
605 }
606
607
608 /*=========================================================================*/
609
WriteDataArgMax(unif01_Gen * gen,char * TestName,long N,long n,int r,long k,long m)610 static void WriteDataArgMax (unif01_Gen * gen, char *TestName,
611 long N, long n, int r, long k, long m)
612 {
613 double x;
614
615 swrite_Head (gen, TestName, N, n, r);
616 printf (", k = %1ld, m = %1ld\n\n", k, m);
617 printf (" Number of balls = n = %1ld\n", n);
618 printf (" Number of urns = k = %1ld\n", k);
619
620 x = n;
621 x = x * x / (2 * k);
622 printf (" Number (approx) of collisions = n^2 / 2k = %g\n\n\n", x);
623 }
624
625 /*-------------------------------------------------------------------------*/
626
svaria_CollisionArgMax_00(unif01_Gen * gen,sres_Chi2 * res,long N,long n,int r,long k,long m)627 static int svaria_CollisionArgMax_00 (unif01_Gen *gen, sres_Chi2 *res,
628 long N, long n, int r, long k, long m)
629 /*
630 * Return 0 if no error, otherwise return != 0.
631 */
632 {
633 double X;
634 double U;
635 double Max;
636 long NbColl;
637 long Indice = -1;
638 long j;
639 long i;
640 long Rep;
641 long Seq;
642 long NbClasses;
643 long *Loc;
644 int *Urne;
645 double V[1];
646 fmass_INFO Q;
647 lebool localRes = FALSE;
648 chrono_Chrono *chro, *Timer;
649 char *TestName = "svaria_CollisionArgMax test";
650 char chaine[LEN1 + 1] = "";
651 char str[LEN2 + 1];
652
653 Timer = chrono_Create ();
654 if (swrite_Basic)
655 WriteDataArgMax (gen, TestName, N, n, r, k, m);
656
657 util_Assert (n <= 4 * k, "svaria_CollisionArgMax: n > 4k");
658 /* util_Assert (m > 2.0 * gofs_MinExpected,
659 "svaria_CollisionArgMax: m <= 2*gofs_MinExpected"); */
660
661 if (res == NULL) {
662 localRes = TRUE;
663 res = sres_CreateChi2 ();
664 }
665 sres_InitChi2 (res, N, n, "svaria_CollisionArgMax");
666 Loc = res->Loc;
667 Urne = util_Calloc ((size_t) k + 1, sizeof (int));
668
669 if (svaria_Timer) {
670 printf ("-----------------------------------------------");
671 printf ("\nCPU time to initialize the collision distribution: ");
672 chro = chrono_Create ();
673 }
674 Q = smultin_CreateCollisions (n, (smultin_CellType) k);
675 if (svaria_Timer) {
676 chrono_Write (chro, chrono_hms);
677 printf ("\n\n");
678 }
679
680 /* Compute the expected numbers of collisions: m*P(j) */
681 for (j = 0; j <= n; j++)
682 res->NbExp[j] = m * smultin_CollisionsTerm (Q, j);
683 smultin_DeleteCollisions (Q);
684
685 res->jmin = 0;
686 res->jmax = n;
687 if (swrite_Classes)
688 gofs_WriteClasses (res->NbExp, Loc, res->jmin, res->jmax, 0);
689
690 gofs_MergeClasses (res->NbExp, Loc, &res->jmin, &res->jmax, &NbClasses);
691
692 if (swrite_Classes)
693 gofs_WriteClasses (res->NbExp, Loc, res->jmin, res->jmax, NbClasses);
694
695 strncpy (chaine, "CollisionArgMax sVal1: chi2 with ", (size_t) LEN1);
696 sprintf (str, "%ld", NbClasses - 1);
697 strncat (chaine, str, (size_t) LEN2);
698 strncat (chaine, " degrees of freedom", (size_t) LEN1);
699 statcoll_SetDesc (res->sVal1, chaine);
700 res->degFree = NbClasses - 1;
701 if (res->degFree < 1) {
702 if (localRes)
703 sres_DeleteChi2 (res);
704 return 1;
705 }
706
707 if (svaria_Timer)
708 chrono_Init (chro);
709
710 for (Seq = 1; Seq <= N; Seq++) {
711 for (j = 0; j <= n; j++)
712 res->Count[j] = 0;
713
714 for (Rep = 1; Rep <= m; Rep++) {
715 for (j = 0; j <= k; j++)
716 Urne[j] = -1;
717
718 NbColl = 0;
719 for (j = 1; j <= n; j++) {
720 Max = -1.0;
721 for (i = 1; i <= k; i++) {
722 U = unif01_StripD (gen, r);
723 if (U > Max) {
724 Max = U;
725 Indice = i;
726 }
727 }
728 if (Urne[Indice] < 0)
729 Urne[Indice] = 1;
730 else
731 ++NbColl;
732 }
733 if (NbColl > res->jmax)
734 ++res->Count[res->jmax];
735 else
736 ++res->Count[Loc[NbColl]];
737 }
738 if (swrite_Counters)
739 tables_WriteTabL (res->Count, res->jmin, res->jmax, 5, 10,
740 "Observed numbers:");
741 X = gofs_Chi2 (res->NbExp, res->Count, res->jmin, res->jmax);
742 statcoll_AddObs (res->sVal1, X);
743 }
744
745 if (svaria_Timer) {
746 printf ("\n----------------------------------------------\n"
747 "CPU time for the test : ");
748 chrono_Write (chro, chrono_hms);
749 printf ("\n\n");
750 chrono_Delete (chro);
751 }
752
753 V[0] = NbClasses - 1;
754 gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_ChiSquare, V,
755 res->sVal2, res->pVal2);
756 res->pVal1->NObs = N;
757 sres_GetChi2SumStat (res);
758
759 if (swrite_Collectors)
760 statcoll_Write (res->sVal1, 5, 14, 4, 3);
761
762 if (swrite_Basic) {
763 swrite_AddStrChi (str, LEN2, res->degFree);
764 gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2, str);
765 swrite_Chi2SumTest (N, res);
766 swrite_Final (gen, Timer);
767 }
768 util_Free (Urne);
769 if (localRes)
770 sres_DeleteChi2 (res);
771 chrono_Delete (Timer);
772 return 0;
773 }
774
775
776 /*-------------------------------------------------------------------------*/
777
svaria_CollisionArgMax(unif01_Gen * gen,sres_Chi2 * res,long N,long n,int r,long k,long m)778 void svaria_CollisionArgMax (unif01_Gen * gen, sres_Chi2 * res,
779 long N, long n, int r, long k, long m)
780 {
781 if (m > 1) {
782 svaria_CollisionArgMax_00 (gen, res, N, n, r, k, m);
783
784 } else if (m == 1) {
785 double ValDelta[] = { -1.0 };
786 smultin_Param *par;
787
788 if (swrite_Basic) {
789 printf (
790 "***********************************************************\n"
791 "Test svaria_CollisionArgMax calling smultin_Multinomial\n\n");
792 }
793 par = smultin_CreateParam (1, ValDelta, smultin_GenerCellMax, -3);
794 if (NULL == res) {
795 smultin_Multinomial (gen, par, NULL, N, n, r, 0, k, TRUE);
796 } else {
797 smultin_Res *resm;
798 resm = smultin_CreateRes (par);
799 smultin_Multinomial (gen, par, resm, N, n, r, 0, k, TRUE);
800 sres_InitChi2 (res, N, -1, "svaria_CollisionArgMax");
801 statcoll_SetDesc (res->sVal1, "CollisionArgMax sVal1");
802 res->sVal1->NObs = resm->Collector[0]->NObs;
803 tables_CopyTabD (resm->Collector[0]->V, res->sVal1->V, 1, N);
804 tables_CopyTabD (resm->sVal2[0], res->sVal2, 0, gofw_NTestTypes - 1);
805 tables_CopyTabD (resm->pVal2[0], res->pVal2, 0, gofw_NTestTypes - 1);
806 smultin_DeleteRes (resm);
807 }
808 smultin_DeleteParam (par);
809 } else {
810 util_Warning (m <= 0, "svaria_CollisionArgMax: m <= 0");
811 }
812 }
813
814
815 /*=========================================================================*/
816
WriteDataSumColl(unif01_Gen * gen,char * TestName,long N,long n,int r,double g)817 static void WriteDataSumColl (unif01_Gen * gen, char *TestName,
818 long N, long n, int r, double g)
819 {
820 swrite_Head (gen, TestName, N, n, r);
821 printf (", g = %g\n\n", g);
822 }
823
824 /*-------------------------------------------------------------------------*/
825
ProbabiliteG(int jmin,int j,double g)826 static double ProbabiliteG (int jmin, int j, double g)
827 /*
828 * Returns the probability that the minimum number of random U(0,1) whose
829 * sum is larger than g is j+1. g cannot be too large because the
830 * calculation here becomes numerically unstable.
831 */
832 {
833 int s;
834 double temp;
835 const double jLR = j;
836 double signe; /* +1 or -1 */
837 double somme;
838
839 signe = 1.0;
840 somme = 0.0;
841 for (s = 0; s <= jmin; s++) {
842 temp = signe * num2_Combination (j + 1, s);
843 temp *= pow (g - s, jLR);
844 somme += temp;
845 signe = -signe;
846 }
847 somme = (jLR + 1.0 - g) * somme / num2_Factorial (j + 1);
848 return somme;
849 }
850
851
852 /*-------------------------------------------------------------------------*/
853
svaria_SumCollector(unif01_Gen * gen,sres_Chi2 * res,long N,long n,int r,double g)854 void svaria_SumCollector (unif01_Gen * gen, sres_Chi2 * res,
855 long N, long n, int r, double g)
856 {
857 const double gmax = 10.0; /* Maximal value of g */
858 const int jmax = 50; /* Maximal number of classes */
859 int j; /* Class index */
860 long Seq;
861 long i;
862 double X;
863 double Y;
864 double Sum;
865 long NbClasses;
866 long *Loc;
867 double V[1];
868 lebool localRes = FALSE;
869 chrono_Chrono *Timer;
870 char *TestName = "svaria_SumCollector test";
871 char chaine[LEN1 + 1] = "";
872 char str[LEN2 + 1];
873
874 Timer = chrono_Create ();
875 if (swrite_Basic)
876 WriteDataSumColl (gen, TestName, N, n, r, g);
877
878 if (g < 1.0 || g > gmax) {
879 util_Error ("svaria_SumCollector: g < 1.0 or g > 10.0");
880 }
881 if (res == NULL) {
882 localRes = TRUE;
883 res = sres_CreateChi2 ();
884 }
885 sres_InitChi2 (res, N, jmax, "svaria_SumCollector");
886 Loc = res->Loc;
887
888 res->jmin = g;
889 res->jmax = jmax;
890 Sum = 0.0;
891 for (j = res->jmin; j < jmax; j++) {
892 res->NbExp[j] = n * ProbabiliteG (res->jmin, j, g);
893 Sum += res->NbExp[j];
894 }
895 res->NbExp[jmax] = util_Max (0.0, n - Sum);
896
897 if (swrite_Classes)
898 gofs_WriteClasses (res->NbExp, Loc, res->jmin, res->jmax, 0);
899 gofs_MergeClasses (res->NbExp, Loc, &res->jmin, &res->jmax, &NbClasses);
900 if (swrite_Classes)
901 gofs_WriteClasses (res->NbExp, Loc, res->jmin, res->jmax, NbClasses);
902
903 strncpy (chaine, "SumCollector sVal1: chi2 with ", (size_t) LEN1);
904 sprintf (str, "%ld", NbClasses - 1);
905 strncat (chaine, str, (size_t) LEN2);
906 strncat (chaine, " degrees of freedom", (size_t) LEN1);
907 statcoll_SetDesc (res->sVal1, chaine);
908 res->degFree = NbClasses - 1;
909 if (res->degFree < 1) {
910 if (localRes)
911 sres_DeleteChi2 (res);
912 return;
913 }
914
915 for (Seq = 1; Seq <= N; Seq++) {
916 for (j = 1; j <= jmax; j++)
917 res->Count[j] = 0;
918
919 for (i = 1; i <= n; i++) {
920 X = 0.0;
921 j = 0;
922 do {
923 X += unif01_StripD (gen, r);
924 ++j;
925 }
926 while (X <= g);
927 if (j > res->jmax)
928 ++res->Count[res->jmax];
929 else
930 ++res->Count[Loc[j - 1]];
931 }
932 if (swrite_Counters)
933 tables_WriteTabL (res->Count, res->jmin, res->jmax, 5, 10,
934 "Observed numbers:");
935 Y = gofs_Chi2 (res->NbExp, res->Count, res->jmin, res->jmax);
936 statcoll_AddObs (res->sVal1, Y);
937 }
938
939 V[0] = NbClasses - 1;
940 gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_ChiSquare, V,
941 res->sVal2, res->pVal2);
942 res->pVal1->NObs = N;
943 sres_GetChi2SumStat (res);
944
945 if (swrite_Collectors)
946 statcoll_Write (res->sVal1, 5, 14, 4, 3);
947
948 if (swrite_Basic) {
949 swrite_AddStrChi (str, LEN2, res->degFree);
950 gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2, str);
951 swrite_Chi2SumTest (N, res);
952 swrite_Final (gen, Timer);
953 }
954 if (localRes)
955 sres_DeleteChi2 (res);
956 chrono_Delete (Timer);
957 }
958
959
960 /*=========================================================================*/
961
InitAppear(int r,int s,int L,long Q,double E[],double KV[])962 static void InitAppear (int r, int s, int L, long Q, double E[], double KV[])
963 {
964 util_Assert (r >= 0, "svaria_AppearanceSpacings: r < 0");
965 util_Assert (s > 0, "svaria_AppearanceSpacings: s <= 0");
966 /* if (L >= s && L % s) {
967 util_Error ("svaria_AppearanceSpacings: L mod s != 0");
968 }*/
969 if (L < s && s % L) {
970 util_Error ("svaria_AppearanceSpacings: s mod L != 0");
971 }
972 util_Warning (Q < 10.0 * num_TwoExp[L],
973 "svaria_AppearanceSpacings: Q < 10 * 2^L");
974
975 /* Theoretical mean E and variance KV for different L given by Maurer */
976 if (L > 16) {
977 /* For L > 16 and near 16, the 6-th decimal of E[L] could be
978 erroneous. For KV [L], the 4-th decimal. */
979 E[L] = L - 8.32746E-1;
980 KV[L] = 3.423715;
981 } else {
982 E [1] = 0.73264948; KV[1] = 0.68977;
983 E [2] = 1.53743829; KV[2] = 1.33774;
984 E [3] = 2.40160681; KV[3] = 1.90133;
985 E [4] = 3.31122472; KV[4] = 2.35774;
986 E [5] = 4.25342659; KV[5] = 2.70455;
987 E [6] = 5.21770525; KV[6] = 2.95403;
988 E [7] = 6.19625065; KV[7] = 3.12539;
989 E [8] = 7.18366555; KV[8] = 3.23866;
990 E [9] = 8.17642476; KV[9] = 3.31120;
991 E [10] = 9.17232431; KV[10] = 3.35646;
992 E [11] = 10.1700323; KV[11] = 3.38409;
993 E [12] = 11.1687649; KV[12] = 3.40065;
994 E [13] = 12.1680703; KV[13] = 3.41043;
995 E [14] = 13.1676926; KV[14] = 3.41614;
996 E [15] = 14.1674884; KV[15] = 3.41943;
997 E [16] = 15.1673788; KV[16] = 3.42130;
998 }
999 }
1000
1001 /*-------------------------------------------------------------------------*/
1002
CalcSigma(int L,long K,double KV[])1003 static double CalcSigma (int L, long K, double KV[])
1004 /*
1005 * Compute the standard deviation for the svaria_AppearanceSpacings test.
1006 */
1007 {
1008 double dCor[AS_DIM + 1]; /* Coron-Naccache factor d */
1009 double eCor[AS_DIM + 1]; /* Coron-Naccache factor e */
1010 double temp;
1011 /* double c; */
1012
1013 #if 0 /* No correction */
1014 return sqrt (KV[L] / K);
1015
1016 #elif 0 /* Maurer's correction c(L, K) */
1017 temp = 3.0 / L * num_Log2 ((double) K);
1018 if (temp >= DBL_MAX_EXP - 1)
1019 temp = 0.0;
1020 else
1021 temp = pow (2.0, -temp);
1022 c = 0.7 - 0.8/L + (4.0 + 32.0/L) * temp / 15.0;
1023 if (L < 3 || L > 16)
1024 c = 1.0;
1025 return c * sqrt (KV[L] / K);
1026
1027 #else /* based on Coron and Naccache exact calculation */
1028 dCor [3] = 0.2732725; eCor [3] = 0.4890883;
1029 dCor [4] = 0.3045101; eCor [4] = 0.4435381;
1030 dCor [5] = 0.3296587; eCor [5] = 0.4137196;
1031 dCor [6] = 0.3489769; eCor [6] = 0.3941338;
1032 dCor [7] = 0.3631815; eCor [7] = 0.3813210;
1033 dCor [8] = 0.3732189; eCor [8] = 0.3730195;
1034 dCor [9] = 0.3800637; eCor [9] = 0.3677118;
1035 dCor [10] = 0.3845867; eCor [10] = 0.3643695;
1036 dCor [11] = 0.3874942; eCor [11] = 0.3622979;
1037 dCor [12] = 0.3893189; eCor [12] = 0.3610336;
1038 dCor [13] = 0.3904405; eCor [13] = 0.3602731;
1039 dCor [14] = 0.3911178; eCor [14] = 0.3598216;
1040 dCor [15] = 0.3915202; eCor [15] = 0.3595571;
1041 dCor [16] = 0.3917561; eCor [16] = 0.3594040;
1042 /* L = infinite */
1043 dCor [0] = 0.3920729; eCor [0] = 0.3592016;
1044
1045 if (L < 3)
1046 return sqrt (KV[L] / K);
1047 if (L > 16)
1048 temp = dCor[0] + eCor [0] * num_TwoExp[L] / K;
1049 else
1050 temp = dCor[L] + eCor [L] * num_TwoExp[L] / K;
1051 return sqrt (temp * KV[L] / K);
1052
1053 #endif
1054 }
1055
1056 /*-------------------------------------------------------------------------*/
1057
WriteDataAppear(unif01_Gen * gen,long N,int r,int s,int L,long Q,long K,double n)1058 static void WriteDataAppear (unif01_Gen * gen,
1059 long N, int r, int s, int L, long Q, long K, double n)
1060 {
1061 printf ("***********************************************************\n");
1062 printf ("HOST = ");
1063 if (swrite_Host) {
1064 gdef_WriteHostName ();
1065 printf ("\n");
1066 } else
1067 printf ("\n\n");
1068 unif01_WriteNameGen (gen);
1069 printf ("\n");
1070 if (swrite_ExperimentName && strcmp (swrite_ExperimentName, "")) {
1071 printf ("%s", swrite_ExperimentName);
1072 printf (":\n\n");
1073 }
1074
1075 printf ("svaria_AppearanceSpacings test:\n"
1076 "-----------------------------------------------\n");
1077
1078 printf (" N = %2ld, Q = %1ld, K = %1ld, r = %1d, s = %1d,"
1079 " L = %1d\n\n", N, Q, K, r, s, L);
1080 printf (" Sequences of n = (K + Q)L = %12.0f bits\n", n);
1081 printf (" Q = %4ld initialization blocks\n", Q);
1082 printf (" K = %4ld blocks for the test\n", K);
1083 printf (" the blocks have L = %2d bits\n\n\n", L);
1084 }
1085
1086
1087 /*-------------------------------------------------------------------------*/
1088
svaria_AppearanceSpacings(unif01_Gen * gen,sres_Basic * res,long N,long Q,long K,int r,int s,int L)1089 void svaria_AppearanceSpacings (unif01_Gen * gen, sres_Basic * res,
1090 long N, long Q, long K, int r, int s, int L)
1091 {
1092 double E[AS_DIM + 1]; /* Theoretical mean of the log (Base2) of
1093 the most recent occurrence of a block */
1094 double KV[AS_DIM + 1]; /* K times the theoretical variance of the
1095 same */
1096 long Seq;
1097 long block;
1098 long Nblocks; /* 2^L = total number of distinct blocks */
1099 long K2;
1100 long Q2;
1101 long i;
1102 long sBits; /* Numerical value of the s given bits */
1103 long d; /* 2^s */
1104 const int SdivL = s / L;
1105 const int LdivS = L / s;
1106 const int LmodS = L % s;
1107 long sd; /* 2^LmodS */
1108 long rang;
1109 double n; /* Total number of bits in a sequence */
1110 double sigma; /* Standard deviation = sqrt (Variance) */
1111 double somme;
1112 double ARang; /* Most recent occurrence of block */
1113 long *Count; /* Index of most recent occurrence of
1114 block */
1115 double FactMoy;
1116 lebool localRes = FALSE;
1117 chrono_Chrono *Timer;
1118
1119 Timer = chrono_Create ();
1120 n = ((double) K + (double) Q) * L;
1121 if (swrite_Basic)
1122 WriteDataAppear (gen, N, r, s, L, Q, K, n);
1123 util_Assert (s < 32, "svaria_AppearanceSpacings: s >= 32");
1124 InitAppear (r, s, L, Q, E, KV);
1125 sigma = CalcSigma (L, K, KV);
1126 d = num_TwoExp[s];
1127 Nblocks = num_TwoExp[L];
1128 FactMoy = 1.0 / (num_Ln2 * K);
1129
1130 if (res == NULL) {
1131 localRes = TRUE;
1132 res = sres_CreateBasic ();
1133 }
1134 sres_InitBasic (res, N, "svaria_AppearanceSpacings");
1135 Count = util_Calloc ((size_t) Nblocks + 2, sizeof (long));
1136
1137 statcoll_SetDesc (res->sVal1,
1138 "AppearanceSpacings sVal1: standard normal");
1139
1140 if (LdivS > 0) {
1141 sd = num_TwoExp[LmodS];
1142
1143 for (Seq = 1; Seq <= N; Seq++) {
1144 for (i = 0; i < Nblocks; i++)
1145 Count[i] = 0;
1146
1147 /* Initialization with Q blocks */
1148 for (rang = 0; rang < Q; rang++) {
1149 block = 0;
1150 for (i = 1; i <= LdivS; i++) {
1151 sBits = unif01_StripB (gen, r, s);
1152 block = block * d + sBits;
1153 }
1154 if (LmodS > 0) {
1155 sBits = unif01_StripB (gen, r, LmodS);
1156 block = block * sd + sBits;
1157 }
1158 Count[block] = rang;
1159 }
1160
1161 /* Test proper with K blocks */
1162 somme = 0.0;
1163 for (rang = Q; rang < Q + K; rang++) {
1164 block = 0;
1165 for (i = 1; i <= LdivS; i++) {
1166 sBits = unif01_StripB (gen, r, s);
1167 block = block * d + sBits;
1168 }
1169 if (LmodS > 0) {
1170 sBits = unif01_StripB (gen, r, LmodS);
1171 block = block * sd + sBits;
1172 }
1173 ARang = rang - Count[block];
1174 somme += log (ARang);
1175 Count[block] = rang;
1176 }
1177 statcoll_AddObs (res->sVal1, (somme * FactMoy - E[L]) / sigma);
1178 }
1179
1180 } else { /* s > L */
1181 Q2 = Q / SdivL;
1182 K2 = K / SdivL;
1183 for (Seq = 1; Seq <= N; Seq++) {
1184 for (i = 0; i < Nblocks; i++)
1185 Count[i] = 0;
1186
1187 /* Initialization: Q blocks */
1188 for (rang = 0; rang < Q2; rang++) {
1189 sBits = unif01_StripB (gen, r, s);
1190 for (i = 0; i < SdivL; i++) {
1191 block = sBits % Nblocks;
1192 Count[block] = SdivL * rang + i;
1193 sBits /= Nblocks;
1194 }
1195 }
1196 /* Test proper with K blocks */
1197 somme = 0.0;
1198 for (rang = Q2; rang < Q2 + K2; rang++) {
1199 sBits = unif01_StripB (gen, r, s);
1200 for (i = 0; i < SdivL; i++) {
1201 block = sBits % Nblocks;
1202 ARang = SdivL * rang + i - Count[block];
1203 somme += log (ARang);
1204 Count[block] = SdivL * rang + i;
1205 sBits /= Nblocks;
1206 }
1207 }
1208 statcoll_AddObs (res->sVal1, (somme * FactMoy - E[L]) / sigma);
1209 }
1210 }
1211
1212 gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_Normal,
1213 (double *) NULL, res->sVal2, res->pVal2);
1214 res->pVal1->NObs = N;
1215 sres_GetNormalSumStat (res);
1216
1217 if (swrite_Collectors)
1218 statcoll_Write (res->sVal1, 5, 12, 4, 3);
1219
1220 if (swrite_Basic) {
1221 gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2,
1222 "Normal statistic :");
1223 swrite_NormalSumTest (N, res);
1224 swrite_Final (gen, Timer);
1225 }
1226 util_Free (Count);
1227 if (localRes)
1228 sres_DeleteBasic (res);
1229 chrono_Delete (Timer);
1230 }
1231