gk.c: Goodman-Kruskal gamma calculator in C

gk.c

/* Copyright (c) 2012 Kyle Gorman
 * 
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to 
 * deal in the Software without restriction, including without limitation the 
 * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or 
 * sell copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
 * IN THE SOFTWARE.
 *
 * gk.c: Goodman-Kruskal gamma calculator
 * Kyle Gorman <[email protected]>
 * 
 * \gamma is defined as
 *
 * (C - D) / (C + D)
 *
 * where C is the number of concordant pairs and D the number of discordant 
 * pairs. The two-tailed statistical test for \gamma and $N$ observations, is
 *
 * \erfc(0.70710678 * |\frac{\gamma}{\sqrt{\frac{4 N + 10}{9 N (N - 1)}}}|
 *
 * where \erfc(x) is
 *
 * \frac{2}{\sqrt{\pi}} \int_x^{\infty} e^{-t^2} dt
 */


// compute two-sided test p-value
double pval(double gamma, int size) {
    return erfc(0.70710678 * fabs(gamma / sqrt((4. * size + 10.) /
                                               (9. * size * (size - 1.)))));
}

int main(int argc, char* argv[]) {

    if (argc != 3) {
        fprintf(stderr, "USAGE:\n\n\tgk C D\n\nAborting.\n");
        return 1;
    }

    int c = atoi(argv[1]);
    int d = atoi(argv[2]);
    int den = c + d;
    if (den < 1) {
        fprintf(stderr, "Invalid values.\n");
        return 1;
    }

    double gamma = (c - d) / (float) den;
    printf("gamma = %.03f\t", gamma);
    double p = pval(gamma, den);
    if (p < 0.001)
        printf("p = %01.01e\n", p);
    else
        printf("p = %01.04f\n", p);

    return 0;
}