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 #include <math.h>
  27 #include <stdlib.h>
  28 #include "xalloc.h"
  29 #include <libpspp/compiler.h>
  30
  31 #define _(msgid) gettext (msgid)
  32
  33 static bool
  34 read_values (struct lexer *lexer, double **values, double **weights, size_t *cnt)
  35 {
  36   size_t cap = 0;
  37
  38   *values = NULL;
  39   *weights = NULL;
  40   *cnt = 0;
  41   while (lex_is_number (lexer))
  42     {
  43       double value = lex_tokval (lexer);
  44       double weight = 1.;
  45       lex_get (lexer);
  46       if (lex_match (lexer, '*'))
  47         {
  48           if (!lex_is_number (lexer))
  49             {
  50               lex_error (lexer, _("expecting weight value"));
  51               return false;
  52             }
  53           weight = lex_tokval (lexer);
  54           lex_get (lexer);
  55         }
  56
  57       if (*cnt >= cap)
  58         {
  59           cap = 2 * (cap + 8);
  60           *values = xnrealloc (*values, cap, sizeof **values);
  61           *weights = xnrealloc (*weights, cap, sizeof **weights);
  62         }
  63
  64       (*values)[*cnt] = value;
  65       (*weights)[*cnt] = weight;
  66       (*cnt)++;
  67     }
  68
  69   return true;
  70 }
  71
  72 int
  73 cmd_debug_moments (struct lexer *lexer, struct dataset *ds UNUSED)
  74 {
  75   int retval = CMD_FAILURE;
  76   double *values = NULL;
  77   double *weights = NULL;
  78   double weight, M[4];
  79   int two_pass = 1;
  80   size_t cnt;
  81   size_t i;
  82
  83   if (lex_match_id (lexer, "ONEPASS"))
  84     two_pass = 0;
  85   if (lex_token (lexer) != '/')
  86     {
  87       lex_force_match (lexer, '/');
  88       goto done;
  89     }
  90   fprintf (stderr, "%s => ", lex_rest_of_line (lexer, NULL));
  91   lex_get (lexer);
  92
  93   if (two_pass)
  94     {
  95       struct moments *m = NULL;
  96
  97       m = moments_create (MOMENT_KURTOSIS);
  98       if (!read_values (lexer, &values, &weights, &cnt))
  99         {
 100           moments_destroy (m);
 101           goto done;
 102         }
 103       for (i = 0; i < cnt; i++)
 104         moments_pass_one (m, values[i], weights[i]);
 105       for (i = 0; i < cnt; i++)
 106         moments_pass_two (m, values[i], weights[i]);
 107       moments_calculate (m, &weight, &M[0], &M[1], &M[2], &M[3]);
 108       moments_destroy (m);
 109     }
 110   else
 111     {
 112       struct moments1 *m = NULL;
 113
 114       m = moments1_create (MOMENT_KURTOSIS);
 115       if (!read_values (lexer, &values, &weights, &cnt))
 116         {
 117           moments1_destroy (m);
 118           goto done;
 119         }
 120       for (i = 0; i < cnt; i++)
 121         moments1_add (m, values[i], weights[i]);
 122       moments1_calculate (m, &weight, &M[0], &M[1], &M[2], &M[3]);
 123       moments1_destroy (m);
 124     }
 125
 126   fprintf (stderr, "W=%.3f", weight);
 127   for (i = 0; i < 4; i++)
 128     {
 129       fprintf (stderr, " M%d=", i + 1);
 130       if (M[i] == SYSMIS)
 131         fprintf (stderr, "sysmis");
 132       else if (fabs (M[i]) <= 0.0005)
 133         fprintf (stderr, "0.000");
 134       else
 135         fprintf (stderr, "%.3f", M[i]);
 136     }
 137   fprintf (stderr, "\n");
 138
 139   retval = lex_end_of_command (lexer);
 140
 141  done:
 142   free (values);
 143   free (weights);
 144   return retval;
 145 }