shapiro1.c

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "local_proto.h"


double *shapiro_wilk(double *x, int n)
{
    static double y[2];
    double a[25], s2, *xcopy;
    double sumb = 0.0, sumx = 0.0, sumx2 = 0.0;
    int i, k;

    if ((xcopy = (double *)malloc(n * sizeof(double))) == NULL) {
        fprintf(stderr, "Memory error in shapiro_wilk\n");
        exit(EXIT_FAILURE);
    }

    k = n / 2;
    for (i = 0; i < n; ++i) {
        xcopy[i] = x[i];
        sumx += x[i];
        sumx2 += x[i] * x[i];
    }
    s2 = sumx2 - sumx * sumx / n;

    qsort(xcopy, n, sizeof(double), dcmp);

    if (n == 3) {
        a[0] = (double).7071;
    }
    else if (n == 4) {
        a[0] = (double).6872;
        a[1] = (double).1677;
    }
    else if (n == 5) {
        a[0] = (double).6646;
        a[1] = (double).2413;
    }
    else if (n == 6) {
        a[0] = (double).6431;
        a[1] = (double).2806;
        a[2] = (double).0875;
    }
    else if (n == 7) {
        a[0] = (double).6233;
        a[1] = (double).3031;
        a[2] = (double).1401;
    }
    else if (n == 8) {
        a[0] = (double).6052;
        a[1] = (double).3164;
        a[2] = (double).1743;
        a[3] = (double).0561;
    }
    else if (n == 9) {
        a[0] = (double).5888;
        a[1] = (double).3244;
        a[2] = (double).1976;
        a[3] = (double).0947;
    }
    else if (n == 10) {
        a[0] = (double).5739;
        a[1] = (double).3291;
        a[2] = (double).2141;
        a[3] = (double).1224;
        a[4] = (double).0399;
    }
    else if (n == 11) {
        a[0] = (double).5601;
        a[1] = (double).3315;
        a[2] = (double).226;
        a[3] = (double).1429;
        a[4] = (double).0695;
    }
    else if (n == 12) {
        a[0] = (double).5475;
        a[1] = (double).3325;
        a[2] = (double).2347;
        a[3] = (double).1586;
        a[4] = (double).0922;
        a[5] = (double).0303;
    }
    else if (n == 13) {
        a[0] = (double).5359;
        a[1] = (double).3325;
        a[2] = (double).2412;
        a[3] = (double).1707;
        a[4] = (double).1099;
        a[5] = (double).0539;
    }
    else if (n == 14) {
        a[0] = (double).5251;
        a[1] = (double).3318;
        a[2] = (double).246;
        a[3] = (double).1802;
        a[4] = (double).124;
        a[5] = (double).0727;
        a[6] = (double).024;
    }
    else if (n == 15) {
        a[0] = (double).515;
        a[1] = (double).3306;
        a[2] = (double).2495;
        a[3] = (double).1878;
        a[4] = (double).1353;
        a[5] = (double).088;
        a[6] = (double).0433;
    }
    else if (n == 16) {
        a[0] = (double).5056;
        a[1] = (double).329;
        a[2] = (double).2521;
        a[3] = (double).1939;
        a[4] = (double).1447;
        a[5] = (double).1005;
        a[6] = (double).0593;
        a[7] = (double).0196;
    }
    else if (n == 17) {
        a[0] = (double).4968;
        a[1] = (double).3273;
        a[2] = (double).254;
        a[3] = (double).1988;
        a[4] = (double).1524;
        a[5] = (double).1109;
        a[6] = (double).0725;
        a[7] = (double).0359;
    }
    else if (n == 18) {
        a[0] = (double).4886;
        a[1] = (double).3253;
        a[2] = (double).2553;
        a[3] = (double).2027;
        a[4] = (double).1587;
        a[5] = (double).1197;
        a[6] = (double).0837;
        a[7] = (double).0496;
        a[8] = (double).0163;
    }
    else if (n == 19) {
        a[0] = (double).4808;
        a[1] = (double).3232;
        a[2] = (double).2561;
        a[3] = (double).2059;
        a[4] = (double).1641;
        a[5] = (double).1271;
        a[6] = (double).0932;
        a[7] = (double).0612;
        a[8] = (double).0303;
    }
    else if (n == 20) {
        a[0] = (double).4734;
        a[1] = (double).3211;
        a[2] = (double).2565;
        a[3] = (double).2085;
        a[4] = (double).1686;
        a[5] = (double).1334;
        a[6] = (double).1013;
        a[7] = (double).0711;
        a[8] = (double).0422;
        a[9] = (double).014;
    }
    else if (n == 21) {
        a[0] = (double).4643;
        a[1] = (double).3185;
        a[2] = (double).2578;
        a[3] = (double).2119;
        a[4] = (double).1736;
        a[5] = (double).1399;
        a[6] = (double).1092;
        a[7] = (double).0804;
        a[8] = (double).053;
        a[9] = (double).0263;
    }
    else if (n == 22) {
        a[0] = (double).459;
        a[1] = (double).3156;
        a[2] = (double).2571;
        a[3] = (double).2131;
        a[4] = (double).1764;
        a[5] = (double).1443;
        a[6] = (double).115;
        a[7] = (double).0878;
        a[8] = (double).0618;
        a[9] = (double).0368;
        a[10] = (double).0122;
    }
    else if (n == 23) {
        a[0] = (double).4542;
        a[1] = (double).3126;
        a[2] = (double).2563;
        a[3] = (double).2139;
        a[4] = (double).1787;
        a[5] = (double).148;
        a[6] = (double).1201;
        a[7] = (double).0941;
        a[8] = (double).0696;
        a[9] = (double).0459;
        a[10] = (double).0228;
    }
    else if (n == 24) {
        a[0] = (double).4493;
        a[1] = (double).3098;
        a[2] = (double).2554;
        a[3] = (double).2145;
        a[4] = (double).1807;
        a[5] = (double).1512;
        a[6] = (double).1245;
        a[7] = (double).0997;
        a[8] = (double).0764;
        a[9] = (double).0539;
        a[10] = (double).0321;
        a[11] = (double).0107;
    }
    else if (n == 25) {
        a[0] = (double).445;
        a[1] = (double).3069;
        a[2] = (double).2543;
        a[3] = (double).2148;
        a[4] = (double).1822;
        a[5] = (double).1539;
        a[6] = (double).1283;
        a[7] = (double).1046;
        a[8] = (double).0823;
        a[9] = (double).061;
        a[10] = (double).0403;
        a[11] = (double).02;
    }
    else if (n == 26) {
        a[0] = (double).4407;
        a[1] = (double).3043;
        a[2] = (double).2533;
        a[3] = (double).2151;
        a[4] = (double).1836;
        a[5] = (double).1563;
        a[6] = (double).1316;
        a[7] = (double).1089;
        a[8] = (double).0876;
        a[9] = (double).0672;
        a[10] = (double).0476;
        a[11] = (double).0284;
        a[12] = (double).0094;
    }
    else if (n == 27) {
        a[0] = (double).4366;
        a[1] = (double).3018;
        a[2] = (double).2522;
        a[3] = (double).2152;
        a[4] = (double).1848;
        a[5] = (double).1584;
        a[6] = (double).1346;
        a[7] = (double).1128;
        a[8] = (double).0923;
        a[9] = (double).0728;
        a[10] = (double).054;
        a[11] = (double).0358;
        a[12] = (double).0178;
    }
    else if (n == 28) {
        a[0] = (double).4328;
        a[1] = (double).2992;
        a[2] = (double).251;
        a[3] = (double).2151;
        a[4] = (double).1857;
        a[5] = (double).1601;
        a[6] = (double).1372;
        a[7] = (double).1162;
        a[8] = (double).0965;
        a[9] = (double).0778;
        a[10] = (double).0598;
        a[11] = (double).0424;
        a[12] = (double).0253;
        a[13] = (double).0084;
    }
    else if (n == 29) {
        a[0] = (double).4291;
        a[1] = (double).2968;
        a[2] = (double).2499;
        a[3] = (double).215;
        a[4] = (double).1864;
        a[5] = (double).1616;
        a[6] = (double).1395;
        a[7] = (double).1192;
        a[8] = (double).1002;
        a[9] = (double).0822;
        a[10] = (double).065;
        a[11] = (double).0483;
        a[12] = (double).032;
        a[13] = (double).0159;
    }
    else if (n == 30) {
        a[0] = (double).4254;
        a[1] = (double).2944;
        a[2] = (double).2487;
        a[3] = (double).2148;
        a[4] = (double).187;
        a[5] = (double).163;
        a[6] = (double).1415;
        a[7] = (double).1219;
        a[8] = (double).1036;
        a[9] = (double).0862;
        a[10] = (double).0697;
        a[11] = (double).0537;
        a[12] = (double).0381;
        a[13] = (double).0227;
        a[14] = (double).0076;
    }
    else if (n == 31) {
        a[0] = (double).422;
        a[1] = (double).2921;
        a[2] = (double).2475;
        a[3] = (double).2145;
        a[4] = (double).1874;
        a[5] = (double).1641;
        a[6] = (double).1433;
        a[7] = (double).1243;
        a[8] = (double).1066;
        a[9] = (double).0899;
        a[10] = (double).0739;
        a[11] = (double).0585;
        a[12] = (double).0435;
        a[13] = (double).0289;
        a[14] = (double).0144;
    }
    else if (n == 32) {
        a[0] = (double).4188;
        a[1] = (double).2898;
        a[2] = (double).2463;
        a[3] = (double).2141;
        a[4] = (double).1878;
        a[5] = (double).1651;
        a[6] = (double).1449;
        a[7] = (double).1265;
        a[8] = (double).1093;
        a[9] = (double).0931;
        a[10] = (double).0777;
        a[11] = (double).0629;
        a[12] = (double).0485;
        a[13] = (double).0344;
        a[14] = (double).0206;
        a[15] = (double).0068;
    }
    else if (n == 33) {
        a[0] = (double).4156;
        a[1] = (double).2876;
        a[2] = (double).2451;
        a[3] = (double).2137;
        a[4] = (double).188;
        a[5] = (double).166;
        a[6] = (double).1463;
        a[7] = (double).1284;
        a[8] = (double).1118;
        a[9] = (double).0961;
        a[10] = (double).0812;
        a[11] = (double).0669;
        a[12] = (double).053;
        a[13] = (double).0395;
        a[14] = (double).0262;
        a[15] = (double).0131;
    }
    else if (n == 34) {
        a[0] = (double).4127;
        a[1] = (double).2854;
        a[2] = (double).2439;
        a[3] = (double).2132;
        a[4] = (double).1882;
        a[5] = (double).1667;
        a[6] = (double).1475;
        a[7] = (double).1301;
        a[8] = (double).114;
        a[9] = (double).0988;
        a[10] = (double).0844;
        a[11] = (double).0706;
        a[12] = (double).0572;
        a[13] = (double).0441;
        a[14] = (double).0314;
        a[15] = (double).0187;
        a[16] = (double).0062;
    }
    else if (n == 35) {
        a[0] = (double).4096;
        a[1] = (double).2834;
        a[2] = (double).2427;
        a[3] = (double).2127;
        a[4] = (double).1883;
        a[5] = (double).1673;
        a[6] = (double).1487;
        a[7] = (double).1317;
        a[8] = (double).116;
        a[9] = (double).1013;
        a[10] = (double).0873;
        a[11] = (double).0739;
        a[12] = (double).061;
        a[13] = (double).0484;
        a[14] = (double).0361;
        a[15] = (double).0239;
        a[16] = (double).0119;
    }
    else if (n == 36) {
        a[0] = (double).4068;
        a[1] = (double).2813;
        a[2] = (double).2415;
        a[3] = (double).2121;
        a[4] = (double).1883;
        a[5] = (double).1678;
        a[6] = (double).1496;
        a[7] = (double).1331;
        a[8] = (double).1179;
        a[9] = (double).1036;
        a[10] = (double).09;
        a[11] = (double).077;
        a[12] = (double).0645;
        a[13] = (double).0523;
        a[14] = (double).0404;
        a[15] = (double).0287;
        a[16] = (double).0172;
        a[17] = (double).0057;
    }
    else if (n == 37) {
        a[0] = (double).404;
        a[1] = (double).2794;
        a[2] = (double).2403;
        a[3] = (double).2116;
        a[4] = (double).1883;
        a[5] = (double).1683;
        a[6] = (double).1505;
        a[7] = (double).1344;
        a[8] = (double).1196;
        a[9] = (double).1056;
        a[10] = (double).0924;
        a[11] = (double).0798;
        a[12] = (double).0677;
        a[13] = (double).0559;
        a[14] = (double).0444;
        a[15] = (double).0331;
        a[16] = (double).022;
        a[17] = (double).011;
    }
    else if (n == 38) {
        a[0] = (double).4015;
        a[1] = (double).2774;
        a[2] = (double).2391;
        a[3] = (double).211;
        a[4] = (double).1881;
        a[5] = (double).1686;
        a[6] = (double).1513;
        a[7] = (double).1356;
        a[8] = (double).1211;
        a[9] = (double).1075;
        a[10] = (double).0947;
        a[11] = (double).0824;
        a[12] = (double).0706;
        a[13] = (double).0592;
        a[14] = (double).0481;
        a[15] = (double).0372;
        a[16] = (double).0264;
        a[17] = (double).0158;
        a[18] = (double).0053;
    }
    else if (n == 39) {
        a[0] = (double).3989;
        a[1] = (double).2755;
        a[2] = (double).238;
        a[3] = (double).2104;
        a[4] = (double).188;
        a[5] = (double).1689;
        a[6] = (double).152;
        a[7] = (double).1366;
        a[8] = (double).1225;
        a[9] = (double).1092;
        a[10] = (double).0967;
        a[11] = (double).0848;
        a[12] = (double).0733;
        a[13] = (double).0622;
        a[14] = (double).0515;
        a[15] = (double).0409;
        a[16] = (double).0305;
        a[17] = (double).0203;
        a[18] = (double).0101;
    }
    else if (n == 40) {
        a[0] = (double).3964;
        a[1] = (double).2737;
        a[2] = (double).2368;
        a[3] = (double).2098;
        a[4] = (double).1878;
        a[5] = (double).1691;
        a[6] = (double).1526;
        a[7] = (double).1376;
        a[8] = (double).1237;
        a[9] = (double).1108;
        a[10] = (double).0986;
        a[11] = (double).087;
        a[12] = (double).0759;
        a[13] = (double).0651;
        a[14] = (double).0546;
        a[15] = (double).0444;
        a[16] = (double).0343;
        a[17] = (double).0244;
        a[18] = (double).0146;
        a[19] = (double).0049;
    }
    else if (n == 41) {
        a[0] = (double).394;
        a[1] = (double).2719;
        a[2] = (double).2357;
        a[3] = (double).2091;
        a[4] = (double).1876;
        a[5] = (double).1693;
        a[6] = (double).1531;
        a[7] = (double).1384;
        a[8] = (double).1249;
        a[9] = (double).1123;
        a[10] = (double).1004;
        a[11] = (double).0891;
        a[12] = (double).0782;
        a[13] = (double).0677;
        a[14] = (double).0575;
        a[15] = (double).0476;
        a[16] = (double).0379;
        a[17] = (double).0283;
        a[18] = (double).0188;
        a[19] = (double).0094;
    }
    else if (n == 42) {
        a[0] = (double).3917;
        a[1] = (double).2701;
        a[2] = (double).2345;
        a[3] = (double).2085;
        a[4] = (double).1874;
        a[5] = (double).1694;
        a[6] = (double).1535;
        a[7] = (double).1392;
        a[8] = (double).1259;
        a[9] = (double).1136;
        a[10] = (double).102;
        a[11] = (double).0909;
        a[12] = (double).0804;
        a[13] = (double).0701;
        a[14] = (double).0602;
        a[15] = (double).0506;
        a[16] = (double).0411;
        a[17] = (double).0318;
        a[18] = (double).0227;
        a[19] = (double).0136;
        a[20] = (double).0045;
    }
    else if (n == 43) {
        a[0] = (double).3894;
        a[1] = (double).2684;
        a[2] = (double).2334;
        a[3] = (double).2078;
        a[4] = (double).1871;
        a[5] = (double).1695;
        a[6] = (double).1539;
        a[7] = (double).1398;
        a[8] = (double).1269;
        a[9] = (double).1149;
        a[10] = (double).1035;
        a[11] = (double).0927;
        a[12] = (double).0824;
        a[13] = (double).0724;
        a[14] = (double).0628;
        a[15] = (double).0534;
        a[16] = (double).0442;
        a[17] = (double).0352;
        a[18] = (double).0263;
        a[19] = (double).0175;
        a[20] = (double).0087;
    }
    else if (n == 44) {
        a[0] = (double).3872;
        a[1] = (double).2667;
        a[2] = (double).2323;
        a[3] = (double).2072;
        a[4] = (double).1868;
        a[5] = (double).1695;
        a[6] = (double).1542;
        a[7] = (double).1405;
        a[8] = (double).1278;
        a[9] = (double).116;
        a[10] = (double).1049;
        a[11] = (double).0943;
        a[12] = (double).0842;
        a[13] = (double).0745;
        a[14] = (double).0651;
        a[15] = (double).056;
        a[16] = (double).0471;
        a[17] = (double).0383;
        a[18] = (double).0296;
        a[19] = (double).0211;
        a[20] = (double).0126;
        a[21] = (double).0042;
    }
    else if (n == 45) {
        a[0] = (double).385;
        a[1] = (double).2651;
        a[2] = (double).2313;
        a[3] = (double).2065;
        a[4] = (double).1865;
        a[5] = (double).1695;
        a[6] = (double).1545;
        a[7] = (double).141;
        a[8] = (double).1286;
        a[9] = (double).117;
        a[10] = (double).1062;
        a[11] = (double).0959;
        a[12] = (double).086;
        a[13] = (double).0765;
        a[14] = (double).0673;
        a[15] = (double).0584;
        a[16] = (double).0497;
        a[17] = (double).0412;
        a[18] = (double).0328;
        a[19] = (double).0245;
        a[20] = (double).0163;
        a[21] = (double).0081;
    }
    else if (n == 46) {
        a[0] = (double).383;
        a[1] = (double).2635;
        a[2] = (double).2302;
        a[3] = (double).2058;
        a[4] = (double).1862;
        a[5] = (double).1695;
        a[6] = (double).1548;
        a[7] = (double).1415;
        a[8] = (double).1293;
        a[9] = (double).118;
        a[10] = (double).1073;
        a[11] = (double).0972;
        a[12] = (double).0876;
        a[13] = (double).0783;
        a[14] = (double).0694;
        a[15] = (double).0607;
        a[16] = (double).0522;
        a[17] = (double).0439;
        a[18] = (double).0357;
        a[19] = (double).0277;
        a[20] = (double).0197;
        a[21] = (double).0118;
        a[22] = (double).0039;
    }
    else if (n == 47) {
        a[0] = (double).3808;
        a[1] = (double).262;
        a[2] = (double).2291;
        a[3] = (double).2052;
        a[4] = (double).1859;
        a[5] = (double).1695;
        a[6] = (double).155;
        a[7] = (double).142;
        a[8] = (double).13;
        a[9] = (double).1189;
        a[10] = (double).1085;
        a[11] = (double).0986;
        a[12] = (double).0892;
        a[13] = (double).0801;
        a[14] = (double).0713;
        a[15] = (double).0628;
        a[16] = (double).0546;
        a[17] = (double).0465;
        a[18] = (double).0385;
        a[19] = (double).0307;
        a[20] = (double).0229;
        a[21] = (double).0153;
        a[22] = (double).0076;
    }
    else if (n == 48) {
        a[0] = (double).3789;
        a[1] = (double).2604;
        a[2] = (double).2281;
        a[3] = (double).2045;
        a[4] = (double).1855;
        a[5] = (double).1693;
        a[6] = (double).1551;
        a[7] = (double).1423;
        a[8] = (double).1306;
        a[9] = (double).1197;
        a[10] = (double).1095;
        a[11] = (double).0998;
        a[12] = (double).0906;
        a[13] = (double).0817;
        a[14] = (double).0731;
        a[15] = (double).0648;
        a[16] = (double).0568;
        a[17] = (double).0489;
        a[18] = (double).0411;
        a[19] = (double).0335;
        a[20] = (double).0259;
        a[21] = (double).0185;
        a[22] = (double).0111;
        a[23] = (double).0037;
    }
    else if (n == 49) {
        a[0] = (double).377;
        a[1] = (double).2589;
        a[2] = (double).2271;
        a[3] = (double).2038;
        a[4] = (double).1851;
        a[5] = (double).1692;
        a[6] = (double).1553;
        a[7] = (double).1427;
        a[8] = (double).1312;
        a[9] = (double).1205;
        a[10] = (double).1105;
        a[11] = (double).101;
        a[12] = (double).0919;
        a[13] = (double).0832;
        a[14] = (double).0748;
        a[15] = (double).0667;
        a[16] = (double).0588;
        a[17] = (double).0511;
        a[18] = (double).0436;
        a[19] = (double).0361;
        a[20] = (double).0288;
        a[21] = (double).0215;
        a[22] = (double).0143;
        a[23] = (double).0071;
    }
    else if (n == 50) {
        a[0] = (double).3751;
        a[1] = (double).2574;
        a[2] = (double).226;
        a[3] = (double).2032;
        a[4] = (double).1847;
        a[5] = (double).1691;
        a[6] = (double).1554;
        a[7] = (double).143;
        a[8] = (double).1317;
        a[9] = (double).1212;
        a[10] = (double).1113;
        a[11] = (double).102;
        a[12] = (double).0932;
        a[13] = (double).0846;
        a[14] = (double).0764;
        a[15] = (double).0685;
        a[16] = (double).0608;
        a[17] = (double).0532;
        a[18] = (double).0459;
        a[19] = (double).0386;
        a[20] = (double).0314;
        a[21] = (double).0244;
        a[22] = (double).0174;
        a[23] = (double).0104;
        a[24] = (double).0035;
    }

    if (n > 50 || n < 3) {
#ifdef NOISY
        fprintf(stdout,
                "  THIS IS THE SHAPIRO-WILK TEST FOR SMALL SAMPLES\n");
        fprintf(stdout,
                "  THE SAMPLE SIZE MUST BE LESS THAN OR EQUAL TO 50\n");
#endif /* NOISY */

        y[0] = y[1] = 0.0;
    }
    else {
        for (i = 1; i <= k; ++i)
            sumb += a[i - 1] * (x[n - i + 1] - x[i]);

        y[0] = sumb * sumb / s2;
        y[1] = s2;

#ifdef NOISY
        fprintf(stdout, "  TEST13 SW(N)  =%10.4f\n", y[0]);
#endif /* NOISY */
    }

    free(xcopy);

    return y;
}