pintos-os.org Git - pspp/blob - src/language/stats/means-calc.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2011, 2012, 2013, 2019 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
  19 #include "data/case.h"
  20 #include "data/format.h"
  21 #include "data/variable.h"
  22
  23 #include "libpspp/bt.h"
  24 #include "libpspp/hmap.h"
  25 #include "libpspp/misc.h"
  26 #include "libpspp/pool.h"
  27
  28 #include "math/moments.h"
  29 #include "output/pivot-table.h"
  30
  31 #include <math.h>
  32
  33 #include "means.h"
  34
  35 #include "gettext.h"
  36 #define _(msgid) gettext (msgid)
  37 #define N_(msgid) (msgid)
  38
  39 /* A base struct for all statistics.  */
  40 struct statistic
  41 {
  42 };
  43
  44 /* Statistics which accumulate a single value.  */
  45 struct statistic_simple
  46 {
  47   struct statistic parent;
  48   double acc;
  49 };
  50
  51 /* Statistics based on moments.  */
  52 struct statistic_moment
  53 {
  54   struct statistic parent;
  55   struct moments1 *mom;
  56 };
  57
  58
  59 static struct statistic *
  60 default_create (struct pool *pool)
  61 {
  62   struct statistic_moment *pvd = pool_alloc (pool, sizeof *pvd);
  63
  64   pvd->mom = moments1_create (MOMENT_KURTOSIS);
  65
  66   return (struct statistic *) pvd;
  67 }
  68
  69 static void
  70 default_update (struct statistic *stat, double w, double x)
  71 {
  72   struct statistic_moment *pvd = (struct statistic_moment *) stat;
  73
  74   moments1_add (pvd->mom, x, w);
  75 }
  76
  77 static void
  78 default_destroy (struct statistic *stat)
  79 {
  80   struct statistic_moment *pvd = (struct statistic_moment *) stat;
  81   moments1_destroy (pvd->mom);
  82 }
  83
  84
  85 /* Simple statistics have nothing to destroy.  */
  86 static void
  87 simple_destroy (struct statistic *stat UNUSED)
  88 {
  89 }
  90
  91 \f
  92
  93 /* HARMONIC MEAN: The reciprocal of the sum of the reciprocals:
  94    1 / (1/(x_0) + 1/(x_1) + ... + 1/(x_{n-1})) */
  95
  96 struct harmonic_mean
  97 {
  98   struct statistic parent;
  99   double rsum;
 100   double n;
 101 };
 102
 103 static struct statistic *
 104 harmonic_create (struct pool *pool)
 105 {
 106   struct harmonic_mean *hm = pool_alloc (pool, sizeof *hm);
 107
 108   hm->rsum = 0;
 109   hm->n = 0;
 110
 111   return (struct statistic *) hm;
 112 }
 113
 114
 115 static void
 116 harmonic_update (struct statistic *stat, double w, double x)
 117 {
 118   struct harmonic_mean *hm = (struct harmonic_mean *) stat;
 119   hm->rsum  += w / x;
 120   hm->n += w;
 121 }
 122
 123
 124 static double
 125 harmonic_get (const struct statistic *pvd)
 126 {
 127   const struct harmonic_mean *hm = (const struct harmonic_mean *) pvd;
 128
 129   return hm->n / hm->rsum;
 130 }
 131
 132 \f
 133
 134 /* GEOMETRIC MEAN:  The nth root of the product of all n observations
 135    pow ((x_0 * x_1 * ... x_{n - 1}), 1/n)  */
 136 struct geometric_mean
 137 {
 138   struct statistic parent;
 139   double prod;
 140   double n;
 141 };
 142
 143 static struct statistic *
 144 geometric_create (struct pool *pool)
 145 {
 146   struct geometric_mean *gm = pool_alloc (pool, sizeof *gm);
 147
 148   gm->prod = 1.0;
 149   gm->n = 0;
 150
 151   return (struct statistic *) gm;
 152 }
 153
 154 static void
 155 geometric_update (struct statistic  *pvd, double w, double x)
 156 {
 157   struct geometric_mean *gm = (struct geometric_mean *)pvd;
 158   gm->prod  *=  pow (x, w);
 159   gm->n += w;
 160 }
 161
 162
 163 static double
 164 geometric_get (const struct statistic *pvd)
 165 {
 166   const struct geometric_mean *gm = (const struct geometric_mean *)pvd;
 167   return pow (gm->prod, 1.0 / gm->n);
 168 }
 169
 170 \f
 171
 172 /* The getters for moment based statistics simply calculate the
 173    moment.    The only exception is Std Dev. which needs to call
 174    sqrt as well.  */
 175
 176 static double
 177 sum_get (const struct statistic *pvd)
 178 {
 179   double n, mean;
 180
 181   moments1_calculate (((struct statistic_moment *)pvd)->mom, &n, &mean, 0, 0, 0);
 182
 183   return mean * n;
 184 }
 185
 186
 187 static double
 188 n_get (const struct statistic *pvd)
 189 {
 190   double n;
 191
 192   moments1_calculate (((struct statistic_moment *)pvd)->mom, &n, 0, 0, 0, 0);
 193
 194   return n;
 195 }
 196
 197 static double
 198 arithmean_get (const struct statistic *pvd)
 199 {
 200   double n, mean;
 201
 202   moments1_calculate (((struct statistic_moment *)pvd)->mom, &n, &mean, 0, 0, 0);
 203
 204   return mean;
 205 }
 206
 207 static double
 208 variance_get (const struct statistic *pvd)
 209 {
 210   double n, mean, variance;
 211
 212   moments1_calculate (((struct statistic_moment *)pvd)->mom, &n, &mean, &variance, 0, 0);
 213
 214   return variance;
 215 }
 216
 217
 218 static double
 219 stddev_get (const struct statistic *pvd)
 220 {
 221   return sqrt (variance_get (pvd));
 222 }
 223
 224 static double
 225 skew_get (const struct statistic *pvd)
 226 {
 227   double skew;
 228
 229   moments1_calculate (((struct statistic_moment *)pvd)->mom, NULL, NULL, NULL, &skew, 0);
 230
 231   return skew;
 232 }
 233
 234 static double
 235 sekurt_get (const struct statistic *pvd)
 236 {
 237   double n;
 238
 239   moments1_calculate (((struct statistic_moment *)pvd)->mom, &n, NULL, NULL, NULL, NULL);
 240
 241   return calc_sekurt (n);
 242 }
 243
 244 static double
 245 seskew_get (const struct statistic *pvd)
 246 {
 247   double n;
 248
 249   moments1_calculate (((struct statistic_moment *)pvd)->mom, &n, NULL, NULL, NULL, NULL);
 250
 251   return calc_seskew (n);
 252 }
 253
 254 static double
 255 kurt_get (const struct statistic *pvd)
 256 {
 257   double kurt;
 258
 259   moments1_calculate (((struct statistic_moment *)pvd)->mom, NULL, NULL, NULL, NULL, &kurt);
 260
 261   return kurt;
 262 }
 263
 264 static double
 265 semean_get (const struct statistic *pvd)
 266 {
 267   double n, var;
 268
 269   moments1_calculate (((struct statistic_moment *)pvd)->mom, &n, NULL, &var, NULL, NULL);
 270
 271   return sqrt (var / n);
 272 }
 273
 274 \f
 275
 276 /* MIN: The smallest (closest to minus infinity) value. */
 277
 278 static struct statistic *
 279 min_create (struct pool *pool)
 280 {
 281   struct statistic_simple *pvd = pool_alloc (pool, sizeof *pvd);
 282
 283   pvd->acc = DBL_MAX;
 284
 285   return (struct statistic *) pvd;
 286 }
 287
 288 static void
 289 min_update (struct statistic *pvd, double w UNUSED, double x)
 290 {
 291   double *r = &((struct statistic_simple *)pvd)->acc;
 292
 293   if (x < *r)
 294     *r = x;
 295 }
 296
 297 static double
 298 min_get (const struct statistic *pvd)
 299 {
 300   double *r = &((struct statistic_simple *)pvd)->acc;
 301
 302   return *r;
 303 }
 304
 305 /* MAX: The largest (closest to plus infinity) value. */
 306
 307 static struct statistic *
 308 max_create (struct pool *pool)
 309 {
 310   struct statistic_simple *pvd = pool_alloc (pool, sizeof *pvd);
 311
 312   pvd->acc = -DBL_MAX;
 313
 314   return (struct statistic *) pvd;
 315 }
 316
 317 static void
 318 max_update (struct statistic *pvd, double w UNUSED, double x)
 319 {
 320   double *r = &((struct statistic_simple *)pvd)->acc;
 321
 322   if (x > *r)
 323     *r = x;
 324 }
 325
 326 static double
 327 max_get (const struct statistic *pvd)
 328 {
 329   double *r = &((struct statistic_simple *)pvd)->acc;
 330
 331   return *r;
 332 }
 333
 334 \f
 335
 336 struct range
 337 {
 338   struct statistic parent;
 339   double min;
 340   double max;
 341 };
 342
 343 /* Initially min and max are set to their most (inverted) extreme possible
 344    values.  */
 345 static struct statistic *
 346 range_create (struct pool *pool)
 347 {
 348   struct range *r = pool_alloc (pool, sizeof *r);
 349
 350   r->min = DBL_MAX;
 351   r->max = -DBL_MAX;
 352
 353   return (struct statistic *) r;
 354 }
 355
 356 /* On each update, set min and max to X or leave unchanged,
 357    as appropriate.  */
 358 static void
 359 range_update (struct statistic *pvd, double w UNUSED, double x)
 360 {
 361   struct range *r = (struct range *) pvd;
 362
 363   if (x > r->max)
 364     r->max = x;
 365
 366   if (x < r->min)
 367     r->min = x;
 368 }
 369
 370 /*  Get the difference between min and max.  */
 371 static double
 372 range_get (const struct statistic *pvd)
 373 {
 374   const struct range *r = (struct range *) pvd;
 375
 376   return r->max - r->min;
 377 }
 378
 379 \f
 380
 381 /* LAST: The last value (the one closest to the end of the file).  */
 382
 383 static struct statistic *
 384 last_create (struct pool *pool)
 385 {
 386   struct statistic_simple *pvd = pool_alloc (pool, sizeof *pvd);
 387
 388   return (struct statistic *) pvd;
 389 }
 390
 391 static void
 392 last_update (struct statistic *pvd, double w UNUSED, double x)
 393 {
 394   struct statistic_simple *stat = (struct statistic_simple *) pvd;
 395
 396   stat->acc = x;
 397 }
 398
 399 static double
 400 last_get (const struct statistic *pvd)
 401 {
 402   const struct statistic_simple *stat = (struct statistic_simple *) pvd;
 403
 404   return stat->acc;
 405 }
 406
 407 /* FIRST: The first value (the one closest to the start of the file).  */
 408
 409 static struct statistic *
 410 first_create (struct pool *pool)
 411 {
 412   struct statistic_simple *pvd = pool_alloc (pool, sizeof *pvd);
 413
 414   pvd->acc = SYSMIS;
 415
 416   return (struct statistic *) pvd;
 417 }
 418
 419 static void
 420 first_update (struct statistic *pvd, double w UNUSED, double x)
 421 {
 422   struct statistic_simple *stat = (struct statistic_simple *) pvd;
 423
 424   if (stat->acc == SYSMIS)
 425     stat->acc = x;
 426 }
 427
 428 static double
 429 first_get (const struct statistic *pvd)
 430 {
 431   const struct statistic_simple *stat = (struct statistic_simple *) pvd;
 432
 433   return stat->acc;
 434 }
 435
 436 /* Table of cell_specs */
 437 const struct cell_spec cell_spec[n_MEANS_STATISTICS] = {
 438   {N_("Mean"),           "MEAN",      NULL          ,   default_create,   default_update,   arithmean_get, default_destroy},
 439   {N_("N"),              "COUNT",     PIVOT_RC_COUNT,   default_create,   default_update,   n_get,         default_destroy},
 440   {N_("Std. Deviation"), "STDDEV",    NULL          ,   default_create,   default_update,   stddev_get,    default_destroy},
 441 #if 0
 442   {N_("Median"),         "MEDIAN",    NULL          ,   default_create,   default_update,   NULL,          default_destroy},
 443   {N_("Group Median"),   "GMEDIAN",   NULL          ,   default_create,   default_update,   NULL,          default_destroy},
 444 #endif
 445   {N_("S.E. Mean"),      "SEMEAN",    NULL          ,   default_create,   default_update,   semean_get,    default_destroy},
 446   {N_("Sum"),            "SUM",       NULL          ,   default_create,   default_update,   sum_get,       default_destroy},
 447   {N_("Minimum"),        "MIN",       NULL          ,   min_create,       min_update,       min_get,       simple_destroy},
 448   {N_("Maximum"),        "MAX",       NULL          ,   max_create,       max_update,       max_get,       simple_destroy},
 449   {N_("Range"),          "RANGE",     NULL          ,   range_create,     range_update,     range_get,     simple_destroy},
 450   {N_("Variance"),       "VARIANCE",  PIVOT_RC_OTHER,   default_create,   default_update,   variance_get,  default_destroy},
 451   {N_("Kurtosis"),       "KURT",      PIVOT_RC_OTHER,   default_create,   default_update,   kurt_get,      default_destroy},
 452   {N_("S.E. Kurt"),      "SEKURT",    PIVOT_RC_OTHER,   default_create,   default_update,   sekurt_get,    default_destroy},
 453   {N_("Skewness"),       "SKEW",      PIVOT_RC_OTHER,   default_create,   default_update,   skew_get,      default_destroy},
 454   {N_("S.E. Skew"),      "SESKEW",    PIVOT_RC_OTHER,   default_create,   default_update,   seskew_get,    default_destroy},
 455   {N_("First"),          "FIRST",     NULL          ,   first_create,     first_update,     first_get,     simple_destroy},
 456   {N_("Last"),           "LAST",      NULL          ,   last_create,      last_update,      last_get,      simple_destroy},
 457 #if 0
 458   {N_("Percent N"),      "NPCT",      PIVOT_RC_PERCENT, default_create,   default_update,   NULL,          default_destroy},
 459   {N_("Percent Sum"),    "SPCT",      PIVOT_RC_PERCENT, default_create,   default_update,   NULL,          default_destroy},
 460 #endif
 461   {N_("Harmonic Mean"),  "HARMONIC",  NULL          ,   harmonic_create,  harmonic_update,  harmonic_get,  simple_destroy},
 462   {N_("Geom. Mean"),     "GEOMETRIC", NULL          ,   geometric_create, geometric_update, geometric_get, simple_destroy}
 463 };