pintos-os.org Git - pspp/blob - src/language/stats/t-test-one-sample.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2011 Free Software Foundation, Inc.
   3
   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.
   8
   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.
  13
  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/>.
  16 */
  17
  18 #include <config.h>
  19
  20
  21 #include "t-test.h"
  22
  23 #include <math.h>
  24 #include <gsl/gsl_cdf.h>
  25
  26 #include "data/variable.h"
  27 #include "data/format.h"
  28 #include "data/casereader.h"
  29 #include "data/dictionary.h"
  30 #include "libpspp/hash-functions.h"
  31 #include "libpspp/hmapx.h"
  32 #include "math/moments.h"
  33
  34 #include "output/pivot-table.h"
  35
  36 #include "gettext.h"
  37 #define N_(msgid) msgid
  38 #define _(msgid) gettext (msgid)
  39
  40
  41 struct per_var_stats
  42 {
  43   const struct variable *var;
  44
  45   /* N, Mean, Variance */
  46   struct moments *mom;
  47
  48   /* Sum of the differences */
  49   double sum_diff;
  50 };
  51
  52
  53 struct one_samp
  54 {
  55   struct per_var_stats *stats;
  56   size_t n_stats;
  57   double testval;
  58 };
  59
  60
  61 static void
  62 one_sample_test (const struct tt *tt, const struct one_samp *os)
  63 {
  64   struct pivot_table *table = pivot_table_create (N_("One-Sample Test"));
  65   pivot_table_set_weight_var (table, tt->wv);
  66
  67   struct pivot_dimension *statistics = pivot_dimension_create (
  68     table, PIVOT_AXIS_COLUMN, N_("Statistics"));
  69   struct pivot_category *group = pivot_category_create_group__ (
  70     statistics->root, pivot_value_new_user_text_nocopy (
  71       xasprintf (_("Test Value = %.*g"), DBL_DIG + 1, os->testval)));
  72   pivot_category_create_leaves (
  73     group,
  74     N_("t"), PIVOT_RC_OTHER,
  75     N_("df"), PIVOT_RC_COUNT,
  76     N_("Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE,
  77     N_("Mean Difference"), PIVOT_RC_OTHER);
  78   struct pivot_category *subgroup = pivot_category_create_group__ (
  79     group, pivot_value_new_user_text_nocopy (
  80       xasprintf (_("%g%% Confidence Interval of the Difference"),
  81                  tt->confidence * 100.0)));
  82   pivot_category_create_leaves (subgroup,
  83                                 N_("Lower"), PIVOT_RC_OTHER,
  84                                 N_("Upper"), PIVOT_RC_OTHER);
  85
  86   struct pivot_dimension *dep_vars = pivot_dimension_create (
  87     table, PIVOT_AXIS_ROW, N_("Dependent Variables"));
  88
  89   for (size_t i = 0; i < os->n_stats; i++)
  90     {
  91       const struct per_var_stats *per_var_stats = &os->stats[i];
  92       const struct moments *m = per_var_stats->mom;
  93
  94       int dep_var_idx = pivot_category_create_leaf (
  95         dep_vars->root, pivot_value_new_variable (per_var_stats->var));
  96
  97       double cc, mean, sigma;
  98       moments_calculate (m, &cc, &mean, &sigma, NULL, NULL);
  99       double tval = (mean - os->testval) * sqrt (cc / sigma);
 100       double mean_diff = per_var_stats->sum_diff / cc;
 101       double se_mean = sqrt (sigma / cc);
 102       double df = cc - 1.0;
 103       double p = gsl_cdf_tdist_P (tval, df);
 104       double q = gsl_cdf_tdist_Q (tval, df);
 105       double sig = 2.0 * (tval > 0 ? q : p);
 106       double tval_qinv = gsl_cdf_tdist_Qinv ((1.0 - tt->confidence) / 2.0, df);
 107       double lower = mean_diff - tval_qinv * se_mean;
 108       double upper = mean_diff + tval_qinv * se_mean;
 109
 110       double entries[] = { tval, df, sig, mean_diff, lower, upper };
 111       for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
 112         pivot_table_put2 (table, j, dep_var_idx,
 113                           pivot_value_new_number (entries[j]));
 114     }
 115
 116   pivot_table_submit (table);
 117 }
 118
 119 static void
 120 one_sample_summary (const struct tt *tt, const struct one_samp *os)
 121 {
 122   struct pivot_table *table = pivot_table_create (N_("One-Sample Statistics"));
 123   pivot_table_set_weight_var (table, tt->wv);
 124
 125   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
 126                           N_("N"), PIVOT_RC_COUNT,
 127                           N_("Mean"), PIVOT_RC_OTHER,
 128                           N_("Std. Deviation"), PIVOT_RC_OTHER,
 129                           N_("S.E. Mean"), PIVOT_RC_OTHER);
 130
 131   struct pivot_dimension *variables = pivot_dimension_create (
 132     table, PIVOT_AXIS_ROW, N_("Variables"));
 133
 134   for (size_t i = 0; i < os->n_stats; i++)
 135     {
 136       const struct per_var_stats *per_var_stats = &os->stats[i];
 137       const struct moments *m = per_var_stats->mom;
 138
 139       int var_idx = pivot_category_create_leaf (
 140         variables->root, pivot_value_new_variable (per_var_stats->var));
 141
 142       double cc, mean, sigma;
 143       moments_calculate (m, &cc, &mean, &sigma, NULL, NULL);
 144
 145       double entries[] = { cc, mean, sqrt (sigma), sqrt (sigma / cc) };
 146       for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
 147         pivot_table_put2 (table, j, var_idx,
 148                           pivot_value_new_number (entries[j]));
 149     }
 150
 151   pivot_table_submit (table);
 152 }
 153
 154 void
 155 one_sample_run (const struct tt *tt, double testval, struct casereader *reader)
 156 {
 157   struct one_samp os;
 158   os.testval = testval;
 159   os.stats = xcalloc (tt->n_vars, sizeof *os.stats);
 160   os.n_stats = tt->n_vars;
 161   for (size_t i = 0; i < tt->n_vars; ++i)
 162     {
 163       struct per_var_stats *per_var_stats = &os.stats[i];
 164       per_var_stats->var = tt->vars[i];
 165       per_var_stats->mom = moments_create (MOMENT_VARIANCE);
 166     }
 167
 168   struct casereader *r = casereader_clone (reader);
 169   struct ccase *c;
 170   for (; (c = casereader_read (r)); case_unref (c))
 171     {
 172       double w = dict_get_case_weight (tt->dict, c, NULL);
 173       for (size_t i = 0; i < os.n_stats; i++)
 174         {
 175           const struct per_var_stats *per_var_stats = &os.stats[i];
 176           const struct variable *var = per_var_stats->var;
 177           const union value *val = case_data (c, var);
 178           if (var_is_value_missing (var, val) & tt->exclude)
 179             continue;
 180
 181           moments_pass_one (per_var_stats->mom, val->f, w);
 182         }
 183     }
 184   casereader_destroy (r);
 185
 186   r = reader;
 187   for (; (c = casereader_read (r)); case_unref (c))
 188     {
 189       double w = dict_get_case_weight (tt->dict, c, NULL);
 190       for (size_t i = 0; i < os.n_stats; i++)
 191         {
 192           struct per_var_stats *per_var_stats = &os.stats[i];
 193           const struct variable *var = per_var_stats->var;
 194           const union value *val = case_data (c, var);
 195           if (var_is_value_missing (var, val) & tt->exclude)
 196             continue;
 197
 198           moments_pass_two (per_var_stats->mom, val->f, w);
 199           per_var_stats->sum_diff += w * (val->f - os.testval);
 200         }
 201     }
 202   casereader_destroy (r);
 203
 204   one_sample_summary (tt, &os);
 205   one_sample_test (tt, &os);
 206
 207   for (size_t i = 0; i < os.n_stats; i++)
 208     {
 209       const struct per_var_stats *per_var_stats = &os.stats[i];
 210       moments_destroy (per_var_stats->mom);
 211     }
 212   free (os.stats);
 213 }
 214