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