1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
|
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "include/externs.h"
#include "include/cephes.h"
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C U M U L A T I V E S U M S T E S T
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
int
CumulativeSums(int n, BitSequence *epsilon)
{
int S, sup, inf, z=1, zrev, k;
double sum1, sum2, p_value;
S = 0;
sup = 0;
inf = 0;
for ( k=0; k<n; k++ ) {
epsilon[k] ? S++ : S--;
if ( S > sup )
sup++;
if ( S < inf )
inf--;
z = (sup > -inf) ? sup : -inf;
zrev = (sup-S > S-inf) ? sup-S : S-inf;
}
// forward
sum1 = 0.0;
for ( k=(-n/z+1)/4; k<=(n/z-1)/4; k++ ) {
sum1 += cephes_normal(((4*k+1)*z)/sqrt(n));
sum1 -= cephes_normal(((4*k-1)*z)/sqrt(n));
}
sum2 = 0.0;
for ( k=(-n/z-3)/4; k<=(n/z-1)/4; k++ ) {
sum2 += cephes_normal(((4*k+3)*z)/sqrt(n));
sum2 -= cephes_normal(((4*k+1)*z)/sqrt(n));
}
p_value = 1.0 - sum1 + sum2;
// fprintf(stats[TEST_CUSUM], "\t\t CUMULATIVE SUMS (FORWARD) TEST\n");
// fprintf(stats[TEST_CUSUM], "\t\t-------------------------------------------\n");
// fprintf(stats[TEST_CUSUM], "\t\tCOMPUTATIONAL INFORMATION:\n");
// fprintf(stats[TEST_CUSUM], "\t\t-------------------------------------------\n");
// fprintf(stats[TEST_CUSUM], "\t\t(a) The maximum partial sum = %d\n", z);
// fprintf(stats[TEST_CUSUM], "\t\t-------------------------------------------\n");
if ( isNegative(p_value) || isGreaterThanOne(p_value) ) {
// fprintf(stats[TEST_CUSUM], "\t\tWARNING: P_VALUE IS OUT OF RANGE\n");
return 0;
}
// fprintf(stats[TEST_CUSUM], "%s\t\tp_value = %f\n\n", p_value < ALPHA ? "FAILURE" : "SUCCESS", p_value);
// fprintf(results[TEST_CUSUM], "%f\n", p_value);
if (p_value < ALPHA) {
return 0;
}
// backwards
sum1 = 0.0;
for ( k=(-n/zrev+1)/4; k<=(n/zrev-1)/4; k++ ) {
sum1 += cephes_normal(((4*k+1)*zrev)/sqrt(n));
sum1 -= cephes_normal(((4*k-1)*zrev)/sqrt(n));
}
sum2 = 0.0;
for ( k=(-n/zrev-3)/4; k<=(n/zrev-1)/4; k++ ) {
sum2 += cephes_normal(((4*k+3)*zrev)/sqrt(n));
sum2 -= cephes_normal(((4*k+1)*zrev)/sqrt(n));
}
p_value = 1.0 - sum1 + sum2;
// fprintf(stats[TEST_CUSUM], "\t\t CUMULATIVE SUMS (REVERSE) TEST\n");
// fprintf(stats[TEST_CUSUM], "\t\t-------------------------------------------\n");
// fprintf(stats[TEST_CUSUM], "\t\tCOMPUTATIONAL INFORMATION:\n");
// fprintf(stats[TEST_CUSUM], "\t\t-------------------------------------------\n");
// fprintf(stats[TEST_CUSUM], "\t\t(a) The maximum partial sum = %d\n", zrev);
// fprintf(stats[TEST_CUSUM], "\t\t-------------------------------------------\n");
if ( isNegative(p_value) || isGreaterThanOne(p_value) ) {
// fprintf(stats[TEST_CUSUM], "\t\tWARNING: P_VALUE IS OUT OF RANGE\n");
return 0;
}
// fprintf(stats[TEST_CUSUM], "%s\t\tp_value = %f\n\n", p_value < ALPHA ? "FAILURE" : "SUCCESS", p_value); fflush(stats[TEST_CUSUM]);
// fprintf(results[TEST_CUSUM], "%f\n", p_value); fflush(results[TEST_CUSUM]);
if (p_value < ALPHA) {
return 0;
} else {
return 1;
}
}
|
