2 Translated from RATFOR lowess code of W. S. Cleveland as obtained from NETLIB
10 static double pow2(x
) double x
; { return(x
* x
); }
11 static double pow3(x
) double x
; { return(x
* x
* x
); }
12 static double fmax(x
,y
) double x
, y
; { return (x
> y
? x
: y
); }
14 lowess(x
, y
, n
, f
, nsteps
, delta
, ys
, rw
, res
)
15 double *x
, *y
, f
, delta
, *ys
, *rw
, *res
;
18 int iter
, ns
, ok
, nleft
, nright
, i
, j
, last
, m1
, m2
;
19 double d1
, d2
, denom
, alpha
, cut
, cmad
, c9
, c1
, r
;
21 if (n
< 2) { ys
[0] = y
[0]; return(1); }
22 ns
= max(min((int) (f
* n
), n
), 2); /* at least two, at most n points */
23 for(iter
= 1; iter
<= nsteps
+ 1; iter
++){ /* robustness iterations */
24 nleft
= 0; nright
= ns
- 1;
25 last
= -1; /* index of prev estimated point */
26 i
= 0; /* index of current point */
28 while(nright
< n
- 1){
29 /* move nleft, nright to right if radius decreases */
31 d2
= x
[nright
+ 1] - x
[i
];
32 /* if d1 <= d2 with x[nright+1] == x[nright], lowest fixes */
34 /* radius will not decrease by move right */
38 lowest(x
, y
, n
, x
[i
], &ys
[i
], nleft
, nright
, res
, (iter
> 1), rw
, &ok
);
39 /* fitted value at x[i] */
40 if (! ok
) ys
[i
] = y
[i
];
41 /* all weights zero - copy over value (all rw==0) */
42 if (last
< i
- 1) { /* skipped points -- interpolate */
43 denom
= x
[i
] - x
[last
]; /* non-zero - proof? */
44 for(j
= last
+ 1; j
< i
; j
= j
+ 1){
45 alpha
= (x
[j
] - x
[last
]) / denom
;
46 ys
[j
] = alpha
* ys
[i
] + (1.0 - alpha
) * ys
[last
];
49 last
= i
; /* last point actually estimated */
50 cut
= x
[last
] + delta
; /* x coord of close points */
51 for(i
=last
+ 1; i
< n
; i
++) { /* find close points */
52 if (x
[i
] > cut
) break; /* i one beyond last pt within cut */
53 if(x
[i
] == x
[last
]) { /* exact match in x */
58 i
= max(last
+ 1,i
- 1);
59 /* back 1 point so interpolation within delta, but always go forward */
60 } while(last
< n
- 1);
61 for (i
= 0; i
< n
; i
++) /* residuals */
62 res
[i
] = y
[i
] - ys
[i
];
63 if (iter
> nsteps
) break; /* compute robustness weights except last time */
64 for (i
= 0; i
< n
; i
++)
67 m1
= 1 + n
/ 2; m2
= n
- m1
+ 1;
68 cmad
= 3.0 * (rw
[m1
] + rw
[m2
]); /* 6 median abs resid */
69 c9
= .999 * cmad
; c1
= .001 * cmad
;
70 for (i
= 0; i
< n
; i
++) {
72 if(r
<= c1
) rw
[i
] = 1.0; /* near 0, avoid underflow */
73 else if(r
> c9
) rw
[i
] = 0.0; /* near 1, avoid underflow */
74 else rw
[i
] = pow2(1.0 - pow2(r
/ cmad
));
81 static lowest(x
, y
, n
, xs
, ys
, nleft
, nright
, w
, userw
, rw
, ok
)
82 double *x
, *y
, *w
, *rw
, xs
, *ys
;
83 int n
, nleft
, nright
, userw
, *ok
;
85 double range
, h
, h1
, h9
, a
, b
, c
, r
;
88 range
= x
[n
- 1] - x
[0];
89 h
= fmax(xs
- x
[nleft
], x
[nright
] - xs
);
93 /* compute weights (pick up all ties on right) */
94 a
= 0.0; /* sum of weights */
95 for(j
= nleft
; j
< n
; j
++) {
98 if (r
<= h9
) { /* small enough for non-zero weight */
99 if (r
> h1
) w
[j
] = pow3(1.0-pow3(r
/h
));
101 if (userw
) w
[j
] = rw
[j
] * w
[j
];
104 else if (x
[j
] > xs
) break; /* get out at first zero wt on right */
106 nrt
= j
- 1; /* rightmost pt (may be greater than nright because of ties) */
107 if (a
<= 0.0) *ok
= FALSE
;
108 else { /* weighted least squares */
111 /* make sum of w[j] == 1 */
112 for (j
= nleft
; j
<= nrt
; j
++) w
[j
] = w
[j
] / a
;
114 if (h
> 0.0) { /* use linear fit */
116 /* find weighted center of x values */
117 for (j
= nleft
, a
= 0.0; j
<= nrt
; j
++) a
+= w
[j
] * x
[j
];
120 for (j
= nleft
, c
= 0.0; j
<= nrt
; j
++)
121 c
+= w
[j
] * (x
[j
] - a
) * (x
[j
] - a
);
123 if(sqrt(c
) > .001 * range
) {
124 /* points are spread out enough to compute slope */
126 for (j
= nleft
; j
<= nrt
; j
++)
127 w
[j
] = w
[j
] * (1.0 + b
*(x
[j
] - a
));
130 for (j
= nleft
, *ys
= 0.0; j
<= nrt
; j
++) *ys
+= w
[j
] * y
[j
];
137 if (*a
< *b
) return(-1);
138 else if (*a
> *b
) return(1);
146 qsort(x
, n
, sizeof(double), compar
);