[quote]
引用第1楼rtist于2007-01-17 19:58发表的“”:
选择窗宽?[/quote]
不像。我贴一下KernSmooth里面叫做bkde的程序
<br />
function (x, kernel = "normal", canonical = FALSE, bandwidth, <br />
gridsize = 401, range.x, truncate = TRUE) <br />
{<br />
n <- length(x)<br />
M <- gridsize<br />
if (kernel == "normal") <br />
del0 <- (1/(4 * pi))^(1/10)<br />
if (kernel == "box") <br />
del0 <- (9/2)^(1/5)<br />
if (kernel == "epanech") <br />
del0 <- 15^(1/5)<br />
if (kernel == "biweight") <br />
del0 <- 35^(1/5)<br />
if (kernel == "triweight") <br />
del0 <- (9450/143)^(1/5)<br />
if (missing(bandwidth)) {<br />
if (canonical) {<br />
bandwidth <- (243/(35 * n))^(1/5) * sqrt(var(x))<br />
}<br />
else {<br />
bandwidth <- del0 * (243/(35 * n))^(1/5) * sqrt(var(x))<br />
}<br />
}<br />
h <- bandwidth<br />
if (canonical) {<br />
if (kernel == "normal") {<br />
tau <- 4 * del0<br />
}<br />
else {<br />
tau <- del0<br />
}<br />
}<br />
else {<br />
if (kernel == "normal") {<br />
tau <- 4<br />
}<br />
else {<br />
tau <- 1<br />
}<br />
}<br />
if (missing(range.x)) {<br />
range.x <- c(min(x) - tau * h, max(x) + tau * h)<br />
}<br />
a <- range.x[1]<br />
b <- range.x[2]<br />
gpoints <- seq(a, b, length = M)<br />
gcounts <- linbin(x, gpoints, truncate)<br />
L <- min(floor(tau * h * (M - 1)/(b - a)), M)<br />
lvec <- (0:L)<br />
delta <- (b - a)/(h * (M - 1))<br />
if (canonical == FALSE) <br />
del0 <- 1<br />
if (kernel == "normal") {<br />
kappa <- dnorm(lvec * delta/del0)/(n * h * del0)<br />
}<br />
else if (kernel == "box") {<br />
kappa <- 0.5 * dbeta(0.5 * (lvec * delta/del0 + 1), 1, <br />
1)/(n * h * del0)<br />
}<br />
else if (kernel == "epanech") {<br />
kappa <- 0.5 * dbeta(0.5 * (lvec * delta/del0 + 1), 2, <br />
2)/(n * h * del0)<br />
}<br />
else if (kernel == "biweight") {<br />
kappa <- 0.5 * dbeta(0.5 * (lvec * delta/del0 + 1), 3, <br />
3)/(n * h * del0)<br />
}<br />
else if (kernel == "triweight") {<br />
kappa <- 0.5 * dbeta(0.5 * (lvec * delta/del0 + 1), 4, <br />
4)/(n * h * del0)<br />
}<br />
P <- 2^(ceiling(log(M + L)/log(2)))<br />
kappa <- c(kappa, rep(0, P - 2 * L - 1), kappa[(L + 1):2])<br />
gcounts <- c(gcounts, rep(0, P - M))<br />
kappa <- fft(kappa)<br />
gcounts <- fft(gcounts)<br />
return(list(x = gpoints, y = (Re(fft(kappa * gcounts, TRUE))/P)[1:M]))<br />
}<br />
<environment: namespace:KernSmooth><br />
<br />