src/language/tests/moments-test.c

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