使用“应用”加速模拟

Question

我有一个矩阵 z (3 x 20000)。将每一行视为一个随机变量，将每一列视为一个模拟。我在R 中编写了以下函数，使用apply 命令查找3 维的经验累积分布函数（EMP.CDF）。这个 k 变量经验 CDF 在 this 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。有什么办法可以对上述函数进行矢量化以减少系统时间？

Answer 1

3 维累积分布函数是 3 个变量的函数。 如果你在一个网格上估计它，它可以表示为一个 3 维数组， 但这将是不精确且巨大的（您的函数返回一个一维数组， 所以它不是它正在计算的东西）。

给定一个点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个八分圆， 递归计算每个八分圆中的点数； 然后您可以使用生成的tree 计算任意点的概率， 通过只检查一小部分点。

这与预处理步骤没有什么不同 用于计算k-nearest neighbours， 或加速n-body simulations。