I want to compute the P-value from the joint cumulative distribution of an n-dimensional order statistic.

One efficient way is using the following recursive formula.

However, the facts are (or would be):

1. I am too stupid to write a recursive function.
2. I didn’t find the efficient formula at first.
3. In other cases, the efficient formula have not been derived yet, or too complicated to derive.

In Statistics Monte Carlo simulation is a “quick” way to compute some complicated formulas. By saying “quick”, I mean I can see the results without knowing or deriving “ugly” Math formulas. It’s actually a very “slow” method in computing aspect.

Anyway, the R function is here.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

# sub function of monte carlo simulation to get the p-value
P_order_stat <- function(ranks) {
  NumRep <- 10000 # number of replicates
  newvec <- sort(ranks) # sort the rank ratios
  pvalue <- 0 # inital pvalue
  for (i in 1:NumRep){

      # generate random uniform distributed data,
      # and then sort the simulated rank ratios
      newx <- sort(runif(length(ranks), min=0, max=1))

      # if all the simulated data is lower than the input,
      # then sucess+1
      judge <- sum(newvec >= newx)
      if (judge == length(ranks)) pvalue<-pvalue+1
  }
  pvalue <- pvalue / NumRep  # get the p-value
  return(pvalue)
}