3
$\begingroup$

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.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
3
$\begingroup$

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
)
Improve this answer
$\endgroup$
5
  • $\begingroup$ 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 $\endgroup$
    – jeannot
    Commented Mar 6, 2012 at 8:44
  • $\begingroup$ 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) ? $\endgroup$
    – jeannot
    Commented Mar 6, 2012 at 9:21
  • $\begingroup$ 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.) $\endgroup$
    – Vincent Zoonekynd
    Commented Mar 6, 2012 at 9:39
  • $\begingroup$ Hi Vincent, could you explain me your "difference" function please? what's the "x" variable? thanks $\endgroup$
    – jeannot
    Commented Mar 22, 2012 at 13:47
  • $\begingroup$ 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). $\endgroup$
    – Vincent Zoonekynd
    Commented Mar 22, 2012 at 22:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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