эффективно найти эмпирическую плотность() для многих произвольных значений выборок (например, dnorm(), но для эмпирического распределения)

Question

Say вы определили эмпирическую плотность (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 распределен нормально с заданным средним и С.Д. , Я хотел бы использовать распределение, определенное эмпирически из произвольной выборки (что не обязательно является нормальным или любым другим хорошо описанным распределением, которое было бы в пакете статистики).

Answer 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()

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