2 Translated from RATFOR lowess code of W. S. Cleveland as obtained from NETLIB
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
); }
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);
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
;
41 range
= x
[n
- 1] - x
[0];
42 h
= fmax(xs
- x
[nleft
], x
[nright
] - xs
);
46 /* compute weights (pick up all ties on right) */
47 a
= 0.0; /* sum of weights */
48 for(j
= nleft
; j
< n
; j
++) {
51 if (r
<= h9
) { /* small enough for non-zero weight */
52 if (r
> h1
) w
[j
] = pow3(1.0-pow3(r
/h
));
54 if (userw
) w
[j
] = rw
[j
] * 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 */
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
];
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 */
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
];
88 sort(double *x
, size_t n
)
92 qsort(x
, n
, sizeof(double), compar
);
98 lowess(double *x
, double *y
, size_t n
,
99 double f
, size_t nsteps
,
100 double delta
, double *ys
, double *rw
, double *res
)
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 */
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 */
119 /* radius will not decrease by move right */
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 */
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
]);
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
++) {
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
));