0

行列 z (3 x 20000) があります。各行を確率変数と見なし、各列を 1 つのシミュレーションと見なします。Rで apply コマンドを使用して次の関数を作成し、3 次元の経験的累積分布関数 (EMP.CDF) を見つけました。この k 変量の経験的 CDF は、 この pdfの 2 ページ目の「多変量 ECDF」のセクションで説明されています。

EMP.CDF=function(z) {
# z is a matrix (3 x 20000) and each row is a realization of a random variable
q1=z[1,];q2=z[2,];q3=z[3,]
# qi = the realization of the ith random variable, i=1,2,3
# Now I am going to evaluate the empirical cumulative distribution function at
# each column of z
# Given each column, the function should return an empirical
# cumulative probability.

d=apply(z,2, function(x) sum(q1<=x[1] & q2<=x[2] & q3<=x[3])/(length(q1)))
return(d)}

> z=matrix(0,3,20000)
> z[1,]=runif(20000,1,2)
> z[2,]=runif(20000,3,5)
> z[3,]=runif(20000,7,9)

> system.time(EMP.CDF(z))
   user  system elapsed 
   30.18    0.01   30.39

上記のコードでは k=3 です。上記の関数をベクトル化してシステム時間を短縮する方法はありますか?

4

1 に答える 1

1

3 次元累積分布関数は、3 つの変数の関数です。グリッドで推定すると、3 次元配列として表すことができますが、不正確で巨大になります (関数は 1 次元配列を返すため、計算対象ではありません)。

ポイント が与えられたとき、xすべての座標が の座標よりも小さいポイントの割合を計算しxます。

z <- matrix(runif(60000), 3, 20000)
emp.cdf <- function(z)
  function(x) mean( apply( z <= x, 2, all ) )
emp.cdf(z)( c(.5,.5,.5) )  # Approximately 1/8

以下は、引用したドキュメントのプロットを再現しています。

n <- 10
z <- matrix(runif(2*n), 2, n)
f <- emp.cdf(z)
g <- function(u,v) f(c(u,v))
persp( outer( sort(z[1,]), sort(z[2,]), Vectorize(g) ) )

x <- seq(0,1,length=100)
persp( outer( x, x, Vectorize(g) ) )

初期点で累積確率分布を評価したい場合は、 を使用できますapply(グリッドで評価したい場合は、 を使用expand.gridして構築できます)。

n <- 100
z <- matrix(runif(3*n), 3, n)
f <- emp.cdf(z)
p <- apply( z, 2, f )

しかし、このアルゴリズムは 2 次です。n計算する確率があり、それらのそれぞれについて、すべての3*n座標を調べます。20,000 ポイントの場合、しばらく時間がかかります。

分割統治法を使用して計算を高速化できますが、簡単ではありません。ポイントをランダムに取得し、それを使用してスペースを 8 つの八分円に分割し、各八分円のポイント数を再帰的に計算します。次に、結果のツリーを使用 して、ポイントの一部のみを調べることで、任意のポイントでの確率を​​計算できます。

これは、 k-nearest neighborsの計算やn-body シミュレーションの高速化に使用される前処理ステップと同じです。

于 2013-03-05T09:08:04.197 に答える