0

假设您已经sample.density为以下示例定义了经验密度 ( ) x.sample：

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是具有指定均值和 sd 的正态分布，我想使用从任意样本凭经验确定的分布（不一定是正态的或任何其他经过充分描述的将在 stats 包中的分布）。

4

1 回答 1

1

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

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 <- 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))]
}

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

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")

在此处输入图像描述

于 2015-05-21T21:08:07.640 回答