I will use an example to explain my problem, but I think that it can be generalised to any randomization testing that does not allow one to create the null distribution by enumerating of all the possible orders of the data (that is the number of possible data permutation is too large to make the complete enumeration impossible), and a sampling of the randomised data should be made.

Example: Differential Correlation

Let be $\mathbf x_1$ and $\mathbf y_1$, both of length $n_1$, a pair of vectors that represents the age and weight of $n_1$ male individuals, and $\mathbf x_2$ and $\mathbf y_2$, both of length $n_2$, a pair of vectors that represents the age and weight of $n_2$ female individuals. The number of males may be different than the number of females.

Let $\rho_1$ be the Pearson's correlation coefficient between $\mathbf x_1$ and $\mathbf y_1$, and $\rho_2$ the Pearson's correlation coefficient between $\mathbf x_2$ and $\mathbf y_2$. To reject the null hypothesis that $\rho_1=\rho_2$ using a permutation test I should create random vectors where paired elements of $\mathbf x_1/\mathbf y_1$ are swapped with paired elements of $\mathbf x_2/\mathbf y_2$, in order to remove the relationship between the measurements and the labels (male/female).

The complete enumeration of all possible order of the data would be equal to all the possible order of the vector $\mathbf x (\mathbf y)$ built by concatenating the two vectors $\mathbf x_1 \mathbf x_2 (\mathbf y_1 \mathbf y_2)$, that is $(n_1+n_2)$! This makes my complete enumeration impossible, and I need to sample $N$ configurations from all the possible configuration of $\mathbf x$, of length $n_1+n_2$.

The problem: How many elements are swapped (in average) in each data configuration?

Each possible configuration will swap a random number $R$ of (paired) elements from $\mathbf x_1/\mathbf y_1$ to $\mathbf x_2/\mathbf y_2$, and I would like to know if $R$ follows any known distribution.

Let then $\mathbf v$ be a vector of length $n_1+n_2$, built by concatenating a vector $\mathbf v_1$, that includes $n_1$ elements set to $0$ representing $\mathbf x_1/\mathbf y_1$, and $\mathbf v_2$ that includes $n_2$ elements set to $1$ representing $\mathbf x_2/\mathbf y_2$. Let then shuffle the element of $\mathbf v$ for $N$ times, and count, at each step $i$ (with $i \in [i, N]$), how many elements $c_i$ of $n_1$ are now set to $1$, that is, how many elements were swapped at that permutation step. Let then build the distribution of my permuted elements $c$.

These steps are performed by the following function, while the histogram will show the distribution:

count.permuted.elements <- function(n1, n2, N=10000) {
    c <- numeric(N)
    for (i in 1:N) {   
        v <- c(rep(0, n1), rep(1, n2))
        v <- sample(v)
        c[i] <- sum(v[1:n1])


hist(count.permuted.elements(n1, n2))

I have run some simulations with different values of $n_1/n_2$, and what I obtained (see plots below) are normal distributions, as also confirmed by the function fitdist (in the R package fitdistrplus). Let me underline that the right tail is limited by the length of the shortest vector, that is, I can't permute more elements than those included in the shorter vector (see the plot on the bottom right for an example).

enter image description here

The problem (restated): How to estimate the parameters of the normal distributions ($\mu$ and $\sigma$) showed above, by using only $n_1$ and $n_2$?

After running several simulation, varying $n_1/n_2$, it seems to me that the mean $\mu$ is a function of the length of the smaller vector (let's say $n_1$) and the square root of the ratio between the two lengths, that is:

$$\mu = n_1 \cdot \left(-0.604 \cdot \sqrt{\frac{n_1}{n_2}} + 1.101\right)$$

However, I can't find a way to calculate $\sigma$ --and I am not sure that what I did to evaluate $\mu$ is very correct.

My question is: does anyone know how to evaluate, in a formal or empirical way, the parameters of the normal distributions ($\mu$ and $\sigma$) build by counting how many elements are swapped (in average) between the two arrays in each data configuration?

Thank you very much for your help!

  • 1
    $\begingroup$ You have a fundamental problem: the permutation test you describe does not work the way you believe. It does not test the hypothesis that the two correlation coefficients are equal. $\endgroup$
    – whuber
    Mar 2, 2016 at 19:18
  • 1
    $\begingroup$ You don't need to know how many are swapped on average to do the randomization; indeed if you're using such a calculation in order to randomize (by swapping that average number) you'll be doing randomization wrong. $\endgroup$
    – Glen_b
    Mar 3, 2016 at 1:21
  • $\begingroup$ @whuber Which part of the question are you referring to? The code (and the main question in the post) is not meant to test if two correlation coefficients are equal, but to simulate how many elements are swapped at each permutation step. $\endgroup$
    – alesssia
    Mar 3, 2016 at 8:20
  • $\begingroup$ @Glen_b I think I got what you mean. But how to implement a correct permutation scheme then? For a number of reasons I can't concatenate my two vectors, as done in the example, and I need to decide how many elements to sample (permute) at each step. $\endgroup$
    – alesssia
    Mar 3, 2016 at 8:22
  • $\begingroup$ Indeed, I think "how to implement a correct permutation scheme?" is what you really need to ask in your question; at present you seem to have an XY-problem $\endgroup$
    – Glen_b
    Mar 3, 2016 at 11:27

1 Answer 1


You have $n_1$ $0$'s and $n_2$ $1$'s, so probability of sampling $1$'s from the whole sample is $p = \frac{n_2}{n_1+n_2}$. You randomly sample from it $n_1$ cases, where probability of sampling $1$'s is $p$, when $n_1 \ll n_2$ this resembles binomial distribution with parameters $n_1$ and $p$.

However your data does not follow binomial distribution since you are sampling without replacement from finite population. Your data follows a hypergeometric distribution. In this case mean and variance are

$$ \mu = n_1 \frac{n_2}{n_1+n_2} \\ \sigma^2 = n_1 \frac{n_2}{n_1+n_2} \frac{n_1}{n_1+n_2} \frac{n_2}{n_1+n_2-1} $$

Below you can see data created with your function with added pmf of binomial (red curve) and hypergeometric (blue curve) distribution.

enter image description here

  • $\begingroup$ Thank you very much. This quite good approximation is more than fine! $\endgroup$
    – alesssia
    Mar 2, 2016 at 12:00

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.