1 /* movstat/test_qqr.c
2  *
3  * Copyright (C) 2018 Patrick Alken
4  *
5  * This program is free software; you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation; either version 3 of the License, or (at
8  * your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful, but
11  * WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13  * General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18  */
19 
20 #include <gsl/gsl_math.h>
21 #include <gsl/gsl_vector.h>
22 #include <gsl/gsl_test.h>
23 #include <gsl/gsl_rng.h>
24 #include <gsl/gsl_movstat.h>
25 
26 /* compute moving QQR by explicitely constructing window and computing QQR */
27 int
slow_movqqr(const gsl_movstat_end_t etype,const double q,const gsl_vector * x,gsl_vector * y,const int H,const int J)28 slow_movqqr(const gsl_movstat_end_t etype, const double q, const gsl_vector * x,
29             gsl_vector * y, const int H, const int J)
30 {
31   const size_t n = x->size;
32   const int K = H + J + 1;
33   double *window = malloc(K * sizeof(double));
34   size_t i;
35 
36   for (i = 0; i < n; ++i)
37     {
38       size_t wsize = gsl_movstat_fill(etype, x, i, H, J, window);
39       double quant1, quant2;
40 
41       gsl_sort(window, 1, wsize);
42       quant1 = gsl_stats_quantile_from_sorted_data(window, 1, wsize, q);
43       quant2 = gsl_stats_quantile_from_sorted_data(window, 1, wsize, 1.0 - q);
44 
45       gsl_vector_set(y, i, quant2 - quant1);
46     }
47 
48   free(window);
49 
50   return GSL_SUCCESS;
51 }
52 
53 static double
func_qqr(const size_t n,double x[],void * params)54 func_qqr(const size_t n, double x[], void * params)
55 {
56   double q = *(double *) params;
57   double quant1, quant2;
58 
59   (void) params;
60 
61   gsl_sort(x, 1, n);
62   quant1 = gsl_stats_quantile_from_sorted_data(x, 1, n, q);
63   quant2 = gsl_stats_quantile_from_sorted_data(x, 1, n, 1.0 - q);
64 
65   return (quant2 - quant1);
66 }
67 
68 static void
test_qqr_proc(const double tol,const double q,const size_t n,const size_t H,const size_t J,const gsl_movstat_end_t etype,gsl_rng * rng_p)69 test_qqr_proc(const double tol, const double q, const size_t n, const size_t H, const size_t J,
70               const gsl_movstat_end_t etype, gsl_rng * rng_p)
71 {
72   gsl_movstat_workspace * w = gsl_movstat_alloc2(H, J);
73   gsl_vector * x = gsl_vector_alloc(n);
74   gsl_vector * y = gsl_vector_alloc(n);
75   gsl_vector * z = gsl_vector_alloc(n);
76   gsl_movstat_function F;
77   char buf[2048];
78 
79   F.function = func_qqr;
80   F.params = (void *) &q;
81 
82   random_vector(x, rng_p);
83 
84   /* y = QQR(x) with slow brute force method */
85   slow_movqqr(etype, q, x, y, H, J);
86 
87   /* y = QQR(x) with fast method */
88   gsl_movstat_qqr(etype, x, q, z, w);
89 
90   /* test y = z */
91   sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u QQR random", n, H, J, etype);
92   compare_vectors(tol, z, y, buf);
93 
94   /* z = QQR(x) in-place */
95   gsl_vector_memcpy(z, x);
96   gsl_movstat_qqr(etype, z, q, z, w);
97 
98   sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u QQR random in-place", n, H, J, etype);
99   compare_vectors(tol, z, y, buf);
100 
101   /* z = QQR(x) with user-defined function */
102   gsl_movstat_apply(etype, &F, x, z, w);
103 
104   sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u QQR user", n, H, J, etype);
105   compare_vectors(tol, z, y, buf);
106 
107   gsl_movstat_free(w);
108   gsl_vector_free(x);
109   gsl_vector_free(y);
110   gsl_vector_free(z);
111 }
112 
113 static void
test_qqr(gsl_rng * rng_p)114 test_qqr(gsl_rng * rng_p)
115 {
116   const double eps = 1.0e-10;
117 
118   test_qqr_proc(eps, 0.1, 100, 0, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
119   test_qqr_proc(eps, 0.1, 1000, 3, 3, GSL_MOVSTAT_END_PADZERO, rng_p);
120   test_qqr_proc(eps, 0.2, 1000, 0, 5, GSL_MOVSTAT_END_PADZERO, rng_p);
121   test_qqr_proc(eps, 0.3, 1000, 5, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
122   test_qqr_proc(eps, 0.25, 2000, 10, 5, GSL_MOVSTAT_END_PADZERO, rng_p);
123   test_qqr_proc(eps, 0.4, 2000, 5, 10, GSL_MOVSTAT_END_PADZERO, rng_p);
124   test_qqr_proc(eps, 0.25, 20, 50, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
125   test_qqr_proc(eps, 0.25, 20, 10, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
126   test_qqr_proc(eps, 0.25, 20, 50, 10, GSL_MOVSTAT_END_PADZERO, rng_p);
127 
128   test_qqr_proc(eps, 0.1, 100, 0, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
129   test_qqr_proc(eps, 0.1, 1000, 3, 3, GSL_MOVSTAT_END_PADVALUE, rng_p);
130   test_qqr_proc(eps, 0.2, 1000, 0, 5, GSL_MOVSTAT_END_PADVALUE, rng_p);
131   test_qqr_proc(eps, 0.3, 1000, 5, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
132   test_qqr_proc(eps, 0.25, 2000, 10, 5, GSL_MOVSTAT_END_PADVALUE, rng_p);
133   test_qqr_proc(eps, 0.4, 2000, 5, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);
134   test_qqr_proc(eps, 0.25, 20, 50, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
135   test_qqr_proc(eps, 0.25, 20, 10, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
136   test_qqr_proc(eps, 0.25, 20, 50, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);
137 
138   test_qqr_proc(eps, 0.1, 100, 0, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
139   test_qqr_proc(eps, 0.1, 1000, 3, 3, GSL_MOVSTAT_END_TRUNCATE, rng_p);
140   test_qqr_proc(eps, 0.2, 1000, 0, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
141   test_qqr_proc(eps, 0.3, 1000, 5, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
142   test_qqr_proc(eps, 0.25, 2000, 10, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
143   test_qqr_proc(eps, 0.4, 2000, 5, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);
144   test_qqr_proc(eps, 0.25, 20, 50, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
145   test_qqr_proc(eps, 0.25, 20, 10, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
146   test_qqr_proc(eps, 0.25, 20, 50, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);
147 }
148