pintos-os.org Git - pspp/blob - src/math/shapiro-wilk.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 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 "math/shapiro-wilk.h"
  20
  21 #include <math.h>
  22 #include <assert.h>
  23 #include <gsl/gsl_cdf.h>
  24 #include "libpspp/cast.h"
  25 #include "libpspp/compiler.h"
  26 #include "libpspp/message.h"
  27
  28 #include "gettext.h"
  29 #define _(msgid) gettext (msgid)
  30 #define N_(msgid) msgid
  31
  32 /* Return the sum of coeff[i] * x^i for all i in the range [0,order).
  33    It is the caller's responsibility to ensure that coeff points to a
  34    large enough array containing the desired coefficients.  */
  35 static double
  36 polynomial (const double *coeff, int order, double x)
  37 {
  38   double result = 0;
  39   for (int i = 0; i < order; ++i)
  40     result += coeff[i] * pow (x, i);
  41
  42   return result;
  43 }
  44
  45 static double
  46 m_i (struct shapiro_wilk *sw, int i)
  47 {
  48   assert (i > 0);
  49   assert (i <= sw->n);
  50   double x = (i - 0.375) / (sw->n + 0.25);
  51
  52   return gsl_cdf_ugaussian_Pinv (x);
  53 }
  54
  55 static double
  56 a_i (struct shapiro_wilk *sw, int i)
  57 {
  58   assert (i > 0);
  59   assert (i <= sw->n);
  60
  61   if (i <  sw->n / 2.0)
  62     return -a_i (sw, sw->n - i + 1);
  63   else if (i == sw->n)
  64     return sw->a_n1;
  65   else if (i == sw->n - 1)
  66     return sw->a_n2;
  67   else
  68     return m_i (sw, i) / sqrt (sw->epsilon);
  69 }
  70
  71 struct ccase;
  72
  73 static void
  74 acc (struct statistic *s, const struct ccase *cx UNUSED, double c,
  75      double cc, double y)
  76 {
  77   struct shapiro_wilk *sw = UP_CAST (s, struct shapiro_wilk, parent.parent);
  78
  79   double int_part, frac_part;
  80   frac_part = modf (c, &int_part);
  81
  82   if (frac_part != 0 && !sw->warned)
  83     {
  84       msg (MW, N_ ("One or more weight values are non-integer."
  85                    "  Fractional parts will be ignored when calculating the Shapiro-Wilk statistic."));
  86       sw->warned = true;
  87     }
  88
  89   for (int i = 0; i < int_part; ++i)
  90   {
  91     double a_ii = a_i (sw, cc - c + i + 1);
  92     double x_ii = y;
  93
  94     sw->numerator += a_ii * x_ii;
  95     sw->denominator += pow (x_ii - sw->mean, 2);
  96   }
  97 }
  98
  99 static void
 100 destroy (struct statistic *s)
 101 {
 102   struct shapiro_wilk *sw = UP_CAST (s, struct shapiro_wilk, parent.parent);
 103   free (sw);
 104 }
 105
 106
 107
 108 double
 109 shapiro_wilk_calculate (const struct shapiro_wilk *sw)
 110 {
 111   return sw->numerator * sw->numerator / sw->denominator;
 112 }
 113
 114 /* Inititialse a SW ready for calculating the statistic for
 115    a dataset of size N.  */
 116 struct shapiro_wilk *
 117 shapiro_wilk_create (int n, double mean)
 118 {
 119   if (n < 3 || n > 5000)
 120     return NULL;
 121
 122   struct shapiro_wilk *sw = XZALLOC (struct shapiro_wilk);
 123   struct order_stats *os = &sw->parent;
 124   struct statistic *stat = &os->parent;
 125
 126   const double c1[] = {0, 0.221157,  -0.147981,
 127                        -2.071190, 4.434685, -2.706056};
 128
 129   const double c2[] = {0, 0.042981, -0.293762,
 130                        -1.752461, 5.682633, -3.582633};
 131
 132   sw->n = n;
 133
 134   const double u = 1.0 / sqrt (sw->n);
 135
 136   double m = 0;
 137   for (int i = 1; i <= sw->n; ++i)
 138     {
 139       double m_ii = m_i (sw, i);
 140       m += m_ii * m_ii;
 141     }
 142
 143   double m_n1 = m_i (sw, sw->n);
 144   double m_n2 = m_i (sw, sw->n - 1);
 145
 146   sw->a_n1 = polynomial (c1, 6, u);
 147   sw->a_n1 += m_n1 / sqrt (m);
 148
 149   sw->a_n2 = polynomial (c2, 6, u);
 150   sw->a_n2 += m_n2 / sqrt (m);
 151   sw->mean = mean;
 152
 153   sw->epsilon =  m - 2 * pow (m_n1, 2) - 2 * pow (m_n2, 2);
 154   sw->epsilon /= 1 - 2 * pow (sw->a_n1, 2) - 2 * pow (sw->a_n2, 2);
 155
 156   sw->warned = false;
 157
 158   stat->accumulate = acc;
 159   stat->destroy = destroy;
 160
 161   return sw;
 162 }
 163
 164 double
 165 shapiro_wilk_significance (double n, double w)
 166 {
 167   const double g[] = {-2.273, 0.459};
 168   const double c3[] = {.544,-.39978,.025054,-6.714e-4};
 169   const double c4[] = {1.3822,-.77857,.062767,-.0020322};
 170   const double c5[] = {-1.5861,-.31082,-.083751,.0038915 };
 171   const double c6[] = {-.4803,-.082676,.0030302};
 172
 173   double m, s;
 174   double y = log (1 - w);
 175   if (n == 3)
 176     {
 177       double pi6 = 6.0 / M_PI;
 178       double stqr = asin(sqrt (3/ 4.0));
 179       double p = pi6 * (asin (sqrt (w)) - stqr);
 180       return (p < 0) ? 0 : p;
 181     }
 182   else if (n <= 11)
 183     {
 184       double gamma = polynomial (g, 2, n);
 185       y = - log (gamma - y);
 186       m = polynomial (c3, 4, n);
 187       s = exp (polynomial (c4, 4, n));
 188     }
 189   else
 190     {
 191       double xx = log(n);
 192       m = polynomial (c5, 4, xx);
 193       s = exp (polynomial (c6, 3, xx));
 194     }
 195
 196   double p = gsl_cdf_gaussian_Q (y - m, s);
 197   return p;
 198 }