home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / gnu / graphics-0.17 / dist-stat / cor.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-11-13  |  9.2 KB  |  439 lines
  1. /*
  2.  * $Header: /files1/home/toy/src/stat/RCS/cor.c,v 1.5 90/09/11 16:45:19 toy Exp $
  3.  * NAME
  4.  *    cor - ordinary correlation coefficient
  5.  *
  6.  * SYNOPSIS
  7.  *    cor [-H] [-P prec] [-c cols] [-F base] [vectors...]
  8.  *    lreg [-H] [-P prec] [-c cols] [-F base] [-ios] [vectors...]
  9.  *
  10.  * DESCRIPTION
  11.  *    Output for cor is the ordinary correlation coefficient between
  12.  *    the base vector and the given vectors.  If the base is not given,
  13.  *    it is assumed to come from stdin.
  14.  *
  15.  *    The correlation is:
  16.  *
  17.  *        E[(x - mx)(y - my)]
  18.  *        -------------------
  19.  *        sqrt(Var x * Var y)
  20.  *
  21.  *    where mx, and my are the means.
  22.  *
  23.  *    For lreg, the output is the slope and intercept from a least
  24.  *    squares linear regression of each vector on a base vector.
  25.  *    If the base is not given, it is assumed to be ascending non-
  26.  *    negative integers from 0.
  27.  *
  28.  * BUGS
  29.  *    The entire vector is read into memory to compute the variance.
  30.  *    This is done for accuracy (but is it really necessary?).
  31.  *
  32.  * HISTORY
  33.  * $Log:    cor.c,v $
  34.  * Revision 1.5  90/09/11  16:45:19  toy
  35.  * If vector is shorter than the base vector, the message should
  36.  * include the file name and the number of elements read.
  37.  *
  38.  * Also recognize "," as an option character that does nothing.
  39.  *
  40.  * Revision 1.4  90/09/05  13:53:54  toy
  41.  * Implemented lreg routine here.
  42.  *
  43.  * Revision 1.3  90/09/04  16:04:03  toy
  44.  * Use print_help_strings to print out the help strings.
  45.  * Get the default precision from the environment variable STAT_PREC.
  46.  *
  47.  * Revision 1.2  90/09/04  14:53:57  toy
  48.  * The correlation coefficient should be computed with
  49.  * the appropriate means subtracted out.
  50.  *
  51.  * Revision 1.1  90/09/02  19:01:47  toy
  52.  * Initial revision
  53.  *
  54.  */
  55.  
  56. #include <stdio.h>
  57. #include <math.h>
  58. #include <string.h>
  59. #include <assert.h>
  60.  
  61. #if    defined(__STDC__) || defined(__GNUC__)
  62. #include <stddef.h>
  63. #include <stdlib.h>
  64. #else
  65. #include <malloc.h>
  66. #endif
  67.  
  68. #include "gps.h"
  69.  
  70. #define    DEF_PROGNAME    "cor|lreg"
  71. #define    COR_NAME    "cor"
  72. #define    LREG_NAME    "lreg"
  73.  
  74. static const char RCSID[] = "@(#) $Id: cor.c,v 1.5 90/09/11 16:45:19 toy Exp $";
  75.  
  76. extern int getopt ();
  77. extern int optind;
  78. extern char *optarg;
  79.  
  80. const char *progname;
  81. const char *vector_name;
  82.  
  83. int intercept_output = TRUE;    /* lreg:  output intercept     */
  84. int slope_output = TRUE;    /* lreg:  output slope         */
  85. int opt_form = FALSE;        /* lreg:  output in option form     */
  86.  
  87. #define    COR_OPT_STRING    ",HP:c:F:"
  88. #define    LREG_OPT_STRING    ",HP:c:F:ios"
  89.  
  90. #if    defined(__STDC__) || defined(__GNUC__)
  91. void
  92. help (void)
  93. #else
  94. void
  95. help ()
  96. #endif
  97. {
  98.   static const char *help_strings[] =
  99.   {
  100.     "\t-H\tThis help\n",
  101.     "\t-P p\tSet number of significant digits to p.  If not specified\n",
  102.     "\t\tuse value of STAT_PREC from environment, if available.  Else\n",
  103.     "\t\tuse default.\n",
  104.     "\t-c c\tSet number of output elements per line to n\n",
  105.     "\t-F F\tBase (x) vector is read from file F.  If not specified, the\n",
  106.     "\t\tdata is read from stdin.\n",
  107.     NULL
  108.   };
  109.   static const char *help_cor[] =
  110.   {
  111.     "\tvec..\tData vectors.  If not given, use stdin.\n",
  112.     "Compute and print the correlation coefficient between the base vector and\n",
  113.     "the given data vectors.  (The coefficient is computed with the mean removed.)\n",
  114.     NULL
  115.   };
  116.   static const char *help_lreg[] =
  117.   {
  118.     "\t-i\tOutput intercept only\n",
  119.     "\t-s\tOutput slope only\n",
  120.     "\t-o\tOutput result in option form (for siline)\n",
  121.     "Compute the slope and intercept of a least-squares linear regression\n",
  122.     "between the base (x) vector and the give data (y) vectors.\n",
  123.     "The intercept and then slope are printed.\n",
  124.     NULL
  125.   };
  126.  
  127.   (void) fprintf (stderr, "%s\n", RCSID);
  128.  
  129.   if (strcmp (progname, COR_NAME) == 0)
  130.     {
  131.       (void) fprintf (stderr, "Usage:  %s [-H] [-P prec] [-c cols] [-F base] [vector ...]\n", progname);
  132.       print_help_strings (help_strings);
  133.       print_help_strings (help_cor);
  134.     }
  135.   else
  136.     {
  137.       (void) fprintf (stderr, "Usage:  %s [-Hios] [-P prec] [-c cols] [-F base] [vector ...]\n", progname);
  138.       print_help_strings (help_strings);
  139.       print_help_strings (help_lreg);
  140.     }
  141. }
  142.  
  143.  
  144. #if    defined(__STDC__) || defined(__GNUC__)
  145. void
  146. do_cor (FILE * fp, double *base_vector, double base_mean, double base_var, int base_rank)
  147. #else
  148. void
  149. do_cor (fp, base_vector, base_mean, base_var, base_rank)
  150.      FILE *fp;
  151.      double *base_vector;
  152.      double base_mean;
  153.      double base_var;
  154.      int base_rank;
  155. #endif
  156. {
  157.   int n;
  158.   int k;
  159.   double cross;
  160.   double var;
  161.   double mean;
  162.   double slope;
  163.   double intercept;
  164.   double diff;
  165.   double *vector;
  166.   double *ptr;
  167.  
  168.   /*
  169.    * Read in the vector
  170.    */
  171.  
  172.   vector = read_vector (fp, &n);
  173.  
  174.   if (n < base_rank)
  175.     {
  176.       message ("insufficient number of elements (%d of %d) in vector `%s'\n",
  177.            n, base_rank, vector_name);
  178.     }
  179.   else
  180.     {
  181.  
  182.       /*
  183.        * compute the cross-correlation and variance of
  184.        * both the data vector
  185.        */
  186.       mean = cross = var = 0;
  187.       ptr = vector;
  188.       for (k = 0; k < base_rank; k++)
  189.     {
  190.       mean += *ptr++;
  191.     }
  192.       mean /= base_rank;
  193.  
  194.       ptr = vector;
  195.       for (k = 0; k < base_rank; k++)
  196.     {
  197.       diff = *ptr++ - mean;
  198.       cross += (*base_vector++ - base_mean) * diff;
  199.       var += diff * diff;
  200.     }
  201.  
  202.       if (strcmp (progname, COR_NAME) == 0)
  203.     {
  204.  
  205.       /*
  206.            * Compute correlation and print it out
  207.            */
  208.       slope = cross / sqrt (base_rank * base_var * var);
  209.       print_number (stdout, slope);
  210.     }
  211.       else
  212.     {
  213.  
  214.       /*
  215.            * Compute the slope and intercept
  216.            */
  217.       slope = cross / base_var / base_rank;
  218.       intercept = mean - slope * base_mean;
  219.       if (opt_form)
  220.         {
  221.           int prec;
  222.           prec = set_precision (-1);    /* Get current precision     */
  223.           if (intercept_output)
  224.         {
  225.           printf ("-i%.*g ", prec, intercept);
  226.         }
  227.           if (slope_output)
  228.         {
  229.           printf ("-s%.*g ", prec, slope);
  230.         }
  231.         }
  232.       else
  233.         {
  234.           if (intercept_output)
  235.         {
  236.           print_number (stdout, intercept);
  237.         }
  238.           if (slope_output)
  239.         {
  240.           print_number (stdout, slope);
  241.         }
  242.         }
  243.     }
  244.     }
  245.  
  246.   free (vector);
  247. }
  248.  
  249. /*
  250.  * Compute the mean and variance of the given vector
  251.  */
  252.  
  253. #if    defined(__STDC__) || defined(__GNUC__)
  254. void
  255. calc_base_var (double *mean, double *var, double *base, int cnt)
  256. #else
  257. void
  258. calc_base_var (mean, var, base, cnt)
  259.      double *mean;
  260.      double *var;
  261.      double *base;
  262.      int cnt;
  263. #endif
  264. {
  265.   int k;
  266.   double diff;
  267.   double *base_ptr;
  268.  
  269.   *var = 0;
  270.   *mean = 0;
  271.   base_ptr = base;
  272.  
  273.   for (k = 0; k < cnt; k++)
  274.     {
  275.       *mean += *base_ptr++;
  276.     }
  277.   *mean /= cnt;
  278.  
  279.   for (k = 0; k < cnt; k++)
  280.     {
  281.       diff = *base++ - *mean;
  282.       *var += diff * diff;
  283.     }
  284.   *var /= cnt;
  285. }
  286.  
  287. int
  288. main (argc, argv)
  289.      int argc;
  290.      char *argv[];
  291. {
  292.   char *base;
  293.   const char *opt_string;
  294.   int option;
  295.   int errcnt;
  296.   int base_rank;
  297.   double base_mean;
  298.   double base_var;
  299.   double *base_vector;
  300.  
  301.   set_def_prec ();
  302.   progname = get_my_name (argv[0], DEF_PROGNAME);
  303.  
  304.   /*
  305.    * Do I know who I am?
  306.    */
  307.  
  308.   if (strcmp (progname, COR_NAME) == 0)
  309.     {
  310.       opt_string = COR_OPT_STRING;
  311.       (void) set_max_columns (1);
  312.     }
  313.   else if (strcmp (progname, LREG_NAME) == 0)
  314.     {
  315.       opt_string = LREG_OPT_STRING;
  316.       (void) set_max_columns (2);
  317.     }
  318.   else
  319.     {
  320.       message ("unknown program name.  Should be cor or lreg.\n");
  321.       exit (1);
  322.     }
  323.  
  324.   errcnt = 0;
  325.   base = NULL;
  326.  
  327.   /*
  328.    * Read options
  329.    */
  330.  
  331.   while ((option = getopt (argc, argv, opt_string)) != -1)
  332.     {
  333.       switch (option)
  334.     {
  335.     case 'H':        /* Help                 */
  336.       help ();
  337.       exit (0);
  338.       break;
  339.     case 'P':        /* Precision             */
  340.       (void) set_precision (atoi (optarg));
  341.       break;
  342.     case 'c':        /* Number of output columns     */
  343.       (void) set_max_columns (atoi (optarg));
  344.       break;
  345.     case 'F':        /* Base vector file         */
  346.       base = optarg;
  347.       break;
  348.     case 'i':        /* lreg:  intercept only     */
  349.       slope_output = FALSE;
  350.       break;
  351.     case 'o':        /* lreg:  output in option form     */
  352.       opt_form = TRUE;
  353.       break;
  354.     case 's':        /* lreg:  slope only         */
  355.       intercept_output = FALSE;
  356.       break;
  357.     case ',':
  358.     default:
  359.       break;
  360.     }            /* endswitch     */
  361.     }                /* endwhile     */
  362.  
  363.   /*
  364.    * Check options
  365.    */
  366.  
  367.   if ((base == NULL) && (optind >= argc))
  368.     {
  369.       message ("the base vector and argument vectors cannot both be from stdin\n");
  370.       errcnt++;
  371.     }
  372.  
  373.   /*
  374.    * If a base vector is given, read it from there.
  375.    * Otherwise read it from stdin.
  376.    */
  377.  
  378.   if (base == NULL)
  379.     {
  380.       base_vector = read_vector (stdin, &base_rank);
  381.     }
  382.   else
  383.     {
  384.       FILE *fp;
  385.       fp = fopen (base, "r");
  386.       if (fp == NULL)
  387.     {
  388.       message ("cannot open base file `%s'\n", base);
  389.       errcnt++;
  390.     }
  391.       else
  392.     {
  393.       base_vector = read_vector (fp, &base_rank);
  394.       (void) fclose (fp);
  395.     }
  396.     }
  397.  
  398.   /*
  399.    * Exit on errors
  400.    */
  401.  
  402.   if (errcnt > 0)
  403.     {
  404.       exit (1);
  405.     }
  406.   calc_base_var (&base_mean, &base_var, base_vector, base_rank);
  407.  
  408.   /*
  409.    * Now correlate the base vector with the given vectors
  410.    * or stdin, as appropriate.
  411.    */
  412.  
  413.   if (optind >= argc)
  414.     {
  415.       vector_name = "standard input";
  416.       do_cor (stdin, base_vector, base_mean, base_var, base_rank);
  417.     }
  418.   else
  419.     {
  420.       FILE *fp;
  421.       for (; optind < argc; optind++)
  422.     {
  423.       fp = fopen (argv[optind], "r");
  424.       if (fp == NULL)
  425.         {
  426.           message ("cannot open vector `%s'\n", argv[optind]);
  427.         }
  428.       else
  429.         {
  430.           vector_name = argv[optind];
  431.           do_cor (fp, base_vector, base_mean, base_var, base_rank);
  432.         }
  433.     }            /* endfor     */
  434.     }
  435.  
  436.   end_column (stdout);
  437.   return 0;
  438. }
  439.