src/c/lowess.c

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