1 /* rstat/rstat.c
2 *
3 * Copyright (C) 2015 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 <config.h>
21 #include <stdlib.h>
22 #include <math.h>
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_rstat.h>
26
27 gsl_rstat_workspace *
gsl_rstat_alloc(void)28 gsl_rstat_alloc(void)
29 {
30 gsl_rstat_workspace *w;
31
32 w = calloc(1, sizeof(gsl_rstat_workspace));
33
34 if (w == 0)
35 {
36 GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
37 }
38
39 w->median_workspace_p = gsl_rstat_quantile_alloc(0.5);
40
41 if (w->median_workspace_p == 0)
42 {
43 GSL_ERROR_NULL ("failed to allocate space for median workspace",
44 GSL_ENOMEM);
45 }
46
47 gsl_rstat_reset(w);
48
49 return w;
50 } /* gsl_rstat_alloc() */
51
52 void
gsl_rstat_free(gsl_rstat_workspace * w)53 gsl_rstat_free(gsl_rstat_workspace *w)
54 {
55 if (w->median_workspace_p)
56 gsl_rstat_quantile_free(w->median_workspace_p);
57
58 free(w);
59 } /* gsl_rstat_free() */
60
61 size_t
gsl_rstat_n(const gsl_rstat_workspace * w)62 gsl_rstat_n(const gsl_rstat_workspace *w)
63 {
64 return w->n;
65 } /* gsl_rstat_n() */
66
67 /* add a data point to the running totals */
68 int
gsl_rstat_add(const double x,gsl_rstat_workspace * w)69 gsl_rstat_add(const double x, gsl_rstat_workspace *w)
70 {
71 double delta = x - w->mean;
72 double delta_n, delta_nsq, term1, n;
73
74 /* update min and max */
75 if (w->n == 0)
76 {
77 w->min = x;
78 w->max = x;
79 }
80 else
81 {
82 if (x < w->min)
83 w->min = x;
84 if (x > w->max)
85 w->max = x;
86 }
87
88 /* update mean and variance */
89 n = (double) ++(w->n);
90 delta_n = delta / n;
91 delta_nsq = delta_n * delta_n;
92 term1 = delta * delta_n * (n - 1.0);
93 w->mean += delta_n;
94 w->M4 += term1 * delta_nsq * (n * n - 3.0 * n + 3.0) +
95 6.0 * delta_nsq * w->M2 - 4.0 * delta_n * w->M3;
96 w->M3 += term1 * delta_n * (n - 2.0) - 3.0 * delta_n * w->M2;
97 w->M2 += term1;
98
99 /* update median */
100 gsl_rstat_quantile_add(x, w->median_workspace_p);
101
102 return GSL_SUCCESS;
103 } /* gsl_rstat_add() */
104
105 double
gsl_rstat_min(const gsl_rstat_workspace * w)106 gsl_rstat_min(const gsl_rstat_workspace *w)
107 {
108 return w->min;
109 } /* gsl_rstat_min() */
110
111 double
gsl_rstat_max(const gsl_rstat_workspace * w)112 gsl_rstat_max(const gsl_rstat_workspace *w)
113 {
114 return w->max;
115 } /* gsl_rstat_max() */
116
117 double
gsl_rstat_mean(const gsl_rstat_workspace * w)118 gsl_rstat_mean(const gsl_rstat_workspace *w)
119 {
120 return w->mean;
121 } /* gsl_rstat_mean() */
122
123 double
gsl_rstat_variance(const gsl_rstat_workspace * w)124 gsl_rstat_variance(const gsl_rstat_workspace *w)
125 {
126 if (w->n > 1)
127 {
128 double n = (double) w->n;
129 return (w->M2 / (n - 1.0));
130 }
131 else
132 return 0.0;
133 } /* gsl_rstat_variance() */
134
135 double
gsl_rstat_sd(const gsl_rstat_workspace * w)136 gsl_rstat_sd(const gsl_rstat_workspace *w)
137 {
138 double var = gsl_rstat_variance(w);
139
140 return (sqrt(var));
141 } /* gsl_rstat_sd() */
142
143 double
gsl_rstat_rms(const gsl_rstat_workspace * w)144 gsl_rstat_rms(const gsl_rstat_workspace *w)
145 {
146 double rms = 0.0;
147
148 if (w->n > 0)
149 {
150 double mean = gsl_rstat_mean(w);
151 double sigma = gsl_rstat_sd(w);
152 double n = (double) w->n;
153 double a = sqrt((n - 1.0) / n);
154 rms = gsl_hypot(mean, a * sigma);
155 }
156
157 return rms;
158 }
159
160 /* standard deviation of the mean: sigma / sqrt(n) */
161 double
gsl_rstat_sd_mean(const gsl_rstat_workspace * w)162 gsl_rstat_sd_mean(const gsl_rstat_workspace *w)
163 {
164 if (w->n > 0)
165 {
166 double sd = gsl_rstat_sd(w);
167 return (sd / sqrt((double) w->n));
168 }
169 else
170 return 0.0;
171 } /* gsl_rstat_sd_mean() */
172
173 double
gsl_rstat_median(gsl_rstat_workspace * w)174 gsl_rstat_median(gsl_rstat_workspace *w)
175 {
176 return gsl_rstat_quantile_get(w->median_workspace_p);
177 }
178
179 double
gsl_rstat_skew(const gsl_rstat_workspace * w)180 gsl_rstat_skew(const gsl_rstat_workspace *w)
181 {
182 if (w->n > 0)
183 {
184 double n = (double) w->n;
185 double fac = pow(n - 1.0, 1.5) / n;
186 return ((fac * w->M3) / pow(w->M2, 1.5));
187 }
188 else
189 return 0.0;
190 } /* gsl_rstat_skew() */
191
192 double
gsl_rstat_kurtosis(const gsl_rstat_workspace * w)193 gsl_rstat_kurtosis(const gsl_rstat_workspace *w)
194 {
195 if (w->n > 0)
196 {
197 double n = (double) w->n;
198 double fac = ((n - 1.0) / n) * (n - 1.0);
199 return ((fac * w->M4) / (w->M2 * w->M2) - 3.0);
200 }
201 else
202 return 0.0;
203 } /* gsl_rstat_kurtosis() */
204
205 int
gsl_rstat_reset(gsl_rstat_workspace * w)206 gsl_rstat_reset(gsl_rstat_workspace *w)
207 {
208 int status;
209
210 w->min = 0.0;
211 w->max = 0.0;
212 w->mean = 0.0;
213 w->M2 = 0.0;
214 w->M3 = 0.0;
215 w->M4 = 0.0;
216 w->n = 0;
217
218 status = gsl_rstat_quantile_reset(w->median_workspace_p);
219
220 return status;
221 } /* gsl_rstat_reset() */
222