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、任意のサンプルから経験的に決定された分布を使用したいと思います統計パッケージに含まれるディストリビューション)。

4

1 に答える 1

1

@MrFlick のコメントは、正しい方向を示してくれたと思います。提案されたapproxfunアプローチとループ アプローチの例に加えて、forを使用できることにも気付きましたmapply。approxfun私の の使用は、を使用する他の 2 つのアプローチによる結果と正確に一致しないことに注意してください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()

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 に答える