4 #define ROOT2PI 2.50662827463100050241
20 kernel(double x
, double y
, double w
, int type
)
30 if (-0.5 < z
&& z
< 0.5) {
31 z
= (1.0 - 4 * z
* z
);
32 k
= 15.0 * z
* z
/ 8.0;
39 k
= exp(- 0.5 * z
* z
) / ROOT2PI
;
44 k
= (fabs(z
) < 0.5) ? 1.0 : 0.0;
47 if (-1.0 <= z
&& z
< 0.0) k
= 1.0 + z
;
48 else if (0.0 <= z
&& z
< 1.0) k
= 1.0 - z
;
51 default: k
= 0.0; break;
61 kernel_smooth(double *x
, double *y
, size_t n
,
62 double width
, double *wts
, double *wds
,
63 double *xs
, double *ys
, size_t ns
, int ktype
)
66 double wsum
, ysum
, lwidth
, lwt
, xmin
, xmax
;
71 for (xmin
= xmax
= x
[0], i
= 1; i
< n
; i
++) {
72 if (xmin
> x
[i
]) xmin
= x
[i
];
73 if (xmax
< x
[i
]) xmax
= x
[i
];
75 width
= (xmax
- xmin
) / (1 + log((double) n
));
78 for (i
= 0; i
< ns
; i
++) {
79 for (j
= 0, wsum
= 0.0, ysum
= 0.0; j
< n
; j
++) {
80 lwidth
= (wds
!= nil
) ? width
* wds
[j
] : width
;
81 lwt
= kernel(xs
[i
], x
[j
], lwidth
, ktype
);
82 if (wts
!= nil
) lwt
*= wts
[j
];
84 if (y
!= nil
) ysum
+= lwt
* y
[j
];
86 if (y
!= nil
) ys
[i
] = (wsum
> 0.0) ? ysum
/ wsum
: 0.0;
87 else ys
[i
] = wsum
/ n
;