0

I have a matrix z (3 x 20000). Consider each row as a random variable and each column as one simulation. I wrote the following function in R using apply command to find the empirical cumulative distribution function (EMP.CDF) in 3 dimensions. This k-variate empirical CDF was explained on the page 2 of this pdf, under the section of "The Multivariate 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

In above code k=3. Is there any way I can vectorize the above function to reduce the system time?

Stat
  • 671
  • 7
  • 19
  • 3
    Can you explain what you're trying to do (in your code)? Why are you comparing all column values with the first, second and third rows(`q1,q2,q3`)? It doesn't seem to make sense to me. – Arun Mar 05 '13 at 07:58
  • Sure, please check out the pdf referenced that I added to the my question. I am trying to calculate $\hat{F_{(X1,X2,X3)}(u1,u2,u3)}$. Here $X_i$ (i=1,2,3) is the ith random variable, i.e. the ith row of z or in my code, q1, q2 and q3. Now the function takes 3 numbers as $u1,u2,u3$. For each $u_i$ (i=1,2,3) it first look at the ith row and return a vector of TRUE and FALSE's. By using the & inside sum, I am taking intersection of these 3 Boolean vectors. By using "sum", it counts the number of TRUE's and then divide it by 20000, i.e the number of columns or # of simulations. – Stat Mar 05 '13 at 12:10

1 Answers1

1

A 3-dimensional cumulative distribution function is a function of 3 variables. If you estimate it on a grid, it could be represented as a 3-dimensional array, but it would be imprecise and huge (your function returns a 1-dimensional array, so it is not what it is computing).

Given a point x, just compute the proportion of points all of whose coordinates are less than those of 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

The following reproduces the plots in the document you cite:

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

If you want to evaluate the cumulative probability distribution on the initial points, you can just use apply (if you wanted to evaluate it on a grid, you could use expand.grid to build it).

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

But this algorithm is quadratic: there are n probabilities to compute, and for each of them, we examine all 3*n coordinates. For your 20,000 points, that will take a while.

You can use a divide-and-conquer approach to speed up the computations, but it is not straightforward: pick up a point at random, use it to separate the space into 8 octants, recursively compute the number of points in each octant; you can then use the resulting tree to compute the probability at any point, by examining only a fraction of the points.

This is not unlike the preprocessing step used to compute the k-nearest neighbours, or to speed up n-body simulations.

Vincent Zoonekynd
  • 31,893
  • 5
  • 69
  • 78
  • @ Vincent: my function is correct. Please check out the pdf I added (2nd page). Your function takes 3 numbers and then compares each number simultaneously with all the elements of the matrix z. This is not what I want. For the first argument (i.e. 0.5 in your code), it should just look at the first row and find how many of them are less than or equal to it. The same thing for 2nd and 3rd argument. Your function compares 0.5 within all 3 rows which is wrong. – Stat Mar 05 '13 at 12:20
  • Also, the function should return a vector and not a number, for each column it should provide me with a cumulative probability. – Stat Mar 05 '13 at 12:27
  • In the document you cite, the function takes k arguments (u1, u2,...,uk), and since it tries to estimate the probability `P(X1<=u1,...,Xk<=uk)`, it returns a number. In my code, since `x` is shorter than `z`, it will be recycled: its first element will be used on the first row, the second on the second, etc. I have added some code to reproduce the plots in the paper. – Vincent Zoonekynd Mar 05 '13 at 13:31
  • Thanks a lot. But I need something like these: g <- function(u,v,t) f(c(u,v,t)) and outer(z[1,], z[2,],z[3,], Vectorize(g)) Can you provide me with an example in 3 dim, I need it in 3 dim and outer is just for 2 dim. – Stat Mar 05 '13 at 14:44
  • Your code in 3-dim uses the "apply" function two times and hence taking so much time to compute. Thanks anyway. – Stat Mar 05 '13 at 16:23