1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2011 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
19 #include "language/stats/ks-one-sample.h"
21 #include <gsl/gsl_cdf.h>
26 #include "math/sort.h"
27 #include "data/case.h"
28 #include "data/casereader.h"
29 #include "data/dataset.h"
30 #include "data/dictionary.h"
31 #include "data/format.h"
32 #include "data/value-labels.h"
33 #include "data/variable.h"
34 #include "language/stats/freq.h"
35 #include "language/stats/npar.h"
36 #include "libpspp/array.h"
37 #include "libpspp/assertion.h"
38 #include "libpspp/cast.h"
39 #include "libpspp/compiler.h"
40 #include "libpspp/hash-functions.h"
41 #include "libpspp/message.h"
42 #include "libpspp/misc.h"
43 #include "output/pivot-table.h"
45 #include "gl/xalloc.h"
48 #define N_(msgid) msgid
49 #define _(msgid) gettext (msgid)
52 /* The per test variable statistics */
69 typedef double theoretical (const struct ks *ks, double x);
70 typedef theoretical *theoreticalfp;
73 theoretical_uniform (const struct ks *ks, double x)
75 return gsl_cdf_flat_P (x, ks->test_min, ks->test_max);
79 theoretical_normal (const struct ks *ks, double x)
81 return gsl_cdf_gaussian_P (x - ks->mu, ks->sigma);
85 theoretical_poisson (const struct ks *ks, double x)
87 return gsl_cdf_poisson_P (x, ks->mu);
91 theoretical_exponential (const struct ks *ks, double x)
93 return gsl_cdf_exponential_P (x, 1/ks->mu);
97 static const theoreticalfp theoreticalf[4] =
102 theoretical_exponential
106 Return the assymptotic approximation to the significance of Z
109 ks_asymp_sig (double z)
119 double q = exp (-1.233701 * pow (z, -2));
120 return 1 - 2.506628 * (q + pow (q, 9) + pow (q, 25))/ z ;
124 double q = exp (-2 * z * z);
125 return 2 * (q - pow (q, 4) + pow (q, 9) - pow (q, 16))/ z ;
129 static void show_results (const struct ks *, const struct ks_one_sample_test *, const struct fmt_spec *);
133 ks_one_sample_execute (const struct dataset *ds,
134 struct casereader *input,
135 enum mv_class exclude,
136 const struct npar_test *test,
137 bool x UNUSED, double y UNUSED)
139 const struct dictionary *dict = dataset_dict (ds);
140 const struct ks_one_sample_test *kst = UP_CAST (test, const struct ks_one_sample_test, parent.parent);
141 const struct one_sample_test *ost = &kst->parent;
143 const struct fmt_spec *wfmt = dict_get_weight_format (dict);
146 struct casereader *r = casereader_clone (input);
148 struct ks *ks = XCALLOC (ost->n_vars, struct ks);
150 for (v = 0; v < ost->n_vars; ++v)
153 ks[v].test_min = DBL_MAX;
154 ks[v].test_max = -DBL_MAX;
155 ks[v].diff_pos = -DBL_MAX;
156 ks[v].diff_neg = DBL_MAX;
161 for (; (c = casereader_read (r)) != NULL; case_unref (c))
163 const double weight = dict_get_case_weight (dict, c, &warn);
165 for (v = 0; v < ost->n_vars; ++v)
167 const struct variable *var = ost->vars[v];
168 const union value *val = case_data (c, var);
170 if (var_is_value_missing (var, val, exclude))
173 minimize (&ks[v].test_min, val->f);
174 maximize (&ks[v].test_max, val->f);
176 ks[v].obs_cc += weight;
178 ks[v].ssq += pow2 (val->f);
181 casereader_destroy (r);
183 for (v = 0; v < ost->n_vars; ++v)
185 const struct variable *var = ost->vars[v];
187 double prev_empirical = 0;
192 if (kst->p[0] != SYSMIS)
193 ks[v].test_min = kst->p[0];
195 if (kst->p[1] != SYSMIS)
196 ks[v].test_max = kst->p[1];
199 if (kst->p[0] != SYSMIS)
200 ks[v].mu = kst->p[0];
202 ks[v].mu = ks[v].sum / ks[v].obs_cc;
204 if (kst->p[1] != SYSMIS)
205 ks[v].sigma = kst->p[1];
208 ks[v].sigma = ks[v].ssq - pow2 (ks[v].sum) / ks[v].obs_cc;
209 ks[v].sigma /= ks[v].obs_cc - 1;
210 ks[v].sigma = sqrt (ks[v].sigma);
216 if (kst->p[0] != SYSMIS)
217 ks[v].mu = ks[v].sigma = kst->p[0];
219 ks[v].mu = ks[v].sigma = ks[v].sum / ks[v].obs_cc;
225 r = sort_execute_1var (casereader_clone (input), var);
226 for (; (c = casereader_read (r)) != NULL; case_unref (c))
228 double theoretical, empirical;
230 const double weight = dict_get_case_weight (dict, c, &warn);
231 const union value *val = case_data (c, var);
233 if (var_is_value_missing (var, val, exclude))
238 empirical = cc / ks[v].obs_cc;
240 theoretical = theoreticalf[kst->dist] (&ks[v], val->f);
242 d = empirical - theoretical;
243 dp = prev_empirical - theoretical;
246 maximize (&ks[v].diff_pos, d);
248 minimize (&ks[v].diff_neg, d);
251 maximize (&ks[v].diff_pos, dp);
253 minimize (&ks[v].diff_neg, dp);
255 prev_empirical = empirical;
258 casereader_destroy (r);
261 show_results (ks, kst, wfmt);
264 casereader_destroy (input);
269 show_results (const struct ks *ks,
270 const struct ks_one_sample_test *kst,
271 const struct fmt_spec *wfmt)
273 struct pivot_table *table = pivot_table_create (
274 N_("One-Sample Kolmogorov-Smirnov Test"));
275 pivot_table_set_weight_format (table, wfmt);
277 struct pivot_dimension *statistics = pivot_dimension_create (
278 table, PIVOT_AXIS_ROW, N_("Statistics"),
279 N_("N"), PIVOT_RC_COUNT);
284 pivot_category_create_group (statistics->root, N_("Uniform Parameters"),
285 N_("Minimum"), N_("Maximum"));
289 pivot_category_create_group (statistics->root, N_("Normal Parameters"),
290 N_("Mean"), N_("Std. Deviation"));
294 pivot_category_create_group (statistics->root, N_("Poisson Parameters"),
299 pivot_category_create_group (statistics->root,
300 N_("Exponential Parameters"), N_("Scale"));
307 pivot_category_create_group (
308 statistics->root, N_("Most Extreme Differences"),
309 N_("Absolute"), N_("Positive"), N_("Negative"));
311 pivot_category_create_leaves (
312 statistics->root, N_("Kolmogorov-Smirnov Z"),
313 _("Asymp. Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE);
315 struct pivot_dimension *variables = pivot_dimension_create (
316 table, PIVOT_AXIS_COLUMN, N_("Variables"));
318 for (size_t i = 0; i < kst->parent.n_vars; ++i)
320 int col = pivot_category_create_leaf (
321 variables->root, pivot_value_new_variable (kst->parent.vars[i]));
326 values[n++] = ks[i].obs_cc;
331 values[n++] = ks[i].test_min;
332 values[n++] = ks[i].test_max;
336 values[n++] = ks[i].mu;
337 values[n++] = ks[i].sigma;
342 values[n++] = ks[i].mu;
349 double abs = ks[i].diff_pos;
350 maximize (&abs, -ks[i].diff_neg);
352 double z = sqrt (ks[i].obs_cc) * abs;
355 values[n++] = ks[i].diff_pos;
356 values[n++] = ks[i].diff_neg;
358 values[n++] = ks_asymp_sig (z);
360 for (size_t j = 0; j < n; j++)
361 pivot_table_put2 (table, j, col, pivot_value_new_number (values[j]));
364 pivot_table_submit (table);