src/language/tests/moments-test.c

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