1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2008, 2009, 2011 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/>. */
19 #include "math/percentiles.h"
21 #include "data/casereader.h"
22 #include "data/val-type.h"
23 #include "data/variable.h"
24 #include "libpspp/assertion.h"
25 #include "libpspp/cast.h"
26 #include "math/order-stats.h"
28 #include "gl/xalloc.h"
30 /* Return the value of the percentile. */
32 percentile_calculate (const struct percentile
*ptl
, enum pc_alg alg
)
34 struct percentile
*mutable = CONST_CAST (struct percentile
*, ptl
);
35 const struct order_stats
*os
= &ptl
->parent
;
37 if (ptl
->g1
== SYSMIS
)
38 mutable->g1
= (os
->k
[0].tc
- os
->k
[0].cc
) / os
->k
[0].c_p1
;
40 if (ptl
->g1_star
== SYSMIS
)
41 mutable->g1_star
= os
->k
[0].tc
- os
->k
[0].cc
;
43 if (ptl
->g2
== SYSMIS
)
46 mutable->g2
= os
->k
[1].tc
/ os
->k
[1].c_p1
;
47 else if (os
->k
[1].c_p1
== 0)
50 mutable->g2
= (os
->k
[1].tc
- os
->k
[1].cc
) / os
->k
[1].c_p1
;
53 if (ptl
->g2_star
== SYSMIS
)
56 mutable->g2_star
= os
->k
[1].tc
;
57 else if (os
->k
[1].c_p1
== 0)
60 mutable->g2_star
= os
->k
[1].tc
- os
->k
[1].cc
;
66 if (ptl
->g1_star
>= 1.0)
70 double a
= (os
->k
[0].y
== SYSMIS
) ? 0 : os
->k
[0].y
;
72 if (os
->k
[0].c_p1
>= 1.0)
73 return (1 - ptl
->g1_star
) * a
+ ptl
->g1_star
* os
->k
[0].y_p1
;
75 return (1 - ptl
->g1
) * a
+ ptl
->g1
* os
->k
[0].y_p1
;
81 double a
= (os
->k
[0].y
== SYSMIS
) ? 0 : os
->k
[0].y
;
83 if (os
->k
[0].c_p1
>= 1.0)
84 return (ptl
->g1_star
< 0.5) ? a
: os
->k
[0].y_p1
;
86 return (ptl
->g1
< 0.5) ? a
: os
->k
[0].y_p1
;
91 if (ptl
->g1_star
== 0)
98 if (ptl
->g2_star
>= 1.0)
100 return os
->k
[1].y_p1
;
104 double a
= (os
->k
[1].y
== SYSMIS
) ? 0 : os
->k
[1].y
;
106 if (os
->k
[1].c_p1
>= 1.0)
108 if (ptl
->g2_star
== 0)
111 return (1 - ptl
->g2_star
) * a
+ ptl
->g2_star
* os
->k
[1].y_p1
;
115 return (1 - ptl
->g2
) * a
+ ptl
->g2
* os
->k
[1].y_p1
;
122 if (ptl
->g1_star
== 0)
123 return (os
->k
[0].y
+ os
->k
[0].y_p1
)/ 2.0;
125 return os
->k
[0].y_p1
;
140 destroy (struct statistic
*stat
)
142 struct percentile
*ptl
= UP_CAST (stat
, struct percentile
, parent
.parent
);
147 /* Create the Pth percentile.
148 W is the total sum of weights in the data set.
151 percentile_create (double p
, double W
)
153 struct percentile
*ptl
= XZALLOC (struct percentile
);
154 struct order_stats
*os
= &ptl
->parent
;
155 struct statistic
*stat
= &os
->parent
;
166 os
->k
[1].tc
= (W
+ 1.0) * p
;
168 ptl
->g1
= ptl
->g1_star
= SYSMIS
;
169 ptl
->g2
= ptl
->g2_star
= SYSMIS
;
171 os
->k
[1].y_p1
= os
->k
[1].y
= SYSMIS
;
172 os
->k
[0].y_p1
= os
->k
[0].y
= SYSMIS
;
174 stat
->destroy
= destroy
;