1 /*
 2 
 3 App : MTT.cpp
 4 
 5 This is a C++ program for
 6 detecting abrupt climatic change
 7 by using moving t-test technique.
 8 
 9 BY : ZXD
10 Date : January 9, 2014
11 
12 */
13 
14 #define _CRT_SECURE_NO_WARNINGS
15 
16 #include <fstream>
17 #include <cmath>
18 
19 #define N_SIZE 51
20 // The N_SIZE is size of sample,
21 // e.i. the number of total units.
22 
23 #define SUB_SIZE 10
24 // The SUB_SIZE is length of sub-series
25 // e.i. n1 = n2 = SUB_SIZE
26 
27 void smtt( int n, int sub, double *x, double *t );
28 
29 int main()
30 {
31     double x[N_SIZE] = { 0 };
32     double t[N_SIZE] = { 0 };
33     int n = N_SIZE;
34     int sub = SUB_SIZE;
35     int tmp;
36 
37     FILE *infp;
38     FILE *outfp;
39 
40     infp = fopen( "data4.txt", "r" );
41     for ( int i = 0; i < n; i++ )
42     {
43         fscanf( infp, "%d %lf", &tmp, &x[i] );
44         //printf( "%10.4lf\n", x[i] );
45     }
46     fclose( infp );
47 
48     smtt( n, sub, x, t );
49 
50 
51     outfp = fopen( "result.txt", "w" );
52     for ( int i = 0; i < n - sub - sub + 1; i++ )
53     {
54         fprintf( outfp, "%10.4lf\n", t[i] );
55         printf( "%10.4lf\n", t[i] );
56     }
57 
58     fclose( outfp );
59 
60     system( "pause" );
61 }
62 
63 void smtt( int n, int sub, double *x, double *t )
64 {
65     double x1Ave[N_SIZE] = { 0 };
66     double x2Ave[N_SIZE] = { 0 };
67     double x1Var[N_SIZE] = { 0 };
68     double x2Var[N_SIZE] = { 0 };
69     double S[N_SIZE] = { 0 };
70 
71     for ( int i = 0; i < n - sub - sub + 1; i++ )
72     {
73         for ( int j = i; j < i + sub; j++ )
74         {
75             x1Ave[i] += x[j] / (double)sub;
76             x2Ave[i] += x[j + sub] / (double)sub;
77         }
78 
79         for ( int j = i; j < i + sub; j++ )
80         {
81             x1Var[i] += ( x[j] - x1Ave[i] ) * ( x[j] - x1Ave[i] );
82             x2Var[i] += ( x[j + sub] - x1Ave[i] ) * ( x[j + sub] - x1Ave[i] );
83         }
84 
85         x1Var[i] /= (double)sub;
86         x2Var[i] /= (double)sub;
87 
88         S[i] = sqrt( ( (double)sub * x1Var[i] + (double)sub * x1Var[i] ) / (double)( sub + sub - 2 ) );
89 
90         t[i] = ( x1Ave[i] - x2Ave[i] ) / ( S[i] * sqrt( (double)2 / (double)sub ) );
91     }
92 
93 }