lib/lowess.c

   1 /*
   2   Translated from RATFOR lowess code of W. S. Cleveland as obtained from NETLIB
   3 */
   4
   5 #include "xmath.h"
   6
   7 #define FALSE 0
   8 #define TRUE 1
   9
  10 extern int max(int, int);
  11 extern int min(int, int);
  12
  13 static double pow2(x) double x; { return(x * x); }
  14 static double pow3(x) double x; { return(x * x * x); }
  15 /* static */ double fmax(x,y) double x, y; { return (x > y ? x : y); }
  16
  17 int
  18 static compar(double *a, double *b)
  19 {
  20   if (*a < *b) return(-1);
  21   else if (*a > *b) return(1);
  22   else return(0);
  23 }
  24
  25 static void
  26 lowest(double *x, double *y, int n, double xs, double *ys, int nleft, int nright,
  27        double *w, int userw, double *rw, int *ok)
  28 {
  29   double range, h, h1, h9, a, b, c, r;
  30   int j, nrt;
  31
  32   range = x[n - 1] - x[0];
  33   h = fmax(xs - x[nleft], x[nright] - xs);
  34   h9 = .999 * h;
  35   h1 = .001 * h;
  36
  37   /* compute weights (pick up all ties on right) */
  38   a = 0.0;        /* sum of weights */
  39   for(j = nleft; j < n; j++) {
  40     w[j]=0.0;
  41     r = fabs(x[j] - xs);
  42     if (r <= h9) {    /* small enough for non-zero weight */
  43       if (r > h1) w[j] = pow3(1.0-pow3(r/h));
  44       else w[j] = 1.0;
  45       if (userw) w[j] = rw[j] * w[j];
  46       a += w[j];
  47     }
  48     else if (x[j] > xs) break;  /* get out at first zero wt on right */
  49   }
  50   nrt = j - 1;  /* rightmost pt (may be greater than nright because of ties) */
  51   if (a <= 0.0) *ok = FALSE;
  52   else { /* weighted least squares */
  53     *ok = TRUE;
  54
  55     /* make sum of w[j] == 1 */
  56     for (j = nleft; j <= nrt; j++) w[j] = w[j] / a;
  57
  58     if (h > 0.0) {     /* use linear fit */
  59
  60       /* find weighted center of x values */
  61       for (j = nleft, a = 0.0; j <= nrt; j++) a += w[j] * x[j];
  62
  63       b = xs - a;
  64       for (j = nleft, c = 0.0; j <= nrt; j++)
  65         c += w[j] * (x[j] - a) * (x[j] - a);
  66
  67       if(sqrt(c) > .001 * range) {
  68         /* points are spread out enough to compute slope */
  69         b = b/c;
  70         for (j = nleft; j <= nrt; j++)
  71           w[j] = w[j] * (1.0 + b*(x[j] - a));
  72       }
  73     }
  74     for (j = nleft, *ys = 0.0; j <= nrt; j++) *ys += w[j] * y[j];
  75   }
  76 }
  77
  78 static void
  79 sort(double *x, int n)
  80 {
  81   extern void qsort();
  82
  83   qsort(x, n, sizeof(double), compar);
  84 }
  85
  86
  87
  88 int
  89 lowess(double *x, double *y, int n,
  90        double f, int nsteps,
  91        double delta, double *ys, double *rw, double *res)
  92 {
  93   int iter, ns, ok, nleft, nright, i, j, last, m1, m2;
  94   double d1, d2, denom, alpha, cut, cmad, c9, c1, r;
  95
  96   if (n < 2) { ys[0] = y[0]; return(1); }
  97   ns = max(min((int) (f * n), n), 2);  /* at least two, at most n points */
  98   for(iter = 1; iter <= nsteps + 1; iter++){      /* robustness iterations */
  99     nleft = 0; nright = ns - 1;
 100     last = -1;        /* index of prev estimated point */
 101     i = 0;   /* index of current point */
 102     do {
 103       while(nright < n - 1){
 104         /* move nleft, nright to right if radius decreases */
 105         d1 = x[i] - x[nleft];
 106         d2 = x[nright + 1] - x[i];
 107         /* if d1 <= d2 with x[nright+1] == x[nright], lowest fixes */
 108         if (d1 <= d2) break;
 109         /* radius will not decrease by move right */
 110         nleft++;
 111         nright++;
 112       }
 113       lowest(x, y,
 114              n, x[i],
 115              &ys[i],
 116              nleft, nright,
 117              res, (iter > 1), rw, &ok);
 118       /* fitted value at x[i] */
 119       if (! ok) ys[i] = y[i];
 120       /* all weights zero - copy over value (all rw==0) */
 121       if (last < i - 1) { /* skipped points -- interpolate */
 122         denom = x[i] - x[last];    /* non-zero - proof? */
 123         for(j = last + 1; j < i; j = j + 1){
 124           alpha = (x[j] - x[last]) / denom;
 125           ys[j] = alpha * ys[i] + (1.0 - alpha) * ys[last];
 126         }
 127       }
 128       last = i;        /* last point actually estimated */
 129       cut = x[last] + delta;     /* x coord of close points */
 130       for(i=last + 1; i < n; i++) {     /* find close points */
 131         if (x[i] > cut) break;     /* i one beyond last pt within cut */
 132         if(x[i] == x[last]) {      /* exact match in x */
 133           ys[i] = ys[last];
 134           last = i;
 135         }
 136       }
 137       i = max(last + 1,i - 1);
 138       /* back 1 point so interpolation within delta, but always go forward */
 139     } while(last < n - 1);
 140     for (i = 0; i < n; i++)      /* residuals */
 141       res[i] = y[i] - ys[i];
 142     if (iter > nsteps) break; /* compute robustness weights except last time */
 143     for (i = 0; i < n; i++)
 144       rw[i] = fabs(res[i]);
 145     sort(rw,n);
 146     m1 = 1 + n / 2; m2 = n - m1 + 1;
 147     cmad = 3.0 * (rw[m1] + rw[m2]);      /* 6 median abs resid */
 148     c9 = .999 * cmad; c1 = .001 * cmad;
 149     for (i = 0; i < n; i++) {
 150       r = fabs(res[i]);
 151       if(r <= c1) rw[i] = 1.0;      /* near 0, avoid underflow */
 152       else if(r > c9) rw[i] = 0.0;  /* near 1, avoid underflow */
 153       else rw[i] = pow2(1.0 - pow2(r / cmad));
 154     }
 155   }
 156   return(0);
 157 }
 158
 159
 160