for gls users

[Matthias Andree]

Greg Louis <glouis at dynamicro.on.ca> writes:

> The recent release (v1.4) of the GNU scientific library contains a new
> section with cumulative distribution function support.  I find it's
> still necessary to catch integration errors:

Ouch.

Here's what I've developed after digging through several maths books,
it's exact but may not be very fast. I've found it compares pretty well
to dcdflib.

# This is a shell archive.  Save it in a file, remove anything before
# this line, and then unpack it by entering "sh file".  Note, it may
# create directories; files and directories will be owned by you and
# have default permissions.
#
# This archive contains:
#
#	macdfchi.c
#	macdfchi.h
#
echo x - macdfchi.c
sed 's/^X//' >macdfchi.c << 'END-of-macdfchi.c'
X#include "macdfchi.h"
X#include <math.h>
X
X/* (C) 2003 by Matthias Andree */
X
X/* calculate complement of cumulative distribution function of chi^2
X * distribution, 1-P(x,nu). This function expects nu to be an even
X * number so that our incomplete Gamma function sum works! */
X
Xdouble macompcdfchi2(double x, int nu) {
X    int k;
X    /* halve arguments */
X    int robn = nu / 2;
X    double gx = x / 2;
X    double lngx = log(gx);
X    double lnqi = -gx;
X    double s = exp(-gx);
X
X    if (nu % 2) {
X	abort();
X    }
X
X    for (k = 1; k <= robn-1; k++) {
X	lnqi += lngx - log(k);
X	s += exp(lnqi);
X    }
X    return s;
X}
X
X#ifdef MAIN
X
X#include <stdio.h>
X#include <stdlib.h>
X
Xint main(int argc, char **argv) {
X    double q;
X
X    if (argc != 3) {
X	fprintf(stderr, "Usage: %s X ROBN\n", argv[0]);
X	exit(1);
X    }
X    q = macdfchi(atof(argv[1]), atof(argv[2]));
X    printf("q=%g\n", q);
X    exit(0);
X}
X#endif
END-of-macdfchi.c
echo x - macdfchi.h
sed 's/^X//' >macdfchi.h << 'END-of-macdfchi.h'
X#ifndef MACDFCHI_H
X#define MACDFCHI_H
X
Xextern double macompcdfchi2(double x, int robn); /* WARNING: robn is 2 * NU! */
X
X#endif
END-of-macdfchi.h
exit

-- 
Matthias Andree