有效地找到许多任意样本值的经验密度（）（如 dnorm（），但用于经验分布）

Question

假设您已经为 x.sample 的样本定义了经验密度 (sample.density)，如下所示：

set.seed(1)
x.sample <- rnorm(100)
sample.density <- density(x.sample)

现在假设我们有一个梯度 G，我们想知道它的预期密度

G <- seq(-2,2, length.out=20)

根据经验分布sample.density，G中每个值的密度是多少？

如果我使用 for() 循环，我可以得到这样的答案：

G.dens <- c()
for(i in 1:length(G)){
    t.G <- G[i]
    G.dens[i] <- sample.density$y[which.min(abs(sample.density$x-t.G))]
}

总体想法是做类似 dnorm() 的事情，但我不想假设 x 是具有指定均值和 s.d 的正态分布，我想使用根据任意样本凭经验确定的分布（不一定是正态分布或统计包中的任何其他描述良好的分布）。

Answer 1

我认为@MrFlick 的评论为我指明了正确的方向。除了建议的 approxfun 方法和我的示例 for 循环方法之外，我还意识到我可以使用 mapply。请注意，我对 approxfun 的使用与使用 which.min 的其他 2 种方法的结果并不完全匹配，但我并不担心输出差异太大，尽管其他人可能会。

First, reproducing the sample data from the question:
set.seed(1)
x.sample <- rnorm(100)
sample.density <- density(x.sample)
G <- seq(-2,2, length.out=20)

现在，为循环版本创建一个函数

循环()

loop <- function(x, ref){
    if(class(ref)!="density"){
        ref <- density(ref)
    }

    ref.y <- ref$y
    ref.x <- ref$x

    G.dens <- c()
    for(i in 1:length(x)){
        t.G <- x[i]
        G.dens[i] <- ref.y[which.min(abs(ref.x-t.G))]
    }

    G.dens
}

接下来，我将使用我提出的使用 mapply

的方法

dsample()

dsample <- function(x, ref){

    if(class(ref)!="density"){
        ref <- density(ref)
    }

    XisY <- function(x,y){ # which of several X values most closely matches a single Y value?
        which.min(abs(y-x))
    }

    ref.y <- ref$y
    ref.x <- ref$x

    # ds <- approxfun(ref)

    # ds(x)

    ref.y[mapply(XisY, x, MoreArgs=list(y=ref.x))]
}

最后，按照@MrFlick 的建议利用 approxfun 的方法：

af()

af <- function(x, ref){
    if(class(ref)!="density"){
        ref <- density(ref)
    }

    # XisY <- function(x,y){ # which of several X values most closely matches a single Y value?
    #   which.min(abs(y-x))
    # }

    ref.y <- ref$y
    ref.x <- ref$x

    ds <- approxfun(ref)

    ds(x)

    # ref.y[mapply(XisY, x, MoreArgs=list(y=ref.x))]
}

现在比较一下他们的速度：

microbenchmark(
    loop(G, sample.density),
    dsample(G, sample.density),
    af(G, sample.density)
)
# Unit: microseconds
#                        expr     min       lq     mean  median       uq      max neval
#     loop(G, sample.density) 221.801 286.6675 360.3698 348.065 409.9785  942.071   100
#  dsample(G, sample.density) 252.641 290.7900 359.0432 368.388 417.1510  520.960   100
#       af(G, sample.density) 201.331 227.8740 276.0425 253.339 273.6225 2545.081   100

现在比较随着 G 大小的增加速度：

speed.loop <- c()
speed.dsample <- c()
speed.af <- c()
lengths <- seq(20, 5E3, by=200)
for(i in 1:length(lengths)){
    G <- seq(-2,2, length.out=lengths[i])
    bm <- microbenchmark(
        loop(G, sample.density),
        dsample(G, sample.density),
        af(G, sample.density), times=10
    )

    means <- aggregate(bm$time, by=list(bm$expr), FUN=mean)[,"x"]/1E6 # in milliseconds
    speed.loop[i] <- means[1]
    speed.dsample[i] <- means[2]
    speed.af[i] <- means[3]


}

speed.ylim <- range(c(speed.loop, speed.dsample, speed.af))
plot(lengths, (speed.loop), ylim=(speed.ylim), type="l", ylab="Time (milliseconds)", xlab="# Elements in G")
lines(lengths, (speed.dsample), col="red")
lines(lengths, (speed.af), col="blue")

有效地找到许多任意样本值的经验密度（）（如 dnorm（），但用于经验分布）

Efficiently find empirical density() for many arbitrary sample values (like dnorm(), but for empirical distribution)

r

kernel-density

probability-density

mapply

循环()

dsample()

af()

现在比较一下他们的速度：

现在比较随着 G 大小的增加速度：