src/math/percentiles.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2008, 2009 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 "percentiles.h"
  19 #include <math/order-stats.h>
  20
  21
  22 #include "gettext.h"
  23 #define _(msgid) gettext (msgid)
  24 #define N_(msgid) msgid
  25
  26 #include <libpspp/assertion.h>
  27 #include <libpspp/cast.h>
  28 #include <data/val-type.h>
  29 #include <gl/xalloc.h>
  30 #include <data/variable.h>
  31 #include <data/casereader.h>
  32
  33
  34 const char *const ptile_alg_desc[] = {
  35   "",
  36   N_("HAverage"),
  37   N_("Weighted Average"),
  38   N_("Rounded"),
  39   N_("Empirical"),
  40   N_("Empirical with averaging")
  41 };
  42
  43
  44
  45 double
  46 percentile_calculate (const struct percentile *ptl, enum pc_alg alg)
  47 {
  48   struct percentile *mutable = CONST_CAST (struct percentile *, ptl);
  49   const struct order_stats *os = &ptl->parent;
  50
  51   assert (os->cc == ptl->w);
  52
  53   if ( ptl->g1 == SYSMIS)
  54     mutable->g1 = (os->k[0].tc - os->k[0].cc) / os->k[0].c_p1;
  55
  56   if ( ptl->g1_star == SYSMIS)
  57     mutable->g1_star = os->k[0].tc - os->k[0].cc;
  58
  59   if ( ptl->g2 == SYSMIS)
  60     {
  61       if ( os->k[1].c == 0 )
  62         mutable->g2 = os->k[1].tc / os->k[1].c_p1;
  63       else if ( os->k[1].c_p1 == 0 )
  64         mutable->g2 = 0;
  65       else
  66         mutable->g2 = (os->k[1].tc - os->k[1].cc) / os->k[1].c_p1;
  67     }
  68
  69   if ( ptl->g2_star == SYSMIS)
  70     {
  71       if ( os->k[1].c == 0 )
  72         mutable->g2_star = os->k[1].tc;
  73       else if ( os->k[1].c_p1 == 0 )
  74         mutable->g2_star = 0;
  75       else
  76         mutable->g2_star = os->k[1].tc - os->k[1].cc;
  77     }
  78
  79   switch (alg)
  80     {
  81     case PC_WAVERAGE:
  82       if ( ptl->g1_star >= 1.0)
  83         return os->k[0].y_p1;
  84       else
  85         {
  86           double a = ( os->k[0].y == SYSMIS ) ? 0 : os->k[0].y;
  87
  88           if (os->k[0].c_p1 >= 1.0)
  89             return (1 - ptl->g1_star) * a + ptl->g1_star * os->k[0].y_p1;
  90           else
  91             return (1 - ptl->g1) * a + ptl->g1 * os->k[0].y_p1;
  92         }
  93       break;
  94
  95     case PC_ROUND:
  96       {
  97         double a = ( os->k[0].y == SYSMIS ) ? 0 : os->k[0].y;
  98
  99         if (os->k[0].c_p1 >= 1.0)
 100           return (ptl->g1_star < 0.5) ? a : os->k[0].y_p1;
 101         else
 102           return (ptl->g1 < 0.5) ? a : os->k[0].y_p1;
 103       }
 104       break;
 105
 106     case PC_EMPIRICAL:
 107       if ( ptl->g1_star == 0 )
 108         return os->k[0].y;
 109       else
 110         return os->k[0].y_p1;
 111       break;
 112
 113     case PC_HAVERAGE:
 114       if ( ptl->g2_star >= 1.0)
 115         {
 116           return os->k[1].y_p1;
 117         }
 118       else
 119         {
 120           double a = ( os->k[1].y == SYSMIS ) ? 0 : os->k[1].y;
 121
 122           if ( os->k[1].c_p1 >= 1.0)
 123             {
 124               if ( ptl->g2_star == 0)
 125                 return os->k[1].y;
 126
 127               return (1 - ptl->g2_star) * a + ptl->g2_star * os->k[1].y_p1;
 128             }
 129           else
 130             {
 131               return (1 - ptl->g2) * a + ptl->g2 * os->k[1].y_p1;
 132             }
 133         }
 134
 135       break;
 136
 137     case PC_AEMPIRICAL:
 138       if ( ptl->g1_star == 0 )
 139         return (os->k[0].y + os->k[0].y_p1)/ 2.0;
 140       else
 141         return os->k[0].y_p1;
 142       break;
 143
 144     default:
 145       NOT_REACHED ();
 146       break;
 147     }
 148
 149   NOT_REACHED ();
 150
 151   return SYSMIS;
 152 }
 153
 154
 155 static void
 156 destroy (struct statistic *stat)
 157 {
 158   struct percentile *ptl = UP_CAST (stat, struct percentile, parent.parent);
 159   struct order_stats *os = &ptl->parent;
 160   free (os->k);
 161   free (ptl);
 162 }
 163
 164
 165 struct percentile *
 166 percentile_create (double p, double W)
 167 {
 168   struct percentile *ptl = xzalloc (sizeof (*ptl));
 169   struct order_stats *os = &ptl->parent;
 170   struct statistic *stat = &os->parent;
 171
 172   assert (p >= 0);
 173   assert (p <= 1.0);
 174
 175   ptl->ptile = p;
 176   ptl->w = W;
 177
 178   os->n_k = 2;
 179   os->k = xcalloc (sizeof (*os->k), 2);
 180   os->k[0].tc = W * p;
 181   os->k[1].tc = (W + 1.0) * p;
 182
 183   ptl->g1 = ptl->g1_star = SYSMIS;
 184   ptl->g2 = ptl->g2_star = SYSMIS;
 185
 186   os->k[1].y_p1 = os->k[1].y = SYSMIS;
 187   os->k[0].y_p1 = os->k[0].y = SYSMIS;
 188
 189   stat->destroy = destroy;
 190
 191   return ptl;
 192 }
 193