home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / misc / volume7 / tot_info < prev    next >
Encoding:
Text File  |  1989-07-08  |  37.8 KB  |  1,755 lines
  1. Newsgroups: comp.sources.misc
  2. From: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc)
  3. Subject: v07i063: information statistic and chi-square for 2-D contingency tables
  4. Reply-To: gwyn@BRL.MIL
  5.  
  6. Posting-number: Volume 7, Issue 63
  7. Submitted-by: gwyn@BRL.MIL
  8. Archive-name: tot_info
  9.  
  10. #!/bin/sh
  11. # Self-unpacking archive format.  To unbundle, sh this file.
  12. echo 'Makefile' 1>&2
  13. cat >'Makefile' <<'END OF Makefile'
  14. #    Makefile -- -- makefile for "tot_info" utility
  15.  
  16. #    last edit:    89/02/06    D A Gwyn
  17.  
  18. #    SCCS ID:    @(#)tot_info.mk    1.1 (edited for publication)
  19.  
  20. PRODUCT = tot_info
  21.  
  22. MAKEFIL    = Makefile
  23. HEADERS = chisq.h gamma.h
  24. CFILES    = $(PRODUCT).c info.c gamma.c chisq.c # chisq.c is a "freebie"
  25. OBJECTS    = $(PRODUCT).o info.o gamma.o # chisq.o
  26. LIBTEST    = g_test
  27. TSRC    = $(LIBTEST).c
  28. TOBJ    = $(LIBTEST).o
  29. TFILES    = $(PRODUCT).in $(PRODUCT).exp
  30. LIBES    = -lm
  31. LLIBES    = $(LIBES)
  32. MANSRCS    = $(PRODUCT).1 chisq.3 gamma.3
  33.  
  34. # "standard" UNIX utilities that may need invocations tweaked:
  35. INS    = cp
  36. LINT    = lint
  37. FLOW    = cflow
  38. XREF    = cxref -c -s -w132
  39.  
  40. # printer configuration:
  41. LPR    = lp
  42. PR    = pr -f
  43.  
  44. # DWB configuration:
  45. TOPT    = -Ti300
  46. TMAC    = -man
  47. EQN    = eqn $(TOPT)
  48. TROFF    = troff $(TOPT) $(TMAC)
  49. TPOST    = | dimp
  50.  
  51.  
  52. $(PRODUCT):    $(OBJECTS)
  53.     $(CC) $(LDFLAGS) -o $@ $(OBJECTS) $(LIBES)
  54.  
  55. chisq.o info.o:    chisq.h
  56.  
  57. gamma.o:    gamma.h
  58.  
  59. print:        $(MAKEFIL) $(HEADERS) $(CFILES) $(TFILES) $(MANSRCS)
  60.     $(PR) $(MAKEFIL) $(HEADERS) $(CFILES) $(TFILES) $(MANSRCS) | $(LPR)
  61.  
  62. typeset:    $(MANSRCS)
  63.     $(EQN) $(PRODUCT).1 | $(TROFF) $(TPOST)
  64.     $(EQN) chisq.3 | $(TROFF) $(TPOST)
  65.     $(EQN) gamma.3 | $(TROFF) $(TPOST)
  66.  
  67. lint:        $(HEADERS) $(CFILES)
  68.     $(LINT) $(CFILES) $(LLIBES) > $@
  69.  
  70. flow:        $(HEADERS) $(CFILES)
  71.     $(FLOW) $(CFILES) > $@
  72.  
  73. xref:        $(HEADERS) $(CFILES)
  74.     $(XREF) $(CFILES) > $@
  75.  
  76. test:        gamma_test tot_test
  77.  
  78. tot_test:    $(PRODUCT) $(TFILES)
  79.     ./$(PRODUCT) < $(PRODUCT).in > $(PRODUCT).out
  80.     @if cmp -s $(PRODUCT).exp $(PRODUCT).out ; \
  81.     then    echo 'Tested okay' ; \
  82.     else    echo '*** Test failed; differences:' ; \
  83.         diff $(PRODUCT).exp $(PRODUCT).out ; \
  84.         exit 1 ; \
  85.     fi
  86.  
  87. gamma_test:    $(LIBTEST)
  88.     ./$(LIBTEST)
  89.  
  90. $(LIBTEST):    $(TOBJ) gamma.o
  91.     $(CC) $(LDFLAGS) -o $@ $(TOBJ) gamma.o $(LIBES)
  92.  
  93. $(LIBTEST).o:    gamma.h
  94.  
  95. clean:
  96.     -rm -f $(OBJECTS) $(TOBJ) $(LIBTEST) $(PRODUCT).out \
  97.     lint flow xref core mon.out
  98.  
  99. clobber:    clean
  100.     -rm -f $(PRODUCT)
  101. END OF Makefile
  102. echo 'chisq.3' 1>&2
  103. cat >'chisq.3' <<'END OF chisq.3'
  104. '\" e
  105. .TH CHISQ 3V VMB
  106. '\"    last edit:    88/09/19    D A Gwyn
  107. '\"    SCCS ID:    @(#)chisq.3    1.2 (edited for publication)
  108. '\"
  109. '\"    This source must be typeset via "eqn | troff -man"
  110. '\"
  111. .EQ
  112. delim @@
  113. .EN
  114. .SH NAME
  115. chisq \- independence tests for 2-way contingency tables
  116. .SH SYNOPSIS
  117. .ds cW (CW\" change to B (without the paren) if you don't have a CW font
  118. \f\*(cW
  119. #include "chisq.h"
  120. .br
  121. /* \fPThe following functions are declared in this header file.\f\*(cW */
  122. .sp
  123. double ChiSqTbl(int r, int c, const long *f, int *df);
  124. .br
  125. double InfoTbl(int r, int c, const long *f, int *df);\fP
  126. .SH DESCRIPTION
  127. \f\*(cWChiSqTbl\fP
  128. returns Pearson's @chi sup 2@ statistic
  129. for the \f\*(cWr\fP-by-\f\*(cWc\fP contingency table
  130. of observed frequencies of occurrence for each combination of
  131. row and column categories stored in row-major order
  132. in the vector starting at \f\*(cWf\fP.
  133. The computed number of degrees of freedom is stored into the location
  134. pointed to by \f\*(cWdf\fP.
  135. .br
  136. .EQ
  137. chi sup 2 ~==~ roman sum from i=0 to r-1 roman sum from j=0 to c-1
  138. { left ( { f sub ij ~-~ { f sub { i. } f sub { .j }} over N } right ) sup 2 }
  139. over { { f sub { i. } f sub { .j }} over N }
  140. .EN
  141. where the row and column sums are denoted
  142. @f sub { i. }~==~roman sum from j=0 to c-1 f sub ij@
  143. and
  144. @f sub { .j }~==~roman sum from i=0 to r-1 f sub ij@
  145. and the total number of occurrences is
  146. @N~==~roman sum from i=0 to r-1 roman sum from j=0 to c-1 f sub ij@.
  147. Terms involving zero row or column sum are omitted and the returned
  148. number of degrees of freedom are
  149. correspondingly reduced from the nominal value
  150. @(r-1)(c-1)@.
  151. .P
  152. \f\*(cWInfoTbl\fP
  153. returns 2 times Kullback's
  154. @"\v'-0.2'^\v'0.2'\h'-\w;^;u'I\^" ( H sub 1 : H sub 2 )@ statistic
  155. for a similar contingency table.
  156. The number of degrees of freedom to be used when relating this statistic
  157. to the @chi sup 2@ distribution is stored into the location
  158. pointed to by \f\*(cWdf\fP.
  159. (See \s-1DISCUSSION\s0 below.)
  160. .br
  161. .EQ
  162. 2 "\v'-0.2'^\v'0.2'\h'-\w;^;u'I\^" ( H sub 1 : H sub 2 ) ~==~
  163. 2 roman sum from i=0 to r-1 roman sum from j=0 to c-1
  164. { f sub ij log { { N f sub ij }
  165. over { f sub { i. } f sub { .j } } } }
  166. .EN
  167. where the row and column sums @f sub { i. }@ and @f sub { .j }@
  168. and the total number of observations @N@
  169. are the same as for Pearson's @chi sup 2@.
  170. .SH DISCUSSION
  171. @2 "\v'-0.2'^\v'0.2'\h'-\w;^;u'I\^" ( H sub 1 : H sub 2 )@
  172. represents twice the estimated information
  173. in favor of hypothesis @H sub 1@ over hypothesis @H sub 2@
  174. contained in the observed frequencies,
  175. where the \fInull hypothesis\fP @H sub 2@ is that the row and column
  176. categorizations are statistically independent and
  177. the \fIalternative hypothesis\fP @H sub 1@ is that they aren't.
  178. This statistic is asymptotically distributed as @chi sup 2@
  179. with @(r-1)(c-1)@ degrees of freedom;
  180. therefore by use of the \f\*(cWQChiSq\fP function (see \fIgamma\fP(3V))
  181. it can be converted to the probability that
  182. the statistic would be as large as was observed
  183. if the categorizations really are independent.
  184. This is of course the traditional use of Pearson's @chi sup 2@ statistic,
  185. which @2 "\v'-0.2'^\v'0.2'\h'-\w;^;u'I\^"@ approaches for large samples
  186. when @H sub 2@ is true.
  187. The main difference is that Pearson's @chi sup 2@ is not useful
  188. for small sample sizes,
  189. whereas @2 "\v'-0.2'^\v'0.2'\h'-\w;^;u'I\^"@ fully takes into account
  190. all available information for even the smallest samples.
  191. .P
  192. @2 "\v'-0.2'^\v'0.2'\h'-\w;^;u'I\^" ( H sub 1 : H sub 2 )@
  193. (along with its corresponding degrees of freedom
  194. for use with \f\*(cWQChiSq\fP) is \fIadditive\fP for independent samples,
  195. so that the information can be accumulated over the course of many
  196. independent experiments in order to improve one's ability
  197. to detect a violation of the null hypothesis.
  198. .SH CAVEATS
  199. Pearson's @chi sup 2@ test is unreliable for low frequencies;
  200. consequently it is usually recommended that
  201. categories be chosen to tally at least
  202. 5 occurrences in each bin.
  203. However, excessive combination of bins causes loss of information and
  204. consequently loss of discriminating power.
  205. Because @2 "\v'-0.2'^\v'0.2'\h'-\w;^;u'I\^"@ is correct for any sample size,
  206. it does not require combination of bins
  207. and is therefore immune from this problem.
  208. .bp
  209. .SH EXAMPLE
  210. \f\*(cW
  211. .ta 8n 16n 24n 32n 40n 48n 56n 64n
  212. .nf
  213. \&
  214. /*
  215.     Example program to read table dimensions, then frequencies,
  216.     and to print both statistics along with significance levels.
  217. */
  218. #include    <stdio.h>
  219. #include    <stdlib.h>            /* for EXIT_SUCCESS */
  220. #include    "chisq.h"
  221. #include    "gamma.h"            /* for QChiSq(\|) */
  222.  
  223. #define MAXTBL    1000
  224. static long    f[MAXTBL];        /* frequency tallies */
  225. static int    r;            /* # of rows */
  226. static int    c;            /* # of columns */
  227. #define x(i,j)    f[(i)*c+(j)]    /* convenient way to access freqs */
  228.  
  229. static void Calc(char *name, double (*func)(int, int, const long *, int *))
  230. {
  231.     int    df;            /* degrees of freedom for chi-square */
  232.     double    stat = (*func)(r, c, f, &df);        /* computed statistic */
  233.  
  234.     /* print results */
  235.     if ( stat >= 0.0 )
  236.         (void)printf("%s = %5.2f\etdf = %2d\etq = %7.4f\en",
  237.                  name, stat, df, QChiSq(stat, df)
  238.                 );
  239.     else
  240.         (void)printf(stat < -3.5 ? "out of memory\en"
  241.                : stat < -2.5 ? "table too small\en"
  242.                : stat < -1.5 ? "negative frequency\en"
  243.                : "too many zeros\en"
  244.                 );
  245. }
  246.  
  247. int main(int argc, char *argv[])
  248. {
  249.     register int    i;        /* row index */
  250.     register int    j;        /* column index */
  251.  
  252.     while ( scanf("%d %d\en", &r, &c) == 2 )    /* start new table */
  253.         {
  254.         /* input tallies */
  255.         for ( i = 0; i < r; ++i )
  256.             for ( j = 0; j < c; ++j )
  257.                 (void)scanf(" %ld", &x(i,j));
  258.  
  259.         /* compute statistics and print results */
  260.         Calc("chisq", ChiSqTbl);
  261.         Calc("2info", InfoTbl);
  262.         }
  263.     return EXIT_SUCCESS;
  264. }\fP
  265. .bp
  266. .SH FILES
  267. .ta 3i
  268. chisq.h    header file
  269. .br
  270. chisq.o, info.o    object modules
  271. .SH REFERENCE
  272. Solomon Kullback, \fIInformation Theory and Statistics\fP (Dover, 1968).
  273. .SH "SEE ALSO"
  274. gamma(3V).
  275. .SH DIAGNOSTICS
  276. Both these functions a return negative value for the statistic
  277. when it cannot be meaningfully computed, as follows:
  278. .RS
  279. .PD 0
  280. .TP
  281. \-1.0
  282. too many 0 entries in the table
  283. (for \f\*(cWChiSqTbl\fP, this means insufficient degrees of freedom;
  284. for \f\*(cWInfoTbl\fP, this means \fIevery\fP entry in the table was 0)
  285. .TP
  286. \-2.0
  287. some frequency was specified as less than 0
  288. .TP
  289. \-3.0
  290. specified number of rows or columns was less than 2
  291. .TP
  292. \-4.0
  293. unable to dynamically allocate enough working storage
  294. .PD
  295. .RE
  296. .SH AUTHOR
  297. Douglas A.\& Gwyn, BRL/VLD-VMB
  298. END OF chisq.3
  299. echo 'chisq.c' 1>&2
  300. cat >'chisq.c' <<'END OF chisq.c'
  301. /*
  302.     ChiSqTbl -- Pearson's chi-square test for a 2-way contingency table
  303.  
  304.     last edit:    88/09/19    D A Gwyn
  305.  
  306.     SCCS ID:    @(#)chisq.c    1.1 (edited for publication)
  307.  
  308.     Special return values:
  309.         -1.0    insufficient degrees of freedom (too many 0 entries)
  310.         -2.0    invalid table entry (frequency less than 0)
  311.         -3.0    invalid table dimensions (r or c less than 2)
  312.         -4.0    unable to allocate enough working storage
  313. */
  314.  
  315. #ifdef __STDC__
  316. #include    <stdlib.h>        /* malloc, free */
  317.  
  318. #include    "std.h"
  319. #else
  320. #include    "std.h"
  321.  
  322. extern pointer    malloc();
  323. extern void    free();
  324. #endif
  325.  
  326. double
  327. ChiSqTbl( r, c, f, pdf )
  328.     int        r;        /* # rows in table */
  329.     int        c;        /* # columns in table */
  330.     const long    *f;        /* -> r*c frequency tallies */
  331.     int        *pdf;        /* -> return # degrees of freedom */
  332.     {
  333. #define    x(i,j)    f[(i)*c+(j)]        /* convenient way to access freqs */
  334.     register int    i;        /* row index */
  335.     register int    j;        /* column index */
  336.     double        N;        /* (double)n */
  337.     double        chisq;        /* accumulates chi-square */
  338.     double        *xi;        /* row sums */
  339.     double        *xj;        /* col sums */
  340.     int        rdf = r - 1;    /* row degrees of freedom */
  341.     int        cdf = c - 1;    /* column degrees of freedom */
  342.  
  343.     if ( rdf <= 0 || cdf <= 0 )
  344.         {
  345.         chisq = -3.0;
  346.         goto ret3;
  347.         }
  348.  
  349.     if ( (xi = (double *)malloc( r * sizeof(double) )) == NULL )
  350.         {
  351.         chisq = -4.0;
  352.         goto ret3;
  353.         }
  354.  
  355.     if ( (xj = (double *)malloc( c * sizeof(double) )) == NULL )
  356.         {
  357.         chisq = -4.0;
  358.         goto ret2;
  359.         }
  360.  
  361.     /* compute row sums and total */
  362.  
  363.     N = 0.0;
  364.  
  365.     for ( i = 0; i < r; ++i )
  366.         {
  367.         double    sum = 0.0;    /* accumulator */
  368.  
  369.         for ( j = 0; j < c; ++j )
  370.             {
  371.             register long    k = x(i,j);
  372.  
  373.             if ( k < 0L )
  374.                 {
  375.                 chisq = -2.0;
  376.                 goto ret1;
  377.                 }
  378.  
  379.             sum += (double)k;
  380.             }
  381.  
  382.         if ( (xi[i] = sum) <= 0.0 )
  383.             --rdf;
  384.         else
  385.             N += sum;
  386.         }
  387.  
  388.     /* compute column sums */
  389.  
  390.     for ( j = 0; j < c; ++j )
  391.         {
  392.         double    sum = 0.0;    /* accumulator */
  393.  
  394.         for ( i = 0; i < r; ++i )
  395.             sum += (double)x(i,j);
  396.  
  397.         if ( (xj[j] = sum) <= 0.0 )
  398.             --cdf;
  399.         }
  400.  
  401.     if ( rdf <= 0 || cdf <= 0 )
  402.         {
  403.         chisq = -1.0;
  404.         goto ret1;
  405.         }
  406.  
  407.     *pdf = rdf * cdf;        /* total degrees of freedom */
  408.  
  409.     /* compute chi-square */
  410.  
  411.     chisq = 0.0;
  412.  
  413.     for ( i = 0; i < r; ++i )
  414.         {
  415.         double    pi = xi[i];    /* row sum */
  416.  
  417.         if ( pi > 0.0 )
  418.             for ( j = 0; j < c; ++j )
  419.                 {
  420.                 double    pj = xj[j];    /* column sum */
  421.  
  422.                 if ( pj > 0.0 )
  423.                     {
  424.                     double    expected = pi * pj / N;
  425.                     double    delta =
  426.                         (double)x(i,j) - expected;
  427.  
  428.                     chisq += delta * delta / expected;
  429.                     }
  430.                 }
  431.         }
  432.  
  433.     ret1:
  434.     free( (pointer)xj );
  435.     ret2:
  436.     free( (pointer)xi );
  437.     ret3:
  438.     return chisq;
  439.     }
  440. END OF chisq.c
  441. echo 'chisq.h' 1>&2
  442. cat >'chisq.h' <<'END OF chisq.h'
  443. /*
  444.     chisq.h -- definitions for contingency-table routines
  445.  
  446.     last edit:    88/09/19    D A Gwyn
  447.  
  448.     SCCS ID:    @(#)chisq.h    1.1 (edited for publication)
  449. */
  450.  
  451. /* library routines: */
  452.  
  453. #ifdef __STDC__
  454. extern double    ChiSqTbl( int r, int c, const long *f, int *pdf );
  455. extern double    InfoTbl( int r, int c, const long *f, int *pdf );
  456. #else
  457. extern double    ChiSqTbl();
  458. extern double    InfoTbl();
  459. #endif
  460. END OF chisq.h
  461. echo 'g_test.c' 1>&2
  462. cat >'g_test.c' <<'END OF g_test.c'
  463. /*
  464.     g_test -- tests for gamma and related functions
  465.  
  466.     last edit:    88/09/09    D A Gwyn
  467.     SCCS ID:    @(#)g_test.c    1.1 (edited for publication)
  468. */
  469.  
  470. #include    <stdio.h>
  471. #ifdef __STDC__
  472. #include    <stdlib.h>        /* for NULL etc. */
  473. #endif
  474.  
  475. #include    "std.h"
  476.  
  477. #include    "gamma.h"
  478.  
  479. #define    Printf    (void)printf
  480.  
  481. #define    TOL    1.0e-5            /* tolerance for checks */
  482.  
  483. static int    errs = 0;        /* tally errors */
  484.  
  485. static double
  486. RelDif( a, b )            /* returns relative difference:    */
  487.     double    a, b;        /* 0.0 if exactly the same,
  488.                    otherwise ratio of difference
  489.                    to the larger of the two    */
  490.     {
  491.     double    c = Abs( a );
  492.     double    d = Abs( b );
  493.  
  494.     d = Max( c, d );
  495.  
  496.     return d == 0.0 ? 0.0 : Abs( a - b ) / d;
  497.     }
  498.  
  499. static void
  500. RCheck( d, r )                /* check real number */
  501.     double    d;            /* real to be checked */
  502.     double    r;            /* expected value */
  503.     {
  504.     if ( RelDif( d, r ) > TOL )
  505.         {
  506.         errs = 1;
  507.         Printf( "value s.b. %g, was %g\n", r, d );
  508.         }
  509.     }
  510.  
  511. void
  512. GTest()
  513.     {
  514.     static struct
  515.         {
  516.         double    in;            /* input values */
  517.         double    exp;            /* expected output values */
  518.         }    tbl[] =         /* table of test values */
  519.         {
  520.         {    1.0,        1.0            },
  521.         {    2.0,        1.0            },
  522.         {    3.0,        2.0            },
  523.         {    4.0,        6.0            },
  524.         {    5.0,        24.0            },
  525.         {    6.0,        120.0            },
  526.         {    0.5,        1.7724538509        },
  527.         {    1.5,        0.8862269255        },
  528.         {    0.25,        3.6256099082        },
  529.         {    0.333333333333,    2.6789385347        },
  530.         {    0.666666666667,    1.3541179394        },
  531.         {    0.75,        1.2254167024        },
  532.         {    10.0,        362880.0        },
  533.         {    20.0,        1.2164510041e+17    },
  534.         };
  535.     register int    i;        /* indexes tbl[] test values */
  536.  
  537.     for ( i = 0; i < sizeof tbl / sizeof tbl[0]; ++i )
  538.         RCheck( Gamma( tbl[i].in ), tbl[i].exp );
  539.     }
  540.  
  541. void
  542. FTest()
  543.     {
  544.     static struct
  545.         {
  546.         int    in;            /* input values */
  547.         double    exp;            /* expected output values */
  548.         }    tbl[] =         /* table of test values */
  549.         {
  550.         {    0,        1.0            },
  551.         {    1,        1.0            },
  552.         {    2,        2.0            },
  553.         {    3,        6.0            },
  554.         {    4,        24.0            },
  555.         {    5,        120.0            },
  556.         {    6,        720.0            },
  557.         {    10,        3628800.0        },
  558.         {    20,        2.4329020082e+18    },
  559.         };
  560.     register int    i;        /* indexes tbl[] test values */
  561.  
  562.     for ( i = 0; i < sizeof tbl / sizeof tbl[0]; ++i )
  563.         RCheck( Factorial( tbl[i].in ), tbl[i].exp );
  564.     }
  565.  
  566. void
  567. BCTest()
  568.     {
  569.     static struct
  570.         {
  571.         int    n;            /* top parts of inputs */
  572.         int    k;            /* bottom parts of inputs */
  573.         double    exp;            /* expected output values */
  574.         }    tbl[] =         /* table of test values */
  575.         {
  576.         {    1,        0,        1.0        },
  577.         {    1,        1,        1.0        },
  578.         {    2,        0,        1.0        },
  579.         {    2,        1,        2.0        },
  580.         {    2,        2,        1.0        },
  581.         {    3,        0,        1.0        },
  582.         {    3,        1,        3.0        },
  583.         {    3,        2,        3.0        },
  584.         {    5,        3,        10.0        },
  585.         {    10,        4,        210.0        },
  586.         {    10,        5,        252.0        },
  587.         {    40,        6,        3838380.0    },
  588.         {    50,        20,        47129212243960.0},
  589.         };
  590.     register int    i;        /* indexes tbl[] test values */
  591.  
  592.     for ( i = 0; i < sizeof tbl / sizeof tbl[0]; ++i )
  593.         RCheck( BCoeff( tbl[i].n, tbl[i].k ), tbl[i].exp );
  594.     }
  595.  
  596. void
  597. ETest()
  598.     {
  599.     static struct
  600.         {
  601.         double    in;            /* input values */
  602.         double    exp;            /* expected output values */
  603.         }    tbl[] =         /* table of test values */
  604.         {
  605.         {    0.0,        0.0            },
  606.         {    0.1,        0.1124629160        },
  607.         {    0.2,        0.2227025892        },
  608.         {    0.5,        0.5204998778        },
  609.         {    0.8,        0.7421009647        },
  610.         {    1.0,        0.8427007929        },
  611.         {    1.5,        0.9661051465        },
  612.         {    2.0,        0.9953222650        },
  613.         };
  614.     register int    i;        /* indexes tbl[] test values */
  615.  
  616.     for ( i = 0; i < sizeof tbl / sizeof tbl[0]; ++i )
  617.         RCheck( Erf( tbl[i].in ), tbl[i].exp );
  618.     }
  619.  
  620. void
  621. QCTest()
  622.     {
  623.     static struct
  624.         {
  625.         double    chisq;            /* chi-square limit */
  626.         int    df;            /* degrees of freedom */
  627.         double    exp;            /* expected output values */
  628.         }    tbl[] =         /* table of test values */
  629.         {
  630.         {    0.001,        1,        0.97477        },
  631.         {    0.01,        1,        0.92034        },
  632.         {    0.01,        2,        0.99501        },
  633.         {    0.05,        1,        0.82306        },
  634.         {    0.1,        1,        0.75183        },
  635.         {    0.1,        2,        0.95123        },
  636.         {    1.0,        1,        0.31731        },
  637.         {    1.0,        2,        0.60653        },
  638.         {    1.0,        3,        0.80125        },
  639.         {    1.0,        4,        0.90980        },
  640.         {    1.0,        5,        0.96257        },
  641.         {    1.5,        2,        0.47237        },
  642.         {    2.0,        1,        0.15730        },
  643.         {    2.0,        3,        0.57241        },
  644.         {    2.0,        5,        0.84915        },
  645.         {    4.0,        6,        0.67668        },
  646.         {    5.0,        5,        0.41588        },
  647.         {    10.0,        7,        0.188573    },
  648.         {    10.0,        10,        0.44049        },
  649.         {    10.0,        20,        0.96817        },
  650.         {    20.0,        30,        0.91654        },
  651.         {    40.0,        30,        0.104864    },
  652.         {    8.26040,    20,        0.99        },
  653.         {    10.8508,    20,        0.95        },
  654.         {    12.4426,    20,        0.90        },
  655.         {    15.4518,    20,        0.75        },
  656.         {    19.3374,    20,        0.50        },
  657.         {    23.8277,    20,        0.25        },
  658.         {    28.4120,    20,        0.10        },
  659.         {    31.4104,    20,        0.05        },
  660.         {    37.5662,    20,        0.01        },
  661.         };
  662.     register int    i;        /* indexes tbl[] test values */
  663.  
  664.     for ( i = 0; i < sizeof tbl / sizeof tbl[0]; ++i )
  665.         RCheck( QChiSq( tbl[i].chisq, tbl[i].df ), tbl[i].exp );
  666.     }
  667.  
  668. /*ARGSUSED*/
  669. main( argc, argv )
  670.     int        argc;
  671.     char        *argv[];
  672.     {
  673.     GTest();
  674.     FTest();
  675.     BCTest();
  676.     ETest();
  677.     QCTest();
  678.  
  679.     return errs;
  680.     }
  681. END OF g_test.c
  682. echo 'gamma.3' 1>&2
  683. cat >'gamma.3' <<'END OF gamma.3'
  684. '\" e
  685. .TH GAMMA 3V VMB
  686. '\"    last edit:    88/09/09    D A Gwyn
  687. '\"    SCCS ID:    @(#)gamma.3    1.2 (edited for publication)
  688. '\"
  689. '\"    This source must be typeset via "eqn | troff -man"
  690. '\"
  691. .EQ
  692. delim @@
  693. .EN
  694. .SH NAME
  695. gamma \- gamma and related functions
  696. .SH SYNOPSIS
  697. .ds cW (CW\" change to B (without the paren) if you don't have a CW font
  698. \f\*(cW
  699. #include "gamma.h"
  700. .br
  701. /* \fPAll the following functions are declared in this header file.\f\*(cW */
  702. .sp
  703. double Gamma(double z);
  704. .br
  705. double LGamma(double z);
  706. .br
  707. double Factorial(int n);
  708. .br
  709. double LFactorial(int n);
  710. .br
  711. double BCoeff(int n, int k);
  712. .br
  713. double Beta(double z, double w);
  714. .br
  715. double PGamma(double a, double x);
  716. .br
  717. double QGamma(double a, double x);
  718. .br
  719. double Erf(double x);
  720. .br
  721. double Erfc(double x);
  722. .br
  723. double CPoisson(double x, int k);
  724. .br
  725. double PChiSq(double chisq, int nu);
  726. .br
  727. double QChiSq(double chisq, int nu);
  728. .ft
  729. .SH DESCRIPTION
  730. \f\*(cWGamma\fP
  731. returns the real gamma function value
  732. @GAMMA (z)~==~ int "" from 0 to inf t sup z-1 e sup -t dt@
  733. for @z>0@.
  734. .P
  735. \f\*(cWLGamma\fP
  736. returns the value
  737. @ln ~ GAMMA (z)@
  738. for @z>0@.
  739. .P
  740. \f\*(cWFactorial\fP
  741. returns the factorial function value
  742. @n! ~=~ GAMMA (n+1)@
  743. for @n>=0@.
  744. .P
  745. \f\*(cWLFactorial\fP
  746. returns the value
  747. @ln ~n!@\&
  748. for @n>=0@.
  749. .P
  750. \f\*(cWBCoeff\fP
  751. returns the binomial coefficient value
  752. @size +8 ( pile { n above k } size +8 ) ~==~n! over k!(n-k)!@\&
  753. for @0<=k<=n@.
  754. .P
  755. \f\*(cWBeta\fP
  756. returns the real beta function value
  757. @B(z,w)~==~ int "" from 0 to 1 t sup z-1 (1-t) sup w-1 dt@
  758. for @z>0@ and @w>0@.
  759. .P
  760. \f\*(cWPGamma\fP
  761. returns the incomplete gamma function value
  762. @P(a,x)~==~ 1 over { GAMMA (a) } int "" from 0 to x e sup -t t sup a-1 dt@
  763. for @a>0@ and @0<=x<= inf@.
  764. .P
  765. \f\*(cWQGamma\fP
  766. returns the value
  767. @1-P(a,x)@.
  768. .P
  769. \f\*(cWErf\fP
  770. returns the error function value
  771. @roman erf (x)~==~2 over sqrt pi int "" from 0 to x e sup {-t sup 2 } dt@.
  772. .P
  773. \f\*(cWErfc\fP
  774. returns the complementary error function value
  775. @1- roman erf (x)@.
  776. .P
  777. \f\*(cWCPoisson\fP
  778. returns the cumulative Poisson probability value
  779. @P sub x (<k)@,
  780. which is the probability that the number of Poisson random events occurring
  781. will be between 0 and @k-1@ inclusive,
  782. where the expected mean is @x@.
  783. .P
  784. \f\*(cWPChiSq\fP
  785. returns the value
  786. @P( chi sup 2 | nu )@,
  787. which is the probability that the observed chi-square for a correct model
  788. will be less than @chi sup 2@,
  789. where @nu@ is the number of degrees of freedom.
  790. .P
  791. \f\*(cWQChiSq\fP
  792. returns the value
  793. @1-P( chi sup 2 | nu )@,
  794. which is the probability that the observed chi-square
  795. will exceed @chi sup 2@ by chance even for a correct model,
  796. where @nu@ is the number of degrees of freedom.
  797. .P
  798. All these functions return
  799. .I approximations
  800. to the exact results,
  801. generally accurate to about seven significant digits.
  802. .SH CAVEATS
  803. Overflow can occur.
  804. .P
  805. Invalid parameters cause a diagnostic and termination of the process.
  806. .SH FILES
  807. .ta 3i
  808. gamma.h    header file
  809. .br
  810. gamma.o    object module
  811. .SH AUTHORS
  812. William H.\& Press,
  813. Brian P.\& Flannery,
  814. Saul A.\& Teukolsky,
  815. and
  816. William T.\& Vetterling
  817. (\fINumerical Recipes in C\fP, Chapter 6)
  818. .br
  819. Douglas A.\& Gwyn, BRL/VLD-VMB
  820. END OF gamma.3
  821. echo 'gamma.c' 1>&2
  822. cat >'gamma.c' <<'END OF gamma.c'
  823. /*
  824.     Gamma -- gamma and related functions
  825.  
  826.     last edit:    88/09/09    D A Gwyn
  827.  
  828.     SCCS ID:    @(#)gamma.c    1.1 (edited for publication)
  829.  
  830. Acknowledgement:
  831.     Code based on that found in "Numerical Methods in C".
  832. */
  833.  
  834. #include    <assert.h>
  835. #include    <math.h>
  836.  
  837. #include    "std.h"
  838.  
  839. double
  840. LGamma( x )
  841.     double            x;
  842.     {
  843.     static const double    cof[6] =
  844.         {
  845.         76.18009173,    -86.50532033,    24.01409822,
  846.         -1.231739516,    0.120858003e-2,    -0.536382e-5
  847.         };
  848.     double            tmp, ser;
  849.     register int        j;
  850.  
  851.     assert(x > 0.0);
  852.  
  853.     if ( --x < 0.0 )    /* use reflection formula for accuracy */
  854.         {
  855.         double    pix = PI * x;
  856.  
  857.         return log( pix / sin( pix ) ) - LGamma( 1.0 - x );
  858.         }
  859.  
  860.     tmp = x + 5.5;
  861.     tmp -= (x + 0.5) * log( tmp );
  862.  
  863.     ser = 1.0;
  864.  
  865.     for ( j = 0; j < 6; ++j )
  866.         ser += cof[j] / ++x;
  867.  
  868.     return -tmp + log( 2.50662827465 * ser );
  869.     }
  870.  
  871. double
  872. Gamma( x )
  873.     double    x;
  874.     {
  875.     return exp( LGamma( x ) );
  876.     }
  877.  
  878. double
  879. Factorial( n )
  880.     register int    n;
  881.     {
  882.     static double    a[33] =
  883.         {
  884.         1.0,    1.0,    2.0,    6.0,    24.0
  885.         };
  886.     static int    ntop = 4;
  887.  
  888.     assert(n >= 0);
  889.  
  890.     if ( n > 32 )
  891.         return Gamma( (double)n + 1.0 );
  892.  
  893.     while ( ntop < n )
  894.         {
  895.         register int    j = ntop++;
  896.  
  897.         a[ntop] = a[j] * (double)ntop;
  898.         }
  899.  
  900.     return a[n];
  901.     }
  902.  
  903. double
  904. LFactorial( n )
  905.     register int    n;
  906.     {
  907.     static double    a[101];        /* init 0.0 */
  908.  
  909.     assert(n >= 0);
  910.  
  911.     if ( n <= 1 )
  912.         return 0.0;
  913.  
  914.     if ( n <= 100 )
  915.         if ( a[n] > 0.0 )
  916.             return a[n];    /* table value already set up */
  917.         else
  918.             return a[n] = LGamma( (double)n + 1.0 );
  919.     else
  920.         return LGamma( (double)n + 1.0 );    /* beyond table */
  921.     }
  922.  
  923. double
  924. BCoeff( n, k )
  925.     register int    n, k;
  926.     {
  927.     assert(k >= 0);
  928.     assert(n >= k);
  929.  
  930.     return Round( exp( LFactorial( n )
  931.              - (LFactorial( k ) + LFactorial( n - k ))
  932.              )
  933.             );
  934.     }
  935.  
  936. double
  937. Beta( z, w )
  938.     double    z, w;
  939.     {
  940.     return exp( LGamma( z ) + LGamma( w ) - LGamma( z + w ) );
  941.     }
  942.  
  943. #define    ITMAX    100
  944. #define    EPS    3.0e-7
  945.  
  946. static double
  947. gser( a, x )
  948.     double        a, x;
  949.     {
  950.     double        ap, del, sum;
  951.     register int    n;
  952.  
  953.     assert(x >= 0.0);
  954.  
  955.     if ( x <= 0.0 )
  956.         return 0.0;
  957.  
  958.     del = sum = 1.0 / (ap = a);
  959.  
  960.     for ( n = 1; n <= ITMAX; ++n )
  961.         {
  962.         sum += del *= x / ++ap;
  963.  
  964.         if ( Abs( del ) < Abs( sum ) * EPS )
  965.             return sum * exp( -x + a * log( x ) - LGamma( a ) );
  966.         }
  967.  
  968.     assert(n <= ITMAX);
  969.     /*NOTREACHED*/
  970.     }
  971.  
  972. static double
  973. gcf( a, x )
  974.     double        a, x;
  975.     {
  976.     register int    n;
  977.     double        gold = 0.0, fac = 1.0, b1 = 1.0,
  978.             b0 = 0.0, a0 = 1.0, a1 = x;
  979.  
  980.     for ( n = 1; n <= ITMAX; ++n )
  981.         {
  982.         double    anf;
  983.         double    an = (double)n;
  984.         double    ana = an - a;
  985.  
  986.         a0 = (a1 + a0 * ana) * fac;
  987.         b0 = (b1 + b0 * ana) * fac;
  988.         anf = an * fac;
  989.         b1 = x * b0 + anf * b1;
  990.         a1 = x * a0 + anf * a1;
  991.  
  992.         if ( a1 != 0.0 )
  993.             {        /* renormalize */
  994.             double    g = b1 * (fac = 1.0 / a1);
  995.  
  996.             gold = g - gold;
  997.  
  998.             if ( Abs( gold ) < EPS * Abs( g ) )
  999.                 return exp( -x + a * log( x ) - LGamma( a ) )
  1000.                     * g;
  1001.  
  1002.             gold = g;
  1003.             }
  1004.         }
  1005.  
  1006.     assert(n <= ITMAX);
  1007.     /*NOTREACHED*/
  1008.     }
  1009.  
  1010. double
  1011. PGamma( a, x )
  1012.     double    a, x;
  1013.     {
  1014.     assert(x >= 0.0);
  1015.     assert(a > 0.0);
  1016.  
  1017.     return x < a + 1.0 ? gser( a, x ) : 1.0 - gcf( a, x );
  1018.     }
  1019.  
  1020. double
  1021. QGamma( a, x )
  1022.     double    a, x;
  1023.     {
  1024.     assert(x >= 0.0);
  1025.     assert(a > 0.0);
  1026.  
  1027.     return x < a + 1.0 ? 1.0 - gser( a, x ) : gcf( a, x );
  1028.     }
  1029.  
  1030. double
  1031. Erf( x )
  1032.     double    x;
  1033.     {
  1034.     return x < 0.0 ? -PGamma( 0.5, x * x ) : PGamma( 0.5, x * x );
  1035.     }
  1036.  
  1037. double
  1038. Erfc( x )
  1039.     double    x;
  1040.     {
  1041.     return x < 0.0 ? 1.0 + PGamma( 0.5, x * x ) : QGamma( 0.5, x * x );
  1042.     }
  1043.  
  1044. double
  1045. CPoisson( x, k )
  1046.     double    x;
  1047.     int    k;
  1048.     {
  1049.     return QGamma( (double)k, x );
  1050.     }
  1051.  
  1052. double
  1053. PChiSq( chisq, df )
  1054.     double    chisq;
  1055.     int    df;
  1056.     {
  1057.     return PGamma( (double)df / 2.0, chisq / 2.0 );
  1058.     }
  1059.  
  1060. double
  1061. QChiSq( chisq, df )
  1062.     double    chisq;
  1063.     int    df;
  1064.     {
  1065.     return QGamma( (double)df / 2.0, chisq / 2.0 );
  1066.     }
  1067. END OF gamma.c
  1068. echo 'gamma.h' 1>&2
  1069. cat >'gamma.h' <<'END OF gamma.h'
  1070. /*
  1071.     <gamma.h> -- definitions for gamma-function routines
  1072.  
  1073.     last edit:    88/09/09    D A Gwyn
  1074.  
  1075.     SCCS ID:    @(#)gamma.h    1.1 (edited for publication)
  1076. */
  1077.  
  1078. /* library routines: */
  1079.  
  1080. #ifdef __STDC__
  1081. extern double    Gamma( double z ), LGamma( double z ),
  1082.         Factorial( int n ), LFactorial( int n ),
  1083.         BCoeff( int n, int k ),
  1084.         Beta( double z, double w ),
  1085.         PGamma( double a, double x ), QGamma( double a, double x ),
  1086.         Erf( double x ), Erfc( double x ),
  1087.         CPoisson( double x, int k ),
  1088.         PChiSq( double chisq, int nu ), QChiSq( double chisq, int nu );
  1089. #else
  1090. extern double    Gamma(), LGamma(), Factorial(), LFactorial(),
  1091.         BCoeff(), Beta(), PGamma(), QGamma(), Erf(), Erfc(),
  1092.         CPoisson(), PChiSq(), QChiSq();
  1093. #endif
  1094. END OF gamma.h
  1095. echo 'info.c' 1>&2
  1096. cat >'info.c' <<'END OF info.c'
  1097. /*
  1098.     InfoTbl -- Kullback's information measure for a 2-way contingency table
  1099.  
  1100.     last edit:    88/09/19    D A Gwyn
  1101.  
  1102.     SCCS ID:    @(#)info.c    1.1 (edited for publication)
  1103.  
  1104.     Special return values:
  1105.         -1.0    entire table consisted of 0 entries
  1106.         -2.0    invalid table entry (frequency less than 0)
  1107.         -3.0    invalid table dimensions (r or c less than 2)
  1108.         -4.0    unable to allocate enough working storage
  1109. */
  1110.  
  1111. #include    <math.h>        /* for log() */
  1112. #if __STDC__
  1113. #include    <stdlib.h>        /* malloc, free */
  1114.  
  1115. #include    "std.h"
  1116. #else
  1117. #include    "std.h"
  1118.  
  1119. extern pointer    malloc();
  1120. extern void    free();
  1121. #endif
  1122.  
  1123. double
  1124. InfoTbl( r, c, f, pdf )
  1125.     int        r;        /* # rows in table */
  1126.     int        c;        /* # columns in table */
  1127.     const long    *f;        /* -> r*c frequency tallies */
  1128.     int        *pdf;        /* -> return # d.f. for chi-square */
  1129.     {
  1130. #define    x(i,j)    f[(i)*c+(j)]        /* convenient way to access freqs */
  1131.     register int    i;        /* row index */
  1132.     register int    j;        /* column index */
  1133.     double        N;        /* (double)n */
  1134.     double        info;        /* accumulates information measure */
  1135.     double        *xi;        /* row sums */
  1136.     double        *xj;        /* col sums */
  1137.     int        rdf = r - 1;    /* row degrees of freedom */
  1138.     int        cdf = c - 1;    /* column degrees of freedom */
  1139.  
  1140.     if ( rdf <= 0 || cdf <= 0 )
  1141.         {
  1142.         info = -3.0;
  1143.         goto ret3;
  1144.         }
  1145.  
  1146.     *pdf = rdf * cdf;        /* total degrees of freedom */
  1147.  
  1148.     if ( (xi = (double *)malloc( r * sizeof(double) )) == NULL )
  1149.         {
  1150.         info = -4.0;
  1151.         goto ret3;
  1152.         }
  1153.  
  1154.     if ( (xj = (double *)malloc( c * sizeof(double) )) == NULL )
  1155.         {
  1156.         info = -4.0;
  1157.         goto ret2;
  1158.         }
  1159.  
  1160.     /* compute row sums and total */
  1161.  
  1162.     N = 0.0;
  1163.  
  1164.     for ( i = 0; i < r; ++i )
  1165.         {
  1166.         double    sum = 0.0;    /* accumulator */
  1167.  
  1168.         for ( j = 0; j < c; ++j )
  1169.             {
  1170.             register long    k = x(i,j);
  1171.  
  1172.             if ( k < 0L )
  1173.                 {
  1174.                 info = -2.0;
  1175.                 goto ret1;
  1176.                 }
  1177.  
  1178.             sum += (double)k;
  1179.             }
  1180.  
  1181.         N += xi[i] = sum;
  1182.         }
  1183.  
  1184.     if ( N <= 0.0 )
  1185.         {
  1186.         info = -1.0;
  1187.         goto ret1;
  1188.         }
  1189.  
  1190.     /* compute column sums */
  1191.  
  1192.     for ( j = 0; j < c; ++j )
  1193.         {
  1194.         double    sum = 0.0;    /* accumulator */
  1195.  
  1196.         for ( i = 0; i < r; ++i )
  1197.             sum += (double)x(i,j);
  1198.  
  1199.         xj[j] = sum;
  1200.         }
  1201.  
  1202.     /* compute information measure (four parts) */
  1203.  
  1204.     info = N * log( N );                    /* part 1 */
  1205.  
  1206.     for ( i = 0; i < r; ++i )
  1207.         {
  1208.         double    pi = xi[i];    /* row sum */
  1209.  
  1210.         if ( pi > 0.0 )
  1211.             info -= pi * log( pi );            /* part 2 */
  1212.  
  1213.         for ( j = 0; j < c; ++j )
  1214.             {
  1215.             double    pij = (double)x(i,j);
  1216.  
  1217.             if ( pij > 0.0 )
  1218.                 info += pij * log( pij );    /* part 3 */
  1219.             }
  1220.         }
  1221.  
  1222.     for ( j = 0; j < c; ++j )
  1223.         {
  1224.         double    pj = xj[j];    /* column sum */
  1225.  
  1226.         if ( pj > 0.0 )
  1227.             info -= pj * log( pj );            /* part 4 */
  1228.         }
  1229.  
  1230.     info *= 2.0;            /* for comparability with chi-square */
  1231.  
  1232.     ret1:
  1233.     free( (pointer)xj );
  1234.     ret2:
  1235.     free( (pointer)xi );
  1236.     ret3:
  1237.     return info;
  1238.     }
  1239. END OF info.c
  1240. echo 'std.h' 1>&2
  1241. cat >'std.h' <<'END OF std.h'
  1242. /*
  1243.     std.h -- Douglas A. Gwyn's standard C programming definitions
  1244.  
  1245.     Prerequisite:    <math.h> (if you invoke Round())
  1246.  
  1247.     last edit:    89/06/18    D A Gwyn
  1248.  
  1249.     SCCS ID:    @(#)std.h    1.30
  1250.  
  1251.     The master source file is to be modified only by Douglas A. Gwyn
  1252.     <Gwyn@BRL.MIL>.  When installing a VLD/VMB software distribution,
  1253.     this file may need to be tailored slightly to fit the target system.
  1254.     Usually this just involves enabling some of the "kludges for deficient
  1255.     C implementations" at the end of this file.
  1256. */
  1257.  
  1258. #ifndef    VLD_STD_H_
  1259. #define    VLD_STD_H_            /* once-only latch */
  1260.  
  1261. /* Extended data types */
  1262.  
  1263. typedef int    bool;            /* Boolean data */
  1264. #define     false    0
  1265. #define     true    1
  1266.  
  1267. /* ANSI C definitions */
  1268.  
  1269. /* Defense against some silly systems defining __STDC__ to random things. */
  1270. #ifdef STD_C
  1271. #undef STD_C
  1272. #endif
  1273. #ifdef __STDC__
  1274. #if __STDC__ > 0
  1275. #define    STD_C    __STDC__        /* use this instead of __STDC__ */
  1276. #endif
  1277. #endif
  1278.  
  1279. #ifdef STD_C
  1280. typedef void    *pointer;        /* generic pointer */
  1281. #else
  1282. typedef char    *pointer;        /* generic pointer */
  1283. #define    const        /* nothing */    /* ANSI C type qualifier */
  1284. /* There really is no substitute for the following, but these might work: */
  1285. #define    signed        /* nothing */    /* ANSI C type specifier */
  1286. #define    volatile    /* nothing */    /* ANSI C type qualifier */
  1287. #endif
  1288.  
  1289. #ifndef EXIT_SUCCESS
  1290. #define    EXIT_SUCCESS    0
  1291. #endif
  1292.  
  1293. #ifndef EXIT_FAILURE
  1294. #define    EXIT_FAILURE    1
  1295. #endif
  1296.  
  1297. #ifndef NULL
  1298. #define NULL    0            /* null pointer constant, all types */
  1299. #endif
  1300.  
  1301. /* Universal constants */
  1302.  
  1303. #define DEGRAD    57.2957795130823208767981548141051703324054724665642
  1304.                     /* degrees per radian */
  1305. #define    E    2.71828182845904523536028747135266249775724709369996
  1306.                     /* base of natural logs */
  1307. #define    GAMMA    0.57721566490153286061
  1308.                     /* Euler's constant */
  1309. #define LOG10E    0.43429448190325182765112891891660508229439700580367
  1310.                     /* log of e to the base 10 */
  1311. #define PHI    1.618033988749894848204586834365638117720309180
  1312.                     /* golden ratio */
  1313. #define PI    3.14159265358979323846264338327950288419716939937511
  1314.                     /* ratio of circumf. to diam. */
  1315.  
  1316. /* Useful macros */
  1317.  
  1318. /*
  1319.     The comment "UNSAFE" means that the macro argument(s) may be evaluated
  1320.     more than once, so the programmer must realize that the macro doesn't
  1321.     quite act like a genuine function.  This matters only when evaluating
  1322.     an argument produces "side effects".
  1323. */
  1324.  
  1325. /* arbitrary numerical arguments and value: */
  1326. #define Abs( x )    ((x) < 0 ? -(x) : (x))            /* UNSAFE */
  1327. #define Max( a, b )    ((a) > (b) ? (a) : (b))            /* UNSAFE */
  1328. #define Min( a, b )    ((a) < (b) ? (a) : (b))            /* UNSAFE */
  1329.  
  1330. /* floating-point arguments and value: */
  1331. #define Round( d )    floor( (d) + 0.5 )        /* requires <math.h> */
  1332.  
  1333. /* arbitrary numerical arguments, integer value: */
  1334. #define    Sgn( x )    ((x) == 0 ? 0 : (x) > 0 ? 1 : -1)    /* UNSAFE */
  1335.  
  1336. /* string arguments, boolean value: */
  1337. #ifdef gould    /* UTX-32 2.0 compiler has problems with "..."[] */
  1338. #define    StrEq( a, b )    (strcmp( a, b ) == 0)            /* UNSAFE */
  1339. #else
  1340. #define    StrEq( a, b )    (*(a) == *(b) && strcmp( a, b ) == 0)    /* UNSAFE */
  1341. #endif
  1342.  
  1343. /* array argument, integer value: */
  1344. #define    Elements( a )    (sizeof a / sizeof a[0])
  1345.  
  1346. /* integer (or character) arguments and value: */
  1347. #define fromhostc( c )    (c)        /* map host char set to ASCII */
  1348. #define tohostc( c )    (c)        /* map ASCII to host char set */
  1349. #define tonumber( c )    ((c) - '0')    /* convt digit char to number */
  1350. #define todigit( n )    ((n) + '0')    /* convt digit number to char */
  1351.  
  1352. /* Kludges for deficient C implementations */
  1353.  
  1354. /*#define    strchr    index        /* 7th Edition UNIX, 4.2BSD */
  1355. /*#define    strrchr    rindex        /* 7th Edition UNIX, 4.2BSD */
  1356. /*#define    void    int        /* K&R Appendix A followers */
  1357.  
  1358. #endif    /* VLD_STD_H_ */
  1359. END OF std.h
  1360. echo 'tot_info.1' 1>&2
  1361. cat >'tot_info.1' <<'END OF tot_info.1'
  1362. '\" e
  1363. .TH TOT_INFO 1V VMB
  1364. '\"    last edit:    89/02/07    D A Gwyn
  1365. '\"    SCCS ID:    @(#)tot_info.1    1.2 (edited for publication)
  1366. '\"
  1367. '\"    This source must be typeset via "eqn | troff -man"
  1368. '\"
  1369. .EQ
  1370. delim @@
  1371. .EN
  1372. .SH NAME
  1373. tot_info \- total information for multiple 2-way contingency tables
  1374. .SH SYNOPSIS
  1375. .ds cW (CW\" change to I (without the paren) if you don't have a CW font
  1376. .ds cB (CB\" change to B (without the paren) if you don't have a CB font
  1377. \f\*(cBtot_info\fP
  1378. .SH DESCRIPTION
  1379. \f\*(cBtot_info\fP
  1380. reads multiple contingency-table data sets from the standard input
  1381. and prints Kullback's information statistic for each set,
  1382. as well as for the aggregate over all sets,
  1383. on the standard output.
  1384. (See \fIchisq\^\fP(3V) for a discussion of this statistic.)
  1385. .SH "INPUT FORMAT"
  1386. Input consists of one or more data sets,
  1387. each constituting a 2-way contingency table
  1388. (not necessarily all of the same size).
  1389. A data set may be preceded by any number of blank or comment lines;
  1390. a comment line has a
  1391. \f\*(cB#\fP
  1392. character as the first non-whitespace character on the line.
  1393. Following the optional comment lines is a header line
  1394. containing the row and column dimensions of the contingency table
  1395. (in that order),
  1396. separated by white space.
  1397. Finally, the contents of the contingency table (frequency counts)
  1398. must be given in ``row major'' order.
  1399. The table may be freely formatted with white-space separators;
  1400. a row need not be on a single line.
  1401. .SH "OUTPUT FORMAT"
  1402. Input comments are copied to the output as they are encountered.
  1403. Otherwise the output consists solely of an information line for
  1404. each data set (or a diagnostic if the data is invalid) and a final
  1405. information line for the aggregate over all data sets (preceded by
  1406. a blank line).
  1407. An information line exhibits twice the value of Kullback's
  1408. @"\v'-0.2'^\v'0.2'\h'-\w;^;u'I\^" ( H sub 1 : H sub 2 )@ statistic,
  1409. the corresponding number of degrees of freedom,
  1410. and the probability that the statistic
  1411. would be as large as was actually observed
  1412. if the row and column categorizations really were independent.
  1413. .P
  1414. The aggregate statistic is valid if the data sets
  1415. represent independent tests of the same hypothesis.
  1416. .SH DIAGNOSTICS
  1417. The diagnostic messages are intended to be self-explanatory.
  1418. .SH "EXIT STATUS"
  1419. \f\*(cBtot_info\fP
  1420. returns a zero exit status if and only if no problems were encountered.
  1421. .SH EXAMPLE
  1422. .RS
  1423. \f\*(cW$ \fP\f\*(cBtot_info
  1424. .br
  1425. # MilCrypI.60 biliteral cryptogram (trial pairings against G)
  1426. .br
  1427. # G vs B
  1428. .br
  1429. 2 10
  1430. .br
  1431. 2 2 2 0 3 0 0 1 0 1
  1432. .br
  1433. 3 1 1 1 1 2 2 1 2 1
  1434. .br
  1435. # G vs D
  1436. .br
  1437. 2 10
  1438. .br
  1439. 2 2 2 0 3 0 0 1 0 1
  1440. .br
  1441. 4 1 4 4 1 1 1 3 4 2
  1442. .br
  1443. # G vs J
  1444. .br
  1445. 2 10
  1446. .br
  1447. 2 2 2 0 3 0 0 1 0 1
  1448. .br
  1449. 1 1 1 1 1 1 2 1 1 1
  1450. .br
  1451. # G vs L
  1452. .br
  1453. 2 10
  1454. .br
  1455. 2 2 2 0 3 0 0 1 0 1
  1456. .br
  1457. 1 4 0 4 3 4 5 3 3 4
  1458. .br
  1459. # G vs N
  1460. .br
  1461. 2 10
  1462. .br
  1463. 2 2 2 0 3 0 0 1 0 1
  1464. .br
  1465. 4 1 4 3 1 1 1 2 3 3
  1466. .br
  1467. # G vs Q
  1468. .br
  1469. 2 10
  1470. .br
  1471. 2 2 2 0 3 0 0 1 0 1
  1472. .br
  1473. 0 2 0 2 1 1 1 0 1 0
  1474. .br
  1475. # G vs S (correct pairing)
  1476. .br
  1477. 2 10
  1478. .br
  1479. 2 2 2 0 3 0 0 1 0 1
  1480. .br
  1481. 1 2 2 0 2 1 0 0 0 1
  1482. .br
  1483. # G vs V
  1484. .br
  1485. 2 10
  1486. .br
  1487. 2 2 2 0 3 0 0 1 0 1
  1488. .br
  1489. 1 4 1 3 4 4 4 3 4 3
  1490. .br
  1491. # G vs X
  1492. .br
  1493. 2 10
  1494. .br
  1495. 2 2 2 0 3 0 0 1 0 1
  1496. .br
  1497. 0 1 0 1 2 1 1 0 2 0
  1498. .br
  1499. # Since most pairings are incorrect,
  1500. .br
  1501. # the aggregate probability is small.
  1502. .br
  1503. ^D\fP\f\*(cW
  1504. .br
  1505. # MilCrypI.60 biliteral cryptogram (trial pairings against G)
  1506. .br
  1507. # G vs B
  1508. .br
  1509. 2info = 11.01    df =  9    q =  0.2748
  1510. .br
  1511. # G vs D
  1512. .br
  1513. 2info = 12.40    df =  9    q =  0.1915
  1514. .br
  1515. # G vs J
  1516. .br
  1517. 2info =  9.00    df =  9    q =  0.4375
  1518. .br
  1519. # G vs L
  1520. .br
  1521. 2info = 19.03    df =  9    q =  0.0250
  1522. .br
  1523. # G vs N
  1524. .br
  1525. 2info = 10.89    df =  9    q =  0.2830
  1526. .br
  1527. # G vs Q
  1528. .br
  1529. 2info = 15.82    df =  9    q =  0.0707
  1530. .br
  1531. # G vs S (correct pairing)
  1532. .br
  1533. 2info =  3.11    df =  9    q =  0.9596
  1534. .br
  1535. # G vs V
  1536. .br
  1537. 2info = 14.47    df =  9    q =  0.1066
  1538. .br
  1539. # G vs X
  1540. .br
  1541. 2info = 15.31    df =  9    q =  0.0826
  1542. .br
  1543. # Since most pairings are incorrect,
  1544. .br
  1545. # the aggregate probability is small.
  1546. .br
  1547. .sp
  1548. total 2info = 111.05    df = 81    q =  0.0150
  1549. \fP
  1550. .RE
  1551. .SH REFERENCES
  1552. Solomon Kullback, \fIInformation Theory and Statistics\fP (Dover, 1968).
  1553. .br
  1554. William F.\& Friedman and Lambros D.\& Callimahos,
  1555. \fIMilitary Cryptanalytics, Part I \(em Volume 1\fP
  1556. (reprinted by Aegean Park Press, 1985).
  1557. .SH "SEE ALSO"
  1558. chisq(3V).
  1559. .SH AUTHOR
  1560. Douglas A.\& Gwyn, BRL/VLD-VMB
  1561. END OF tot_info.1
  1562. echo 'tot_info.c' 1>&2
  1563. cat >'tot_info.c' <<'END OF tot_info.c'
  1564. /*
  1565.     tot_info -- combine information statistics for multiple tables
  1566.  
  1567.     last edit:    89/02/06    D A Gwyn
  1568.  
  1569.     SCCS ID:    @(#)tot_info.c    1.1 (edited for publication)
  1570. */
  1571.  
  1572. #include    <ctype.h>
  1573. #include    <stdio.h>
  1574.  
  1575. #include    "std.h"
  1576.  
  1577. #include    "chisq.h"
  1578. #include    "gamma.h"        /* for QChiSq() */
  1579.  
  1580. #ifndef MAXLINE
  1581. #define    MAXLINE    256
  1582. #endif
  1583.  
  1584. #ifndef MAXTBL
  1585. #define    MAXTBL    1000
  1586. #endif
  1587.  
  1588. static char    line[MAXLINE];        /* row/column header input line */
  1589. static long    f[MAXTBL];        /* frequency tallies */
  1590. static int    r;            /* # of rows */
  1591. static int    c;            /* # of columns */
  1592.  
  1593. #define    x(i,j)    f[(i)*c+(j)]        /* convenient way to access freqs */
  1594.  
  1595. #define    COMMENT    '#'            /* comment character */
  1596.  
  1597. /*ARGSUSED*/
  1598. int
  1599. main( argc, argv )
  1600.     int        argc;
  1601.     char        *argv[];
  1602.     {
  1603.     register char    *p;        /* input line scan location */
  1604.     register int    i;        /* row index */
  1605.     register int    j;        /* column index */
  1606.     double        info;        /* computed information measure */
  1607.     int        infodf;        /* degrees of freedom for information */
  1608.     double        totinfo = 0.0;    /* accumulated information */
  1609.     int        totdf = 0;    /* accumulated degrees of freedom */
  1610.  
  1611.     while ( fgets( line, MAXLINE, stdin ) != NULL )    /* start new table */
  1612.         {
  1613.         for ( p = line; *p != '\0' && isspace( (int)*p ); ++p )
  1614.             ;
  1615.  
  1616.         if ( *p == '\0' )
  1617.             continue;    /* skip blank line */
  1618.  
  1619.         if ( *p == COMMENT )
  1620.             {        /* copy comment through */
  1621.             (void)fputs( line, stdout );
  1622.             continue;
  1623.             }
  1624.  
  1625.         if ( sscanf( p, "%d %d\n", &r, &c ) != 2 )
  1626.             {
  1627.             (void)fputs( "* invalid row/column line *\n", stderr );
  1628.             return EXIT_FAILURE;
  1629.             }
  1630.  
  1631.         if ( r * c > MAXTBL )
  1632.             {
  1633.             (void)fputs( "* table too large *\n", stderr );
  1634.             return EXIT_FAILURE;
  1635.             }
  1636.  
  1637.         /* input tallies */
  1638.  
  1639.         for ( i = 0; i < r; ++i )
  1640.             for ( j = 0; j < c; ++j )
  1641.                 if ( scanf( " %ld", &x(i,j) ) != 1 )
  1642.                     {
  1643.                     (void)fputs( "* EOF in table *\n",
  1644.                              stderr
  1645.                            );
  1646.                     return EXIT_FAILURE;
  1647.                     }
  1648.  
  1649.         /* compute statistic */
  1650.  
  1651.         info = InfoTbl( r, c, f, &infodf );
  1652.  
  1653.         /* print results */
  1654.  
  1655.         if ( info >= 0.0 )
  1656.             {
  1657.             (void)printf( "2info = %5.2f\tdf = %2d\tq = %7.4f\n",
  1658.                       info, infodf,
  1659.                       QChiSq( info, infodf )
  1660.                     );
  1661.             totinfo += info;
  1662.             totdf += infodf;
  1663.             }
  1664.         else
  1665.             (void)fputs( info < -3.5 ? "out of memory\n"
  1666.                    : info < -2.5 ? "table too small\n"
  1667.                    : info < -1.5 ? "negative freq\n"
  1668.                    : "table all zeros\n",
  1669.                      stdout
  1670.                    );
  1671.         }
  1672.  
  1673.     if ( totdf <= 0 )
  1674.         {
  1675.         (void)fputs( "\n*** no information accumulated ***\n", stdout );
  1676.         return EXIT_FAILURE;
  1677.         }
  1678.  
  1679.     (void)printf( "\ntotal 2info = %5.2f\tdf = %2d\tq = %7.4f\n",
  1680.               totinfo, totdf,
  1681.               QChiSq( totinfo, totdf )
  1682.             );
  1683.     return EXIT_SUCCESS;
  1684.     }
  1685. END OF tot_info.c
  1686. echo 'tot_info.exp' 1>&2
  1687. cat >'tot_info.exp' <<'END OF tot_info.exp'
  1688. # MilCrypI.60 biliteral cryptogram (trial pairings against G)
  1689. # G vs B
  1690. 2info = 11.01    df =  9    q =  0.2748
  1691. # G vs D
  1692. 2info = 12.40    df =  9    q =  0.1915
  1693. # G vs J
  1694. 2info =  9.00    df =  9    q =  0.4375
  1695. # G vs L
  1696. 2info = 19.03    df =  9    q =  0.0250
  1697. # G vs N
  1698. 2info = 10.89    df =  9    q =  0.2830
  1699. # G vs Q
  1700. 2info = 15.82    df =  9    q =  0.0707
  1701. # G vs S (correct pairing)
  1702. 2info =  3.11    df =  9    q =  0.9596
  1703. # G vs V
  1704. 2info = 14.47    df =  9    q =  0.1066
  1705. # G vs X
  1706. 2info = 15.31    df =  9    q =  0.0826
  1707. # Since most pairings are incorrect,
  1708. # the aggregate probability is small.
  1709.  
  1710. total 2info = 111.05    df = 81    q =  0.0150
  1711. END OF tot_info.exp
  1712. echo 'tot_info.in' 1>&2
  1713. cat >'tot_info.in' <<'END OF tot_info.in'
  1714. # MilCrypI.60 biliteral cryptogram (trial pairings against G)
  1715. # G vs B
  1716. 2 10
  1717. 2 2 2 0 3 0 0 1 0 1
  1718. 3 1 1 1 1 2 2 1 2 1
  1719. # G vs D
  1720. 2 10
  1721. 2 2 2 0 3 0 0 1 0 1
  1722. 4 1 4 4 1 1 1 3 4 2
  1723. # G vs J
  1724. 2 10
  1725. 2 2 2 0 3 0 0 1 0 1
  1726. 1 1 1 1 1 1 2 1 1 1
  1727. # G vs L
  1728. 2 10
  1729. 2 2 2 0 3 0 0 1 0 1
  1730. 1 4 0 4 3 4 5 3 3 4
  1731. # G vs N
  1732. 2 10
  1733. 2 2 2 0 3 0 0 1 0 1
  1734. 4 1 4 3 1 1 1 2 3 3
  1735. # G vs Q
  1736. 2 10
  1737. 2 2 2 0 3 0 0 1 0 1
  1738. 0 2 0 2 1 1 1 0 1 0
  1739. # G vs S (correct pairing)
  1740. 2 10
  1741. 2 2 2 0 3 0 0 1 0 1
  1742. 1 2 2 0 2 1 0 0 0 1
  1743. # G vs V
  1744. 2 10
  1745. 2 2 2 0 3 0 0 1 0 1
  1746. 1 4 1 3 4 4 4 3 4 3
  1747. # G vs X
  1748. 2 10
  1749. 2 2 2 0 3 0 0 1 0 1
  1750. 0 1 0 1 2 1 1 0 2 0
  1751. # Since most pairings are incorrect,
  1752. # the aggregate probability is small.
  1753. END OF tot_info.in
  1754.  
  1755.