0

我有一个矩阵 z (3 x 20000)。将每一行视为一个随机变量，将每一列视为一个模拟。我使用 apply 命令在R中编写了以下函数，以查找 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 )

但是这个算法是二次的：有n概率要计算，对于每个概率，我们检查所有的3*n坐标。对于您的 20,000 点，这将需要一段时间。

您可以使用分而治之的方法来加快计算速度，但这并不简单：随机选取一个点，用它将空间分成 8 个八分圆，递归计算每个八分圆中的点数；然后，您可以使用生成的树 来计算任何点的概率，只需检查一小部分点。

这与用于计算k 最近邻或加速n 体模拟的预处理步骤没有什么不同。

于 2013-03-05T09:08:04.197 回答