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