home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / math / ols / pretty.c < prev    next >
C/C++ Source or Header  |  1993-07-28  |  6KB  |  226 lines
  1. /* (Not quite) Routines for doing pretty things ...  */
  2.  
  3. #include <stdio.h>
  4. #include <math.h>
  5. #include <stdlib.h>
  6.  
  7. #include "utils.h"
  8. #include "linstats.h"
  9. #include "distribs.h"
  10.  
  11. /* Prototypes: */
  12.  
  13. void GenAwk (double *beta, int K,
  14.          char **labels, int *maptable, int labelcount);
  15. void GenEPP (char **labels, int K, double *beta, double *S, double sigma2);
  16.  
  17. /* ols(...) is the bottom line doer of all work, called by main() */
  18.  
  19. int 
  20. ols (char **labels, float **data, int K, int ObsNum,
  21.      int purebeta, int pscript, int epp,
  22.      int *maptable, int labelcount)
  23.      /* calls olsengine for computations and does display. */
  24. {
  25.   double *beta, *S;
  26.   float *Yhat;
  27.   double sigma2, R2, probv, F;
  28.   int noinference, i, error, dof1, dof2;
  29.   char *DASHLINE =
  30.   "-------------+-----------------------------------------------";
  31.  
  32.   if ((beta = (double *) malloc (K * sizeof (double))) == NULL)
  33.     {
  34.       fprintf (stderr, "Out of memory when allocating beta.\n");
  35.       return 1;
  36.     }
  37.   if ((S = (double *) malloc (K * sizeof (double))) == NULL)
  38.     {
  39.       fprintf (stderr, "Out of memory when allocating S.\n");
  40.       return 1;
  41.     }
  42.   if ((Yhat = (float *) malloc (ObsNum * sizeof (float))) == NULL)
  43.     {
  44.       fprintf (stderr, "Out of memory when allocating Yhat.\n");
  45.       return 1;
  46.     }
  47.  
  48.   /* if pscript is on, you don't need to do inference. */
  49.   noinference = FALSE;
  50.   if (pscript)
  51.     noinference = TRUE;
  52.   if (1 == (error = olsengine (noinference, data, K, ObsNum,
  53.                    beta, S, &sigma2, &R2, &F, Yhat)))
  54.     {
  55.       fprintf (stderr, "Errors in OLS calculations.\n");
  56.       return 1;
  57.     }
  58.  
  59.   if (purebeta)
  60.     {
  61.       for (i = 0; i < K; i++)
  62.     printf (" %f ", beta[i]);
  63.       printf ("%f\n", sqrt (sigma2));
  64.       return 0;
  65.     }
  66.   if (pscript)
  67.     {
  68.       GenAwk (beta, K, labels, maptable, labelcount);
  69.       return 0;
  70.     }
  71.   if (epp)
  72.     {
  73.       GenEPP (labels, K, beta, S, sigma2);
  74.       return 0;
  75.     }
  76.  
  77.   /* Now start printing results. */
  78.   printf ("%d observations, %d rhs regressors.\n", ObsNum, K);
  79.   printf ("%s\n", DASHLINE);
  80.   printf ("%12s | %12s    %12s  %8s  %8s\n",
  81.       "Variable",
  82.       "Coefficient",
  83.       "Std.Error",
  84.       "t   ",
  85.       "Prob>|t|");
  86.   printf ("%s\n", DASHLINE);
  87.  
  88.   for (i = 0; i < K; i++)
  89.     {
  90.       if (S[i] == 0.0)
  91.     probv = 0.0;
  92.       else
  93.     {
  94.       if (1 == tprob ((double) -fabs (beta[i] / S[i]), ObsNum - K, &probv))
  95.         {
  96.           fprintf (stderr, "Fatal error in computing t probability.\n");
  97.           return 1;
  98.         }
  99.     }
  100.       printf ("%12s | %12.6f  %12.6f  %8.3f  %8.3f\n",
  101.           labels[i],
  102.           beta[i],
  103.           S[i],
  104.           (S[i] == 0.0 ? 0.0 : beta[i] / S[i]),
  105.           (probv * 2.0)
  106.     );
  107.     }
  108.   printf ("%s\n", DASHLINE);
  109.  
  110.   dof1 = K - 1;
  111.   dof2 = ObsNum - K;        /* F(dof1, dof2) */
  112.   printf ("R2 = %5.3f, F(%d, %d) = %3.1f, Prob>F = %5.3f\n",
  113.       R2, dof1, dof2, F, pF (F, dof1, dof2));
  114.   printf ("RMSE = %5.3f\n", sqrt (sigma2));
  115.   return 0;
  116. }
  117.  
  118. void 
  119. GenAwk (double *beta, int K,
  120.     char **labels, int *maptable, int labelcount)
  121.      /* Recall layout of maptable:
  122.     entries are defined from 1..K (not 1..labelcount),
  123.     (maptable[i] == j) ==>
  124.         i'th field of beta comes from j'th field of input line.
  125. */
  126. {
  127.   int i;
  128.   char *ylabel;
  129.  
  130.   printf ("\n#\n# This awk program was generated by ols.\n");
  131.   printf ("# It reads data from StdIn and prints predictions to StdOut.\n");
  132.   printf ("#\n");
  133.   printf ("# It is built for the data layout used when estimating.\n");
  134.   printf ("# However, this should not be hard to modify.\n");
  135.   printf ("#\n\n");
  136.  
  137.   printf ("{\n");
  138.   for (i = 0; i < K; i++)
  139.     printf ("\t%s = $%d;\n", labels[i], (1 + maptable[i]));
  140.   ylabel = labels[K];
  141.   printf ("\t%s = $%d;\n", ylabel, (1 + maptable[K]));
  142.  
  143.   printf ("\n\tpredicted = ");
  144.   for (i = 0; i < K; i++)
  145.     {
  146.       printf ("(%f*%s)", beta[i], labels[i]);
  147.       if (i != (K - 1))
  148.     printf (" \\\n\t\t+ ");
  149.       else
  150.     printf (";\n");
  151.     }
  152.  
  153.   /* Now to find y. */
  154.   printf ("\terror = %s - predicted;\n", ylabel);
  155.   printf ("\tprint predicted \" \" %s \" \" error;\n", ylabel);
  156.   printf ("\tSSE += (error*error);\n");
  157.   printf ("}\n");
  158.  
  159.   printf ("\nEND {\n");
  160.   printf ("\tMSE = SSE/NR;\n");
  161.   printf ("\tprint \"MSE = \" MSE \", SSE = \" SSE > \"/dev/tty\";\n");
  162.   printf ("}\n\n");
  163.   printf ("# End of generated script.\n\n");
  164. }
  165.  
  166. /* Example of generated awk, when labels are unknown:
  167.  * #
  168.  * # This awk program was generated by ols.
  169.  * # It reads data from StdIn and prints predictions to StdOut.
  170.  * #
  171.  * # It is built for the data layout used when estimating.
  172.  * # However, this should not be hard to modify.
  173.  * #
  174.  *
  175.  * {
  176.  *     $1 = $1;
  177.  *     $2 = $2;
  178.  *     $3 = $3;
  179.  *         # In out-of-sample prediction, $3 defaults to 0.0.
  180.  *
  181.  *     predicted = (1.500000*$1) \
  182.  *         + (0.700000*$2);
  183.  *     error = $3 - predicted;
  184.  *     print predicted " " $3 " " error;
  185.  *     SSE += (error*error);
  186.  * }
  187.  *
  188.  * END {
  189.  *     MSE = SSE/NR;
  190.  *     print "MSE = " MSE ", SSE = " SSE > "/dev/tty";
  191.  * }
  192.  *
  193.  * # End of script.
  194.  */
  195.  
  196. void 
  197. GenEPP (char **labels, int K, double *beta, double *S, double sigma2)
  198.      /* Given estimation results, generate EPP.  This can be
  199. post-processed into TeX for sure. */
  200. {
  201.   int i;
  202.  
  203.   printf ("caption Results of OLS Regression for %s\n", labels[K]);
  204.   printf ("comment Regression Coefficients\n");
  205.   for (i = 0; i < K; i++)
  206.     printf ("estimate %s yes %12.6f %8.3f\n",
  207.         labels[i], beta[i], S[i]);
  208.   printf ("blank\n");
  209.   printf ("comment Regression Standard Error\n");
  210.   printf ("estimate $\\sigma$ yes %5.3f not-estimated\n", sqrt (sigma2));
  211. }
  212.  
  213. /* Example:
  214.  * caption Typical display of results after OLS
  215.  * comment Regression Coefficients
  216.  * estimate $\beta$ yes -9.209 0.0286
  217.  * estimate $\tau$ yes -6.000 0.234
  218.  * estimate $\rho$ yes -1.078 0.0288
  219.  * estimate $\delta$ yes -2.303 0.21
  220.  * estimate $\theta$ yes -9.000 0.113
  221.  * estimate $\xi$ yes -9.000 1.00
  222.  * blank
  223.  * comment Regression Standard Error
  224.  * estimate $\sigma_v$ no 0.867
  225.  */
  226.