home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 7 / 07.iso / c / c129 / 1.ddi / MULREG.C < prev    next >
Encoding:
C/C++ Source or Header  |  1989-09-15  |  3.2 KB  |  136 lines
  1. #include "gj.h"
  2. #include <math.h>
  3. #include <stdlib.h>
  4.  
  5. void  MatTxTiX(float *aryin,int numrow,int numcol,float *aryout)
  6. {
  7.    int i;
  8.    int j;
  9.    int k;
  10.  
  11.    for ( i = 0; i <= numcol-1; ++i ) {
  12.       for ( j = 0; j <= numcol-1; ++j ) {
  13.          aryout[i*numcol+j] = 0.0;
  14.          for ( k = 0; k <= numrow-1; ++k ) {
  15.             aryout[i*numcol+j] = aryout[i*numcol+j] +
  16.                                  aryin[k*numcol+j] * aryin[k*numcol+i];
  17.          }
  18.       }
  19.    }
  20. }
  21.  
  22. void MatYTiX(float *aryy,float *aryx,int numrow,
  23.             int numcol,float *aryout)
  24. {
  25.    int i;
  26.    int j;
  27.    int k;
  28.  
  29.    for ( i = 0; i <= numcol-1; ++i ) {
  30.       aryout[i] = 0.0;
  31.       for ( k = 0; k <= numrow-1; ++k ) {
  32.          aryout[i] = aryout[i] + aryy[k] * aryx[k*numcol+i];
  33.       }
  34.    }
  35. }
  36.  
  37. void ResAnalysis(float *regmat,
  38.                 float *ymat,float *regcoef,float *aryinv,
  39.                 int numrow,int numcol,float *yest,float *resid,
  40.                 float *see,float *coefsig,float *rsq,float *r)
  41. {
  42.    int i;
  43.    int j;
  44.    int nxx;
  45.    float sumressq;
  46.    float sumy;
  47.    float sumysq;
  48.    char c;
  49.  
  50.    sumressq = 0.0;
  51.    sumy = 0.0;
  52.    sumysq = 0.0;
  53.    for ( i = 0; i <= numrow-1; ++i ) {
  54.       yest[i] = 0.0;
  55.       for ( j = 0; j <= numcol-1; ++j ) {
  56.          yest[i] = yest[i] + regcoef[j] * regmat[i*numcol+j];
  57.  
  58.       }
  59.       resid[i] = yest[i] - ymat[i];
  60.       sumressq = sumressq + resid[i]*resid[i];
  61.       sumy = sumy + ymat[i];
  62.       sumysq = sumysq + ymat[i]*ymat[i];
  63.    }
  64.    if ( numrow == numcol ) {
  65.       nxx = 1;
  66.    }
  67.    else {
  68.       nxx = numrow - numcol;
  69.    }
  70.    (*see) = sqrt(sumressq / nxx);
  71.    for ( i = 0; i <= numcol-1; ++i ) {
  72.       coefsig[i] = (*see) * sqrt(fabs(aryinv[i*numcol+i]));
  73.    }
  74.    (*rsq) = (1.0 - sumressq / (sumysq - (sumy*sumy) / numrow));
  75.    (*r) = sqrt((*rsq));
  76. }
  77.  
  78. void MultipleReg(float *indvardat,float *depvardat,int numiv,int numobs,
  79.             float *regcoef,float *yest,float *resid,float *see,
  80.             float *coefsig,float *rsq,float *r,char *regerror)
  81. {
  82.  
  83.    char   regerrorxx;
  84.    float  rxx;
  85.    float  rsqxx;
  86.    float  seexx;
  87.  
  88.    float  *regmat;
  89.    float  *aryinv;
  90.    float  *scmat;
  91.    float  *arya;
  92.    float  *aryg;
  93.  
  94.    int numcol;
  95.    int i;
  96.    int j;
  97.    int k;
  98.    float regdet;
  99.    char c;
  100.    
  101.    numcol = numiv + 1;
  102.  
  103.    regmat = (float *) calloc((numobs)*numcol,4);
  104.    aryinv = (float *) calloc(numcol*numcol,4);
  105.    arya = (float *) calloc(numcol*numcol,4);
  106.    aryg = (float *) calloc(numcol,4);
  107.  
  108.    for ( i = 0; i <= numobs-1; ++i ) {
  109.       for ( j = 0; j <= numiv-1; ++j ) {
  110.          regmat[i*numcol + j + 1] = indvardat[i*numiv + j];
  111.       }
  112.       regmat[i*numcol] = 1;
  113.    }
  114.    MatTxTiX(regmat,numobs,numcol,arya);
  115.    MatYTiX(depvardat,regmat,numobs,numcol,aryg);
  116.  
  117.    GaussJordan(arya,aryg,numcol,regcoef,aryinv,®det);
  118.    
  119.    ResAnalysis(regmat,depvardat,regcoef,aryinv,numobs,
  120.                         numcol,yest,resid,&seexx,coefsig,&rsqxx,&rxx);
  121.                         
  122.  
  123.  
  124.    *see = seexx;
  125.    *rsq = rsqxx;
  126.    *r   = rxx;
  127.    if ( regdet < 1.0e-7 ) {
  128.       (*regerror) = 1;
  129.    }
  130.    else {
  131.       (*regerror) = 0;
  132.    }
  133.    free(regmat); free(aryinv); free(arya); free(aryg);
  134. }
  135.  
  136.