Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

I have a distribution of observed measurements and I want to compare it to sampled distributions, using R. I have a program that samples distributions, according to a certain low of probability / constraint.

Suppose I have:

  • $m$ observed measurements
  • two sets of $n$ samples of $m$ measurements

I'd like to compute:

  • the average cumulative distribution function (a-cdf) of the first set of samples (its value at each evaluation points is the mean value of the cdf of each sample at this evaluation point, and the evaluation points are the pooled $n\cdot m$ measurements)
  • the 95% envelope of the first set of samples
  • the spatial distribution index (sdi) of my observed measurements using the second set of $n$-samples. To compute the sdi, we compute all the maximum differences between the cdf of each new sample and the a-cdf, and the maximum difference between the observed cdf and a-cdf. That makes $n+1$ points, and the sdi is the rank of the observed difference among the $n+1$ differences.

I'd like to know if thoses function exists in R, and if not, what would be the smartest way to implement them. Hope i'm clear enough.

share|improve this question
add comment

migrated from stackoverflow.com Mar 5 '12 at 14:50

This question came from our site for professional and enthusiast programmers.

1 Answer

up vote 2 down vote accepted

You can use the ecdf, quantile and rank functions to do this.

# Sample data
n <- 10
m <- 50
observed <- rnorm(m)
set1 <- matrix(rnorm(n*m),nc=n)
set2 <- matrix(rnorm(n*m),nc=n)

# Average cumulated distribution function of the first set of samples
ecdfs <- apply(set1, 2, ecdf)
acdf <- function(u) mean(sapply(ecdfs, function(f) f(u)))
acdf <- Vectorize(acdf)

# 95% envelope
alpha <- .025
x <- unique(sort(set1))
y <- sapply(ecdfs, function(f) f(x))
y <- apply(y, 1, quantile, c(alpha, 1-alpha))
lower <- approxfun( x, y[1,], method="constant", rule=2 )
upper <- approxfun( x, y[2,], method="constant", rule=2 )
curve(upper(x), xlim=c(-3,3))
curve(lower(x), add=TRUE)
curve(acdf(x), add=TRUE, lty=3)

# Rank of the observed sample 
# among the maximum differences
difference <- function(u) {
  max( abs( ecdf(u)(c(x,u)) - acdf(c(x,u)) ) )
}
sdi <- apply(set2, 2, difference)
sdi <- c( difference(observed), sdi )
sdi <- rank(sdi)[1]

# If you were to write a function to do this,
# it would return the following
list(
  acdf  = acdf,
  upper = upper,
  lower = lower, 
  sdi = sdi
)
share|improve this answer
    
Thank you Vincent, this is a great answer, very complete and useful. Would you know by the way if the Ripley's K function (computed from all couple of distances) has been implemented in R? Thanks again. Jeannot –  jeannot Mar 6 '12 at 8:44
    
I have a small question about your answer: why do you use 2 differents methods for computing the envoloppe and the acdf? would'nt it be equivalent to do : lower = function(u) quantile(sapply(ecdfs, function(f) f(u)), alpha) lower = Vectorize(lower) ? –  jeannot Mar 6 '12 at 9:21
    
There is no particular reason: pre-computing the quantiles may be slightly faster, but the two approaches should be equivalent. (After adding the missing method="constant", they are now equivalent.) –  Vincent Zoonekynd Mar 6 '12 at 9:39
    
Hi Vincent, could you explain me your "difference" function please? what's the "x" variable? thanks –  jeannot Mar 22 '12 at 13:47
    
x contains the data from set1. We want to compute the maximum difference between ecdf(u) and acdf (these are functions, not vectors). Since these are stepwise constant functions, it is sufficient to evaluate them at each jump: x contains the jump points of acdf and u those of ecdf(u), so I evaluate both functions at c(x,u). –  Vincent Zoonekynd Mar 22 '12 at 22:40
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.