1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2008 Free Software Foundation, Inc.
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.
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.
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/>. */
18 #include "percentiles.h"
19 #include <math/order-stats.h>
23 #define _(msgid) gettext (msgid)
24 #define N_(msgid) msgid
26 #include <libpspp/assertion.h>
27 #include <data/val-type.h>
28 #include <gl/xalloc.h>
29 #include <data/variable.h>
30 #include <data/casereader.h>
33 const char *const ptile_alg_desc[] = {
36 N_("Weighted Average"),
39 N_("Empirical with averaging")
45 percentile_calculate (const struct percentile *ptl, enum pc_alg alg)
47 struct percentile *mutable = (struct percentile *) ptl;
48 const struct order_stats *os = &ptl->parent;
50 if ( ptl->g1 == SYSMIS)
51 mutable->g1 = (os->k[0].tc - os->k[0].cc) / os->k[0].c_p1;
53 if ( ptl->g1_star == SYSMIS)
54 mutable->g1_star = os->k[0].tc - os->k[0].cc;
56 if ( ptl->g2 == SYSMIS)
58 if ( os->k[1].c == 0 )
59 mutable->g2 = os->k[1].tc / os->k[1].c_p1;
60 else if ( os->k[1].c_p1 == 0 )
63 mutable->g2 = (os->k[1].tc - os->k[1].cc) / os->k[1].c_p1;
66 if ( ptl->g2_star == SYSMIS)
68 if ( os->k[1].c == 0 )
69 mutable->g2_star = os->k[1].tc;
70 else if ( os->k[1].c_p1 == 0 )
73 mutable->g2_star = os->k[1].tc - os->k[1].cc;
79 if ( ptl->g1_star >= 1.0)
83 double a = ( os->k[0].y == SYSMIS ) ? 0 : os->k[0].y;
85 if (os->k[0].c_p1 >= 1.0)
86 return (1 - ptl->g1_star) * a + ptl->g1_star * os->k[0].y_p1;
88 return (1 - ptl->g1) * a + ptl->g1 * os->k[0].y_p1;
94 double a = ( os->k[0].y == SYSMIS ) ? 0 : os->k[0].y;
96 if (os->k[0].c_p1 >= 1.0)
97 return (ptl->g1_star < 0.5) ? a : os->k[0].y_p1;
99 return (ptl->g1 < 0.5) ? a : os->k[0].y_p1;
104 if ( ptl->g1_star == 0 )
107 return os->k[0].y_p1;
111 if ( ptl->g2_star >= 1.0)
113 return os->k[1].y_p1;
117 double a = ( os->k[1].y == SYSMIS ) ? 0 : os->k[1].y;
119 if ( os->k[1].c_p1 >= 1.0)
121 if ( ptl->g2_star == 0)
124 return (1 - ptl->g2_star) * a + ptl->g2_star * os->k[1].y_p1;
128 return (1 - ptl->g2) * a + ptl->g2 * os->k[1].y_p1;
135 if ( ptl->g1_star == 0 )
136 return (os->k[0].y + os->k[0].y_p1)/ 2.0;
138 return os->k[0].y_p1;
153 destroy (struct statistic *stat)
155 struct order_stats *os = (struct order_stats *) stat;
162 percentile_create (double p, double W)
164 struct percentile *ptl = xzalloc (sizeof (*ptl));
165 struct order_stats *os = (struct order_stats *) ptl;
166 struct statistic *stat = (struct statistic *) ptl;
174 os->k = xcalloc (sizeof (*os->k), 2);
176 os->k[1].tc = (W + 1.0) * p;
178 ptl->g1 = ptl->g1_star = SYSMIS;
179 ptl->g2 = ptl->g2_star = SYSMIS;
181 os->k[1].y_p1 = os->k[1].y = SYSMIS;
182 os->k[0].y_p1 = os->k[0].y = SYSMIS;
184 stat->destroy = destroy;
191 percentile_dump (const struct percentile *ptl)
193 printf ("Percentile %g:\n\tk1: ", ptl->ptile);
195 dump_os_k1 ((const struct os *)ptl);
197 dump_os_k2 ((const struct os *)ptl);