src/tools/perftool/linreg.c

   1 /*
   2  *------------------------------------------------------------------
   3  * Copyright (c) 2006-2016 Cisco and/or its affiliates.
   4  * Licensed under the Apache License, Version 2.0 (the "License");
   5  * you may not use this file except in compliance with the License.
   6  * You may obtain a copy of the License at:
   7  *
   8  *     http://www.apache.org/licenses/LICENSE-2.0
   9  *
  10  * Unless required by applicable law or agreed to in writing, software
  11  * distributed under the License is distributed on an "AS IS" BASIS,
  12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  13  * See the License for the specific language governing permissions and
  14  * limitations under the License.
  15  */
  16
  17 /* see "Numerical Recipies in C, 2nd ed." p 665 */
  18
  19 #include <stdio.h>
  20 #include <math.h>
  21
  22 /*
  23  * linreg
  24  * Linear regression of (xi, yi), returns parameters for least-squares
  25  * fit y = a + bx.  Also, compute Pearson's R.
  26  */
  27 void linreg (double *x, double *y, int nitems, double *a, double *b,
  28              double *minp, double *maxp, double *meanp, double *r)
  29 {
  30     double sx = 0.0;
  31     double sy = 0.0;
  32     double st2 = 0.0;
  33     double min = y[0];
  34     double max = 0.0;
  35     double ss, meanx, meany, t;
  36     double errx, erry, prodsum, sqerrx, sqerry;
  37     int i;
  38
  39     *b = 0.0;
  40
  41     for (i = 0; i < nitems; i++) {
  42         sx += x[i];
  43         sy += y[i];
  44         if (y[i] < min)
  45             min = y[i];
  46         if (y[i] > max)
  47             max = y[i];
  48     }
  49     ss = nitems;
  50     meanx = sx / ss;
  51     meany = *meanp = sy / ss;
  52     *minp = min;
  53     *maxp = max;
  54
  55     for (i = 0; i < nitems; i++) {
  56         t = x[i] - meanx;
  57         st2 += t*t;
  58         *b += t*y[i];
  59     }
  60
  61     *b /= st2;
  62     *a = (sy-sx*(*b))/ss;
  63
  64     prodsum = 0.0;
  65     sqerrx = 0.0;
  66     sqerry = 0.0;
  67
  68     /* Compute numerator of Pearson's R */
  69     for (i = 0; i < nitems; i++) {
  70         errx = x[i] - meanx;
  71         erry = y[i] - meany;
  72         prodsum += errx * erry;
  73         sqerrx += errx*errx;
  74         sqerry += erry*erry;
  75     }
  76
  77     *r = prodsum / (sqrt(sqerrx)*sqrt(sqerry));
  78 }