刘思喆 @ 贝吉塔行星

function (x, y, h, n = 25, lims = c(range(x), range(y)))
{
    nx <- length(x)
    if (length(y) != nx)
        stop("data vectors must be the same length")
    gx <- seq(lims[1], lims[2], length = n)
    gy <- seq(lims[3], lims[4], length = n)
    if (missing(h))
        h <- c(bandwidth.nrd(x), bandwidth.nrd(y))
    h <- h/4
    ax <- outer(gx, x, "-")/h[1]
    ay <- outer(gy, y, "-")/h[2]
    z <- matrix(dnorm(ax), n, nx) %*% t(matrix(dnorm(ay), n,
        nx))/(nx * h[1] * h[2])
    return(list(x = gx, y = gy, z = z))
}

bandwidth.nrd(x) <-
function (x)
{
    r <- quantile(x, c(0.25, 0.75))
    h <- (r[2] – r[1])/1.34
    4 * 1.06 * min(sqrt(var(x)), h) * length(x)^(-1/5)
}