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/tab.h"
  35
  36 #include "gettext.h"
  37 #define _(msgid) gettext (msgid)
  38
  39
  40 struct per_var_stats
  41 {
  42   const struct variable *var;
  43
  44   /* The position for reporting purposes */
  45   int posn;
  46
  47   /* N, Mean, Variance */
  48   struct moments *mom;
  49
  50   /* Sum of the differences */
  51   double sum_diff;
  52 };
  53
  54
  55 struct one_samp
  56 {
  57   struct hmapx hmap;
  58   double testval;
  59 };
  60
  61
  62 static void
  63 one_sample_test (const struct tt *tt, const struct one_samp *os)
  64 {
  65   struct hmapx_node *node;
  66   struct per_var_stats *per_var_stats;
  67
  68   const int heading_rows = 3;
  69   const size_t rows = heading_rows + tt->n_vars;
  70   const size_t cols = 7;
  71   const struct fmt_spec *wfmt = tt->wv ? var_get_print_format (tt->wv) : & F_8_0;
  72
  73   struct tab_table *t = tab_create (cols, rows);
  74   tab_set_format (t, RC_WEIGHT, wfmt);
  75
  76   tab_headers (t, 0, 0, heading_rows, 0);
  77   tab_box (t, TAL_2, TAL_2, TAL_0, TAL_0, 0, 0, cols - 1, rows - 1);
  78   tab_hline (t, TAL_2, 0, cols - 1, 3);
  79
  80   tab_title (t, _("One-Sample Test"));
  81   tab_hline (t, TAL_1, 1, cols - 1, 1);
  82   tab_vline (t, TAL_2, 1, 0, rows - 1);
  83
  84   tab_joint_text_format (t, 1, 0, cols - 1, 0, TAB_CENTER,
  85                          _("Test Value = %f"), os->testval);
  86
  87   tab_box (t, -1, -1, -1, TAL_1, 1, 1, cols - 1, rows - 1);
  88
  89   tab_joint_text_format (t, 5, 1, 6, 1, TAB_CENTER,
  90                          _("%g%% Confidence Interval of the Difference"),
  91                          tt->confidence * 100.0);
  92
  93   tab_vline (t, TAL_GAP, 6, 1, 1);
  94   tab_hline (t, TAL_1, 5, 6, 2);
  95   tab_text (t, 1, 2, TAB_CENTER | TAT_TITLE, _("t"));
  96   tab_text (t, 2, 2, TAB_CENTER | TAT_TITLE, _("df"));
  97   tab_text (t, 3, 2, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
  98   tab_text (t, 4, 2, TAB_CENTER | TAT_TITLE, _("Mean Difference"));
  99   tab_text (t, 5, 2, TAB_CENTER | TAT_TITLE, _("Lower"));
 100   tab_text (t, 6, 2, TAB_CENTER | TAT_TITLE, _("Upper"));
 101
 102   HMAPX_FOR_EACH (per_var_stats, node, &os->hmap)
 103     {
 104       const struct moments *m = per_var_stats->mom;
 105       double cc, mean, sigma;
 106       double tval, df;
 107       double p, q;
 108       double mean_diff;
 109       double se_mean ;
 110       const int v = per_var_stats->posn;
 111
 112       moments_calculate (m, &cc, &mean, &sigma, NULL, NULL);
 113
 114       tval = (mean - os->testval) * sqrt (cc / sigma);
 115
 116       mean_diff = per_var_stats->sum_diff / cc;
 117       se_mean = sqrt (sigma / cc);
 118       df = cc - 1.0;
 119       p = gsl_cdf_tdist_P (tval, df);
 120       q = gsl_cdf_tdist_Q (tval, df);
 121
 122       tab_text (t, 0, v + heading_rows, TAB_LEFT, var_to_string (per_var_stats->var));
 123       tab_double (t, 1, v + heading_rows, TAB_RIGHT, tval, NULL, RC_OTHER);
 124       tab_double (t, 2, v + heading_rows, TAB_RIGHT, df, NULL, RC_WEIGHT);
 125
 126       /* Multiply by 2 to get 2-tailed significance, makeing sure we've got
 127          the correct tail*/
 128       tab_double (t, 3, v + heading_rows, TAB_RIGHT, 2.0 * (tval > 0 ? q : p), NULL, RC_PVALUE);
 129
 130       tab_double (t, 4, v + heading_rows, TAB_RIGHT, mean_diff,  NULL, RC_OTHER);
 131
 132       tval = gsl_cdf_tdist_Qinv ( (1.0 - tt->confidence) / 2.0, df);
 133
 134       tab_double (t, 5, v + heading_rows, TAB_RIGHT, mean_diff - tval * se_mean, NULL, RC_OTHER);
 135       tab_double (t, 6, v + heading_rows, TAB_RIGHT, mean_diff + tval * se_mean, NULL, RC_OTHER);
 136     }
 137
 138   tab_submit (t);
 139 }
 140
 141 static void
 142 one_sample_summary (const struct tt *tt, const struct one_samp *os)
 143 {
 144   struct hmapx_node *node;
 145   struct per_var_stats *per_var_stats;
 146
 147   const int cols = 5;
 148   const int heading_rows = 1;
 149   const int rows = tt->n_vars + heading_rows;
 150   struct tab_table *t = tab_create (cols, rows);
 151   const struct fmt_spec *wfmt = tt->wv ? var_get_print_format (tt->wv) : & F_8_0;
 152   tab_set_format (t, RC_WEIGHT, wfmt);
 153   tab_headers (t, 0, 0, heading_rows, 0);
 154   tab_box (t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, cols - 1, rows - 1);
 155   tab_hline (t, TAL_2, 0, cols - 1, 1);
 156
 157   tab_title (t, _("One-Sample Statistics"));
 158   tab_vline (t, TAL_2, 1, 0, rows - 1);
 159   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("N"));
 160   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
 161   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
 162   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("S.E. Mean"));
 163
 164   HMAPX_FOR_EACH (per_var_stats, node, &os->hmap)
 165     {
 166       const struct moments *m = per_var_stats->mom;
 167       const int v = per_var_stats->posn;
 168       double cc, mean, sigma;
 169       moments_calculate (m, &cc, &mean, &sigma, NULL, NULL);
 170
 171       tab_text (t, 0, v + heading_rows, TAB_LEFT, var_to_string (per_var_stats->var));
 172       tab_double (t, 1, v + heading_rows, TAB_RIGHT, cc, NULL, RC_WEIGHT);
 173       tab_double (t, 2, v + heading_rows, TAB_RIGHT, mean, NULL, RC_OTHER);
 174       tab_double (t, 3, v + heading_rows, TAB_RIGHT, sqrt (sigma), NULL, RC_OTHER);
 175       tab_double (t, 4, v + heading_rows, TAB_RIGHT, sqrt (sigma / cc), NULL, RC_OTHER);
 176     }
 177
 178   tab_submit (t);
 179 }
 180
 181 void
 182 one_sample_run (const struct tt *tt, double testval, struct casereader *reader)
 183 {
 184   int i;
 185   struct ccase *c;
 186   struct one_samp os;
 187   struct casereader *r;
 188   struct hmapx_node *node;
 189   struct per_var_stats *per_var_stats;
 190
 191   os.testval = testval;
 192   hmapx_init (&os.hmap);
 193
 194   /* Insert all the variables into the map */
 195   for (i = 0; i < tt->n_vars; ++i)
 196     {
 197       struct per_var_stats *per_var_stats = xzalloc (sizeof (*per_var_stats));
 198
 199       per_var_stats->posn = i;
 200       per_var_stats->var = tt->vars[i];
 201       per_var_stats->mom = moments_create (MOMENT_VARIANCE);
 202
 203       hmapx_insert (&os.hmap, per_var_stats, hash_pointer (per_var_stats->var, 0));
 204     }
 205
 206   r = casereader_clone (reader);
 207   for ( ; (c = casereader_read (r) ); case_unref (c))
 208     {
 209       double w = dict_get_case_weight (tt->dict, c, NULL);
 210       struct hmapx_node *node;
 211       struct per_var_stats *per_var_stats;
 212       HMAPX_FOR_EACH (per_var_stats, node, &os.hmap)
 213         {
 214           const struct variable *var = per_var_stats->var;
 215           const union value *val = case_data (c, var);
 216           if (var_is_value_missing (var, val, tt->exclude))
 217             continue;
 218
 219           moments_pass_one (per_var_stats->mom, val->f, w);
 220         }
 221     }
 222   casereader_destroy (r);
 223
 224   r = reader;
 225   for ( ; (c = casereader_read (r) ); case_unref (c))
 226     {
 227       double w = dict_get_case_weight (tt->dict, c, NULL);
 228       struct hmapx_node *node;
 229       struct per_var_stats *per_var_stats;
 230       HMAPX_FOR_EACH (per_var_stats, node, &os.hmap)
 231         {
 232           const struct variable *var = per_var_stats->var;
 233           const union value *val = case_data (c, var);
 234           if (var_is_value_missing (var, val, tt->exclude))
 235             continue;
 236
 237           moments_pass_two (per_var_stats->mom, val->f, w);
 238           per_var_stats->sum_diff += w * (val->f - os.testval);
 239         }
 240     }
 241   casereader_destroy (r);
 242
 243   one_sample_summary (tt, &os);
 244   one_sample_test (tt, &os);
 245
 246   HMAPX_FOR_EACH (per_var_stats, node, &os.hmap)
 247     {
 248       moments_destroy (per_var_stats->mom);
 249       free (per_var_stats);
 250     }
 251
 252   hmapx_destroy (&os.hmap);
 253 }
 254