2 Translated from RATFOR lowess code of W. S. Cleveland as obtained from NETLIB
10 extern int max(int, int);
11 extern int min(int, int);
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
); }
18 static compar(double *a
, double *b
)
20 if (*a
< *b
) return(-1);
21 else if (*a
> *b
) return(1);
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
)
29 double range
, h
, h1
, h9
, a
, b
, c
, r
;
32 range
= x
[n
- 1] - x
[0];
33 h
= fmax(xs
- x
[nleft
], x
[nright
] - xs
);
37 /* compute weights (pick up all ties on right) */
38 a
= 0.0; /* sum of weights */
39 for(j
= nleft
; j
< n
; j
++) {
42 if (r
<= h9
) { /* small enough for non-zero weight */
43 if (r
> h1
) w
[j
] = pow3(1.0-pow3(r
/h
));
45 if (userw
) w
[j
] = rw
[j
] * w
[j
];
48 else if (x
[j
] > xs
) break; /* get out at first zero wt on right */
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 */
55 /* make sum of w[j] == 1 */
56 for (j
= nleft
; j
<= nrt
; j
++) w
[j
] = w
[j
] / a
;
58 if (h
> 0.0) { /* use linear fit */
60 /* find weighted center of x values */
61 for (j
= nleft
, a
= 0.0; j
<= nrt
; j
++) a
+= w
[j
] * x
[j
];
64 for (j
= nleft
, c
= 0.0; j
<= nrt
; j
++)
65 c
+= w
[j
] * (x
[j
] - a
) * (x
[j
] - a
);
67 if(sqrt(c
) > .001 * range
) {
68 /* points are spread out enough to compute slope */
70 for (j
= nleft
; j
<= nrt
; j
++)
71 w
[j
] = w
[j
] * (1.0 + b
*(x
[j
] - a
));
74 for (j
= nleft
, *ys
= 0.0; j
<= nrt
; j
++) *ys
+= w
[j
] * y
[j
];
79 sort(double *x
, int n
)
83 qsort(x
, n
, sizeof(double), compar
);
89 lowess(double *x
, double *y
, int n
,
91 double delta
, double *ys
, double *rw
, double *res
)
93 int iter
, ns
, ok
, nleft
, nright
, i
, j
, last
, m1
, m2
;
94 double d1
, d2
, denom
, alpha
, cut
, cmad
, c9
, c1
, r
;
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 */
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 */
109 /* radius will not decrease by move right */
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
];
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 */
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
]);
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
++) {
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
));