src/language/tests/moments-test.c

   1 /* PSPP - computes sample statistics.
   2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
   3    Written by Ben Pfaff <blp@gnu.org>.
   4
   5    This program is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU General Public License as
   7    published by the Free Software Foundation; either version 2 of the
   8    License, or (at 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
  18    02110-1301, USA. */
  19
  20 #include <config.h>
  21 #include <stdio.h>
  22 #include "gettext.h"
  23 #include <language/command.h>
  24 #include <language/lexer/lexer.h>
  25 #include <math/moments.h>
  26
  27 #define _(msgid) gettext (msgid)
  28
  29 static int
  30 read_values (double **values, double **weights, size_t *cnt)
  31 {
  32   size_t cap = 0;
  33
  34   *values = NULL;
  35   *weights = NULL;
  36   *cnt = 0;
  37   while (lex_is_number ())
  38     {
  39       double value = tokval;
  40       double weight = 1.;
  41       lex_get ();
  42       if (lex_match ('*'))
  43         {
  44           if (!lex_is_number ())
  45             {
  46               lex_error (_("expecting weight value"));
  47               return 0;
  48             }
  49           weight = tokval;
  50           lex_get ();
  51         }
  52
  53       if (*cnt >= cap)
  54         {
  55           cap = 2 * (cap + 8);
  56           *values = xnrealloc (*values, cap, sizeof **values);
  57           *weights = xnrealloc (*weights, cap, sizeof **weights);
  58         }
  59
  60       (*values)[*cnt] = value;
  61       (*weights)[*cnt] = weight;
  62       (*cnt)++;
  63     }
  64
  65   return 1;
  66 }
  67
  68 int
  69 cmd_debug_moments (void)
  70 {
  71   int retval = CMD_FAILURE;
  72   double *values = NULL;
  73   double *weights = NULL;
  74   double weight, M[4];
  75   int two_pass = 1;
  76   size_t cnt;
  77   size_t i;
  78
  79   if (lex_match_id ("ONEPASS"))
  80     two_pass = 0;
  81   if (token != '/')
  82     {
  83       lex_force_match ('/');
  84       goto done;
  85     }
  86   fprintf (stderr, "%s => ", lex_rest_of_line (NULL));
  87   lex_get ();
  88
  89   if (two_pass)
  90     {
  91       struct moments *m = NULL;
  92
  93       m = moments_create (MOMENT_KURTOSIS);
  94       if (!read_values (&values, &weights, &cnt))
  95         {
  96           moments_destroy (m);
  97           goto done;
  98         }
  99       for (i = 0; i < cnt; i++)
 100         moments_pass_one (m, values[i], weights[i]);
 101       for (i = 0; i < cnt; i++)
 102         moments_pass_two (m, values[i], weights[i]);
 103       moments_calculate (m, &weight, &M[0], &M[1], &M[2], &M[3]);
 104       moments_destroy (m);
 105     }
 106   else
 107     {
 108       struct moments1 *m = NULL;
 109
 110       m = moments1_create (MOMENT_KURTOSIS);
 111       if (!read_values (&values, &weights, &cnt))
 112         {
 113           moments1_destroy (m);
 114           goto done;
 115         }
 116       for (i = 0; i < cnt; i++)
 117         moments1_add (m, values[i], weights[i]);
 118       moments1_calculate (m, &weight, &M[0], &M[1], &M[2], &M[3]);
 119       moments1_destroy (m);
 120     }
 121
 122   fprintf (stderr, "W=%.3f", weight);
 123   for (i = 0; i < 4; i++)
 124     {
 125       fprintf (stderr, " M%d=", i + 1);
 126       if (M[i] == SYSMIS)
 127         fprintf (stderr, "sysmis");
 128       else if (fabs (M[i]) <= 0.0005)
 129         fprintf (stderr, "0.000");
 130       else
 131         fprintf (stderr, "%.3f", M[i]);
 132     }
 133   fprintf (stderr, "\n");
 134
 135   retval = lex_end_of_command ();
 136
 137  done:
 138   free (values);
 139   free (weights);
 140   return retval;
 141 }