# simulate rank correlations of random vectors of length 9
sim9ranks <- sapply(1:10000, function(i) {
  x <- rnorm(9)
  y <- rnorm(9)
  cor(rank(x), rank(y))
})

# preliminaries for plotting
delta <- min(diff(sort(unique(sim9ranks))))
hist.breaks <- seq(min(sim9ranks)-delta/4, max(sim9ranks)+delta, by=delta)
xx <- seq(-1, 1, by=0.01)

# show that the beta distribution fits the observations
# Note that we have to shift and scale to move from the
# interval [0,1] to the interval [-1, 1]
hist(sim9ranks, prob=TRUE, breaks=hist.breaks)
lines(xx, dbeta((xx+1)/2, 7/2, 7/2)/2)

# function to generate correlated multivariate normal data
rmvn <- function(n, mu, sigma) {
	y <- chol(sigma)
	# yields a matrix that satisfies t(y) %*% y = sigma
	m <- length(mu)
	spread <- t(matrix(rnorm(m*n), nrow=m, ncol=n)) %*% y	
	sweep(spread, 2, mu, '+')
}

# simulated rank correlations when the vectors are correlated
sim9relatedranks <- sapply(1:1000, function(i) {
  rho <- runif(1, 0.3, 0.9) # choose a correlation coefficient
  sigma <- matrix(rho, 10, 10) + diag(rep(1-rho, 10)) # covariance matrix
  x <- rmvn(9, rep(0,10), sigma) # set of ten correlated vectors
  y <- rnorm(9)
  sapply(1:10, function(i) {cor(rank(x[,i]), rank(y))})
})
sim9rel <- as.vector(sim9relatedranks) # convert from matrix to vector

# preliminaries for plotting
delta <- min(diff(sort(unique(sim9rel))))
hist.breaks <- seq(min(sim9rel)-delta/4, max(sim9rel)+delta, by=delta)
xx <- seq(-1, 1, by=0.01)

# show that the beta distribution still fits the data
hist(as.vector(sim9rel), prob=TRUE, breaks=hist.breaks)
lines(xx, dbeta((xx+1)/2, 7/2, 7/2)/2)