1 /* PSPP - computes sample statistics.
2 Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
3 Written by Ben Pfaff <blp@gnu.org>.
5 This program is free software; you can redistribute it and/or
6 modify it under the terms of the GNU General Public License as
7 published by the Free Software Foundation; either version 2 of the
8 License, or (at your option) any later version.
10 This program is distributed in the hope that it will be useful, but
11 WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 General Public License for more details.
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
24 /* Returns the fourth power of its argument. */
32 /* Returns the cube of its argument. */
39 /* Returns the square of its argument. */
47 * kurtosis = [(n+1){n*sum(X**4) - 4*sum(X)*sum(X**3)
48 * + 6*sum(X)**2*sum(X**2)/n - 3*sum(X)**4/n**2}]
49 * /[(n-1)(n-2)(n-3)*(variance)**2]
53 * This and other formulas from _Biometry_, Sokal and Rohlf,
54 * W. H. Freeman and Company, 1969. See pages 117 and 136 especially.
57 calc_kurt (const double d[4], double n, double variance)
62 + 6.0 * sqr (d[0]) * d[1] / n
63 - 3.0 * hypercube (d[0]) / sqr (n)))
64 / ((n - 1.0) * (n - 2.0) * (n - 3.0) * sqr (variance))
65 - (3.0 * sqr (n - 1.0))
66 / ((n - 2.0) * (n - 3.)));
70 * standard error of kurtosis = sqrt([24n((n-1)**2)]/[(n-3)(n-2)(n+3)(n+5)])
73 calc_sekurt (double n)
75 return sqrt ((24.0 * n * sqr (n - 1.0))
76 / ((n - 3.0) * (n - 2.0) * (n + 3.0) * (n + 5.0)));
80 * skewness = [n*sum(X**3) - 3*sum(X)*sum(X**2) + 2*sum(X)**3/n]/
81 * /[(n-1)(n-2)*(variance)**3]
84 calc_skew (const double d[3], double n, double stddev)
87 ((n * d[2] - 3.0 * d[0] * d[1] + 2.0 * cube (d[0]) / n)
88 / ((n - 1.0) * (n - 2.0) * cube (stddev)));
92 * standard error of skewness = sqrt([6n(n-1)] / [(n-2)(n+1)(n+3)])
95 calc_seskew (double n)
98 sqrt ((6.0 * n * (n - 1.0))
99 / ((n - 2.0) * (n + 1.0) * (n + 3.0)));
102 /* Returns one-sided significance level corresponding to standard
103 normal deviate X. Algorithm from _SPSS Statistical Algorithms_,
107 normal_sig (double x)
109 const double a1 = .070523078;
110 const double a2 = .0422820123;
111 const double a3 = .0092705272;
112 const double a4 = .0001520143;
113 const double a5 = .0002765672;
114 const double a6 = .0000430638;
116 const double z = fabs (x) <= 14.14 ? 0.7071067812 * fabs (x) : 10.;
119 r = 1. + z * (a1 + z * (a2 + z * (a3 + z * (a4 + z * (a5 + z * a6)))));
122 r *= r; /* r ** 16 */
127 /* Taken from _BASIC Statistics: An Introduction to Problem Solving
128 with Your Personal Computer_, Jerry W. O'Dell, TAB 1984, page 314-5. */
130 normal_sig (double z)
134 h = 1 + 0.0498673470 * z;
136 h += 0.0211410061 * z;
138 h += 0.0032776263 * z;
140 h += 0.0000380036 * z;
142 h += 0.0000488906 * z;
144 h += 0.0000053830 * z;
145 return pow (h, -16.) / 2.;
149 /* Algorithm from _Turbo Pascal Programmer's Toolkit_, Rugg and
150 Feldman, Que 1989. Returns the significance level of chi-square
151 value CHISQ with DF degrees of freedom, correct to at least 7
154 chisq_sig (double x, int k)
156 if (x <= 0. || k < 1)
159 return 2. * normal_sig (sqrt (x));
162 double z, z_partial, term, denom, numerator, value;
170 z_partial *= x / term;
171 if (z_partial >= 10000000.)
175 while (z_partial >= 1.e-7);
176 denom = term = 2 - k % 2;
184 value = ((k + 1) / 2) * log (x) - x / 2.;
185 numerator = exp (value) * sqrt (2. / x / PI);
189 value = k / 2. * log (x) - x / 2.;
190 numerator = exp (value);
192 return 1. - numerator * z / denom;
196 double term, numer, norm_x;
199 numer = pow (x / k, 1. / 3.);
200 norm_x = numer / sqrt (term);
201 return 1.0 - normal_sig (norm_x);