+++ /dev/null
-/* Pspp - a program for statistical analysis.
- Copyright (C) 2012 Free Software Foundation, Inc.
-
- This program is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-
-#include <config.h>
-
-#include "jonckheere-terpstra.h"
-
-#include <gsl/gsl_cdf.h>
-#include <math.h>
-
-#include "data/casegrouper.h"
-#include "data/casereader.h"
-#include "data/casewriter.h"
-#include "data/dataset.h"
-#include "data/dictionary.h"
-#include "data/format.h"
-#include "data/subcase.h"
-#include "data/variable.h"
-#include "libpspp/assertion.h"
-#include "libpspp/hmap.h"
-#include "libpspp/message.h"
-#include "libpspp/misc.h"
-#include "math/sort.h"
-#include "output/pivot-table.h"
-
-#include "gl/minmax.h"
-#include "gl/xalloc.h"
-
-#include "gettext.h"
-#define N_(msgid) msgid
-#define _(msgid) gettext (msgid)
-
-
-/* Returns true iff the independent variable lies in the
- between val1 and val2. Regardless of which is the greater value.
-*/
-static bool
-include_func_bi (const struct ccase *c, void *aux)
-{
- const struct n_sample_test *nst = aux;
- const union value *bigger = NULL;
- const union value *smaller = NULL;
-
- if (0 > value_compare_3way (&nst->val1, &nst->val2, var_get_width (nst->indep_var)))
- {
- bigger = &nst->val2;
- smaller = &nst->val1;
- }
- else
- {
- smaller = &nst->val2;
- bigger = &nst->val1;
- }
-
- if (0 < value_compare_3way (smaller, case_data (c, nst->indep_var), var_get_width (nst->indep_var)))
- return false;
-
- if (0 > value_compare_3way (bigger, case_data (c, nst->indep_var), var_get_width (nst->indep_var)))
- return false;
-
- return true;
-}
-
-struct group_data
-{
- /* The total of the caseweights in the group */
- double cc;
-
- /* A casereader containing the group data.
- This casereader contains just two values:
- 0: The raw value of the data
- 1: The cumulative caseweight
- */
- struct casereader *reader;
-};
-
-
-static double
-u (const struct group_data *grp0, const struct group_data *grp1)
-{
- struct ccase *c0;
-
- struct casereader *r0 = casereader_clone (grp0->reader);
- double usum = 0;
- double prev_cc0 = 0.0;
- for (; (c0 = casereader_read (r0)); case_unref (c0))
- {
- struct ccase *c1;
- struct casereader *r1 = casereader_clone (grp1->reader);
- double x0 = case_num_idx (c0, 0);
- double cc0 = case_num_idx (c0, 1);
- double w0 = cc0 - prev_cc0;
-
- double prev_cc1 = 0;
-
- for (; (c1 = casereader_read (r1)); case_unref (c1))
- {
- double x1 = case_num_idx (c1, 0);
- double cc1 = case_num_idx (c1, 1);
-
- if (x0 > x1)
- {
- /* Do nothing */
- }
- else if (x0 < x1)
- {
- usum += w0 * (grp1->cc - prev_cc1);
- case_unref (c1);
- break;
- }
- else
- {
-#if 1
- usum += w0 * ((grp1->cc - prev_cc1) / 2.0);
-#else
- usum += w0 * (grp1->cc - (prev_cc1 + cc1) / 2.0);
-#endif
- case_unref (c1);
- break;
- }
-
- prev_cc1 = cc1;
- }
- casereader_destroy (r1);
- prev_cc0 = cc0;
- }
- casereader_destroy (r0);
-
- return usum;
-}
-
-
-typedef double func_f (double e_l);
-
-/*
- These 3 functions are used repeatedly in the calculation of the
- variance of the JT statistic.
- Having them explicitly defined makes the variance calculation
- a lot simpler.
-*/
-static double
-ff1 (double e)
-{
- return e * (e - 1) * (2*e + 5);
-}
-
-static double
-ff2 (double e)
-{
- return e * (e - 1) * (e - 2);
-}
-
-static double
-ff3 (double e)
-{
- return e * (e - 1) ;
-}
-
-static func_f *mff[3] =
- {
- ff1, ff2, ff3
- };
-
-
-/*
- This function does the following:
- It creates an ordered set of *distinct* values from IR.
- For each case in that set, it calls f[0..N] passing it the caseweight.
- It returns the sum of f[j] in result[j].
-
- result and f must be allocated prior to calling this function.
- */
-static
-void variance_calculation (struct casereader *ir, const struct variable *var,
- const struct dictionary *dict,
- func_f **f, double *result, size_t n)
-{
- int i;
- struct casereader *r = casereader_clone (ir);
- struct ccase *c;
- const struct variable *wv = dict_get_weight (dict);
- const int w_idx = wv ?
- var_get_case_index (wv) :
- caseproto_get_n_widths (casereader_get_proto (r)) ;
-
- r = sort_execute_1var (r, var);
-
- r = casereader_create_distinct (r, var, dict_get_weight (dict));
-
- for (; (c = casereader_read (r)); case_unref (c))
- {
- double w = case_num_idx (c, w_idx);
-
- for (i = 0; i < n; ++i)
- result[i] += f[i] (w);
- }
-
- casereader_destroy (r);
-}
-
-struct jt
-{
- int levels;
- double n;
- double obs;
- double mean;
- double stddev;
-};
-
-static void show_jt (const struct n_sample_test *, const struct jt *,
- const struct fmt_spec *wfmt);
-
-
-void
-jonckheere_terpstra_execute (const struct dataset *ds,
- struct casereader *input,
- enum mv_class exclude,
- const struct npar_test *test,
- bool exact UNUSED,
- double timer UNUSED)
-{
- int v;
- bool warn = true;
- const struct dictionary *dict = dataset_dict (ds);
- const struct n_sample_test *nst = UP_CAST (test, const struct n_sample_test, parent);
-
- struct caseproto *proto = caseproto_create ();
- proto = caseproto_add_width (proto, 0);
- proto = caseproto_add_width (proto, 0);
-
- /* If the independent variable is missing, then we ignore the case */
- input = casereader_create_filter_missing (input,
- &nst->indep_var, 1,
- exclude,
- NULL, NULL);
-
- /* Remove cases with invalid weigths */
- input = casereader_create_filter_weight (input, dict, &warn, NULL);
-
- /* Remove all those cases which are outside the range (val1, val2) */
- input = casereader_create_filter_func (input, include_func_bi, NULL,
- CONST_CAST (struct n_sample_test *, nst), NULL);
-
- /* Sort the data by the independent variable */
- input = sort_execute_1var (input, nst->indep_var);
-
- for (v = 0; v < nst->n_vars ; ++v)
- {
- struct jt jt;
- double variance;
- int g0;
- double nn = 0;
- int i;
- double sums[3] = {0,0,0};
- double e_sum[3] = {0,0,0};
-
- struct group_data *grp = NULL;
- double ccsq_sum = 0;
-
- struct casegrouper *grouper;
- struct casereader *group;
- struct casereader *vreader= casereader_clone (input);
-
- /* Get a few values into e_sum - we'll be needing these later */
- variance_calculation (vreader, nst->vars[v], dict, mff, e_sum, 3);
-
- grouper =
- casegrouper_create_vars (vreader, &nst->indep_var, 1);
-
- jt.obs = 0;
- jt.levels = 0;
- jt.n = 0;
- for (; casegrouper_get_next_group (grouper, &group);
- casereader_destroy (group))
- {
- struct casewriter *writer = autopaging_writer_create (proto);
- struct ccase *c;
- double cc = 0;
-
- group = sort_execute_1var (group, nst->vars[v]);
- for (; (c = casereader_read (group)); case_unref (c))
- {
- struct ccase *c_out = case_create (proto);
-
- *case_num_rw_idx (c_out, 0) = case_num (c, nst->vars[v]);
-
- cc += dict_get_case_weight (dict, c, &warn);
- *case_num_rw_idx (c_out, 1) = cc;
- casewriter_write (writer, c_out);
- }
-
- grp = xrealloc (grp, sizeof *grp * (jt.levels + 1));
-
- grp[jt.levels].reader = casewriter_make_reader (writer);
- grp[jt.levels].cc = cc;
-
- jt.levels++;
- jt.n += cc;
- ccsq_sum += pow2 (cc);
- }
-
- casegrouper_destroy (grouper);
-
- for (g0 = 0; g0 < jt.levels; ++g0)
- {
- int g1;
- for (g1 = g0 +1 ; g1 < jt.levels; ++g1)
- {
- double uu = u (&grp[g0], &grp[g1]);
- jt.obs += uu;
- }
- nn += pow2 (grp[g0].cc) * (2 * grp[g0].cc + 3);
-
- for (i = 0; i < 3; ++i)
- sums[i] += mff[i] (grp[g0].cc);
-
- casereader_destroy (grp[g0].reader);
- }
-
- free (grp);
-
- variance = (mff[0](jt.n) - sums[0] - e_sum[0]) / 72.0;
- variance += sums[1] * e_sum[1] / (36.0 * mff[1] (jt.n));
- variance += sums[2] * e_sum[2] / (8.0 * mff[2] (jt.n));
-
- jt.stddev = sqrt (variance);
-
- jt.mean = (pow2 (jt.n) - ccsq_sum) / 4.0;
-
- show_jt (nst, &jt, dict_get_weight_format (dict));
- }
-
- casereader_destroy (input);
- caseproto_unref (proto);
-}
-\f
-static void
-show_jt (const struct n_sample_test *nst, const struct jt *jt,
- const struct fmt_spec *wfmt)
-{
- struct pivot_table *table = pivot_table_create (
- N_("Jonckheere-Terpstra Test"));
- pivot_table_set_weight_format (table, wfmt);
-
- struct pivot_dimension *statistics = pivot_dimension_create (
- table, PIVOT_AXIS_COLUMN, N_("Statistics"));
- pivot_category_create_leaf_rc (
- statistics->root,
- pivot_value_new_text_format (N_("Number of levels in %s"),
- var_to_string (nst->indep_var)),
- PIVOT_RC_INTEGER);
- pivot_category_create_leaves (
- statistics->root,
- N_("N"), PIVOT_RC_COUNT,
- N_("Observed J-T Statistic"), PIVOT_RC_OTHER,
- N_("Mean J-T Statistic"), PIVOT_RC_OTHER,
- N_("Std. Deviation of J-T Statistic"), PIVOT_RC_OTHER,
- N_("Std. J-T Statistic"), PIVOT_RC_OTHER,
- N_("Asymp. Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE);
-
- struct pivot_dimension *variables = pivot_dimension_create (
- table, PIVOT_AXIS_ROW, N_("Variable"));
-
- for (size_t i = 0; i < nst->n_vars; ++i)
- {
- int row = pivot_category_create_leaf (
- variables->root, pivot_value_new_variable (nst->vars[i]));
-
- double std_jt = (jt[0].obs - jt[0].mean) / jt[0].stddev;
- double sig = (2.0 * (std_jt > 0
- ? gsl_cdf_ugaussian_Q (std_jt)
- : gsl_cdf_ugaussian_P (std_jt)));
- double entries[] = {
- jt[0].levels,
- jt[0].n,
- jt[0].obs,
- jt[0].mean,
- jt[0].stddev,
- std_jt,
- sig,
- };
- for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
- pivot_table_put2 (table, j, row, pivot_value_new_number (entries[j]));
- }
-
- pivot_table_submit (table);
-}
projects / pspp / blobdiff