if((ptr)==NULL) \
( RA_error ("Cannot allocate memory") )

00421                                        :  (1/(X'X))X'   */
00422   vector y,                   /* vector of measured data   */
00423   vector * coef               /* vector of regression parameters */
00424 )
00425 
00426 {
00427 
00428   /*----- calculate regression coefficients -----*/
00429   vector_multiply (xtxinvxt, y, coef);
00430 
00431 }

00384 {
00385   const float EPSILON = 1.0e-5;      /* protection against divide by zero */
00386   float mspe;                 /* mean square error due to pure error */
00387   float sslf;                 /* error sum of squares due to lack of fit */
00388   float mslf;                 /* mean square error due to lack of fit */
00389   float flof;                 /* F-statistic for lack of fit */
00390 
00391 
00392   /*----- calculate mean sum of squares due to pure error -----*/
00393   mspe = sspe / (n - c);
00394 
00395 
00396   /*----- calculate mean sum of squares due to lack of fit -----*/
00397   sslf = sse - sspe;
00398   mslf = sslf / (c - p);
00399 
00400 
00401   /*----- calculate F-statistic for lack of fit -----*/
00402   if (mspe < EPSILON)
00403     flof = 0.0;
00404   else
00405     flof = mslf / mspe;
00406 
00407 
00408   /*----- return F-statistic for lack of fit -----*/
00409   return (flof);
00410 
00411 }

00553 {
00554   const float MAXF = 1000.0;         /* maximum value for F-statistic */
00555   const float EPSILON = 1.0e-5;      /* protection against divide by zero */
00556   float msef;                 /* mean square error for the full model */
00557   float msreg;                /* mean square due to the regression */
00558   float freg;                 /* F-statistic for the full regression model */
00559 
00560 
00561   /*----- Order of reduced model = order of full model ??? -----*/
00562   if (p <= q)   return (0.0);
00563 
00564 
00565   /*----- Calculate numerator and denominator of F-statistic -----*/
00566   msreg = (sser - ssef) / (p - q);    if (msreg < 0.0)  msreg = 0.0;
00567   msef   = ssef / (n - p);            if (msef  < 0.0)  msef  = 0.0;
00568 
00569   if (msreg > MAXF*msef)  freg = MAXF;
00570   else 
00571     if (msef < EPSILON)
00572       freg = 0.0;
00573     else
00574       freg = msreg / msef;
00575 
00576 
00577   /*----- Limit range of values for F-statistic -----*/
00578        if (freg < 0.0)   freg = 0.0;
00579   else if (freg > MAXF)  freg = MAXF;
00580 
00581 
00582   /*----- Return F-statistic for significance of the regression -----*/
00583   return (freg);
00584 
00585 }

00120                                        :  1/(X'X)  */
00121   matrix c,                   /* matrix representing GLT linear constraints */
00122   matrix * a,                 /* constant matrix for use later */
00123   matrix * cxtxinvct          /* matrix: C(1/(X'X))C' for GLT */
00124 )
00125 
00126 {
00127   matrix ct, xtxinvct, t1, t2;  /* temporary matrix calculation results */
00128   int ok;                       /* flag for successful matrix inversion */
00129 
00130 
00131 ENTRY("calc_glt_matrix") ;
00132 
00133   /*----- initialize matrices -----*/
00134   matrix_initialize (&ct);
00135   matrix_initialize (&xtxinvct);
00136   matrix_initialize (&t1);
00137   matrix_initialize (&t2);
00138 
00139 
00140   /*----- calculate the constant matrix which will be needed later -----*/
00141   matrix_transpose (c, &ct);                 /* [c]' */
00142   matrix_multiply (xtxinv, ct, &xtxinvct);   /* inv{[x]'[x]} [c]' */
00143   matrix_multiply (c, xtxinvct, cxtxinvct);  /* [c] inv{[x]'[x]} [c]' */
00144   ok = matrix_inverse_dsc (*cxtxinvct, &t2); /* inv{ [c] inv{[x]'[x]} [c]' } */
00145   if (ok)
00146     {
00147       matrix_multiply (xtxinvct, t2, &t1);
00148       matrix_multiply (t1, c, &t2);
00149       matrix_identity (xtxinv.rows, &t1);
00150       matrix_subtract (t1, t2, a);
00151     }
00152   else
00153     RA_error ("Improper C matrix  ( cannot invert C(1/(X'X))C' ) ");
00154 
00155 
00156   /*----- dispose of matrices -----*/
00157   matrix_destroy (&ct);
00158   matrix_destroy (&xtxinvct);
00159   matrix_destroy (&t1);
00160   matrix_destroy (&t2);
00161 
00162   RETURN (ok);
00163 }

00466 {
00467 
00468   /*----- calculate linear combinations -----*/
00469   vector_multiply (c, coef, lcoef);
00470 
00471 }

00078                                          :  1/(X'X)                */
00079   matrix * xtxinvxt             /* matrix:  (1/(X'X))X'            */
00080 )
00081 
00082 {
00083   matrix xt, xtx;               /* temporary matrix calculation results */
00084   int ok;                       /* flag for successful matrix inversion */
00085 
00086 ENTRY("calc_matrices") ;
00087 
00088   /*----- extract the independent variable matrix X -----*/
00089   matrix_extract (xdata, p, plist, x);
00090 
00091 
00092   /*----- calculate various matrices which will be needed later -----*/
00093   if( p <= 1 || !use_psinv ){
00094     matrix_initialize (&xt);
00095     matrix_initialize (&xtx);
00096     matrix_transpose (*x, &xt);
00097     matrix_multiply (xt, *x, &xtx);
00098     ok = matrix_inverse_dsc (xtx, xtxinv);
00099     if (ok)
00100       matrix_multiply (*xtxinv, xt, xtxinvxt);
00101     else
00102       RA_error ("Improper X matrix  (cannot invert X'X) ");
00103     matrix_destroy (&xtx);
00104     matrix_destroy (&xt);
00105   } else {
00106     matrix_psinv  (*x, xtxinv , xtxinvxt ); ok = 1 ;  /* 19 Jul 2004 */
00107   }
00108 
00109   RETURN (ok);
00110 }

00446 {
00447 
00448   /*----- calculate reduced model regression coefficients -----*/
00449   vector_multiply (a, coef, rcoef);
00450 
00451 }

00229 {
00230 #if 0
00231   vector yhat;               /* product Xb */
00232 #endif
00233   float sse;                 /* error sum of squares */
00234 
00235 
00236   /*----- initialize vectors -----*/
00237 #if 0
00238   vector_initialize (&yhat);
00239 #endif
00240 
00241 
00242   /*----- calculate the error sum of squares -----*/
00243 #if 0
00244   vector_multiply (x, b, &yhat);
00245   vector_subtract (y, yhat, e);
00246   sse = vector_dot (*e, *e);
00247 #else
00248   sse = vector_multiply_subtract( x , b , y , e ) ;
00249 #endif
00250 
00251 
00252   /*----- dispose of vectors -----*/
00253 #if 0
00254   vector_destroy (&yhat);
00255 #endif
00256 
00257 
00258   /*----- return SSE -----*/
00259   return (sse);
00260  
00261 }

00599 {
00600   const float EPSILON = 1.0e-5;     /* protection against divide by zero */
00601   float rsqr;                       /* coeff. of multiple determination R^2  */
00602 
00603 
00604   /*----- coefficient of multiple determination R^2 -----*/
00605   if (sser < EPSILON)
00606     rsqr = 0.0;
00607   else
00608     rsqr = (sser - ssef) / sser;
00609 
00610 
00611   /*----- Limit range of values for R^2 -----*/
00612        if (rsqr < 0.0)   rsqr = 0.0;
00613   else if (rsqr > 1.0)   rsqr = 1.0;
00614 
00615 
00616   /*----- Return coefficient of multiple determination R^2 -----*/
00617   return (rsqr);
00618 }

00178 {
00179 #if 0
00180   vector yhat;               /* product Xb */
00181 #endif
00182   vector e;                  /* vector of residuals */
00183   float sse;                 /* error sum of squares */
00184 
00185 
00186   /*----- initialize vectors -----*/
00187 #if 0
00188   vector_initialize (&yhat);
00189 #endif
00190   vector_initialize (&e);
00191 
00192 
00193   /*----- calculate the error sum of squares -----*/
00194 #if 0
00195   vector_multiply (x, b, &yhat);
00196   vector_subtract (y, yhat, &e);
00197   sse = vector_dot (e, e);
00198 #else
00199   sse = vector_multiply_subtract( x , b , y , &e ) ;
00200 #endif
00201 
00202 
00203   /*----- dispose of vectors -----*/
00204   vector_destroy (&e);
00205 #if 0
00206   vector_destroy (&yhat);
00207 #endif
00208 
00209 
00210   /*----- return SSE -----*/
00211   return (sse);
00212  
00213 }

00279 {
00280   vector yhat;               /* product Xb */
00281   vector e;                  /* vector of residuals */
00282   float sse;                 /* error sum of squares */
00283   int it;                    /* time point index */
00284 
00285 
00286   /*----- initialize vectors -----*/
00287   vector_initialize (&yhat);
00288   vector_initialize (&e);
00289 
00290 
00291   /*----- calculate the error sum of squares -----*/
00292   vector_multiply (x, b, &yhat);
00293   vector_subtract (y, yhat, &e);
00294   sse = vector_dotself( e );
00295 
00296   /*----- save the fitted time series and residual errors -----*/
00297   for (it = 0;  it < x.rows;  it++)
00298     {
00299       fitts[it] = yhat.elts[it];
00300       errts[it] = e.elts[it];
00301     }
00302 
00303 
00304   /*----- dispose of vectors -----*/
00305   vector_destroy (&e);
00306   vector_destroy (&yhat);
00307 
00308 
00309   /*----- return SSE -----*/
00310   return (sse);
00311  
00312 }

00328 {
00329   register int i, j;                  /* indices */
00330   register float * sum = NULL;        /* sum of observations at each level */
00331   register float diff;                /* difference between observation and average */
00332   register float sspe;                /* pure error sum of squares */
00333 
00334 
00335   /*----- initialize sum -----*/
00336   sum = (float *) malloc (sizeof(float) * c);
00337   MTEST (sum);
00338   
00339   for (j = 0;  j < c;  j++)
00340     sum[j] = 0.0;
00341 
00342 
00343   /*----- accumulate sum for each level -----*/
00344   for (i = 0;  i < y.dim;  i++)
00345     {
00346       j = levels[i];
00347       sum[j] += y.elts[i];
00348     }
00349 
00350 
00351   /*----- calculate SSPE -----*/
00352   sspe = 0.0;
00353   for (i = 0;  i < y.dim;  i++)
00354     {
00355       j = levels[i];
00356       diff = y.elts[i] - (sum[j]/counts[j]);
00357       sspe += diff * diff;
00358     }
00359 
00360   
00361   free (sum);   sum = NULL;
00362 
00363 
00364   /*----- return SSPE -----*/
00365   return (sspe);
00366  
00367 }

00485                                        :  1/(Xf'Xf)        */
00486   vector coef,                /* vector of regression parameters */
00487   vector * scoef,             /* std. devs. for regression parameters */
00488   vector * tcoef              /* t-statistics for regression parameters */
00489 )
00490 
00491 {
00492   const float MAXT = 1000.0;         /* maximum value for t-statistic */
00493   const float EPSILON = 1.0e-5;      /* protection against divide by zero */
00494   int df;                     /* error degrees of freedom */
00495   float mse;                  /* mean square error */
00496   register int i;             /* parameter index */
00497   float stddev;               /* standard deviation for parameter estimate */
00498   float tstat;                /* t-statistic for parameter estimate */
00499   float num;                  /* numerator of t-statistic */
00500   float var;                  /* variance for parameter estimate */
00501 
00502 
00503   /*----- Create vectors -----*/
00504   vector_create (p, scoef);
00505   vector_create (p, tcoef);
00506 
00507 
00508   /*----- Calculate mean square error -----*/
00509   df = n - p;
00510   mse = sse / df;
00511 
00512 
00513   for (i = 0;  i < xtxinv.rows;  i++)
00514     {
00515       /*----- Calculate standard deviation for regression parameters -----*/
00516       var = mse * xtxinv.elts[i][i];
00517       if (var <= 0.0) stddev = 0.0;
00518       else            stddev = sqrt (var);
00519       scoef->elts[i] = stddev;
00520 
00521 
00522       /*----- Calculate t-statistic for regression parameters -----*/
00523       num = coef.elts[i];
00524            if (num >  MAXT*stddev) tstat = MAXT;
00525       else if (num < -MAXT*stddev) tstat = -MAXT;
00526            else if (stddev < EPSILON)   tstat = 0.0;
00527            else                         tstat = num / stddev;
00528 
00529 
00530       /*----- Limit range of values for t-statistic -----*/
00531       if (tstat < -MAXT)  tstat = -MAXT;
00532       if (tstat >  MAXT)  tstat = MAXT;
00533       
00534       tcoef->elts[i] = tstat;
00535     }
00536 }

00183 {
00184   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00185   exit(1);
00186 }