# How to generate a large full-rank random correlation matrix with some strong correlations present?

I would like to generate a random correlation matrix $\mathbf C$ of $n \times n$ size such that there are some moderately strong correlations present:

• square real symmetric matrix of $n \times n$ size, with e.g. $n=100$;
• positive-definite, i.e. with all eigenvalues real and positive;
• full rank;
• all diagonal elements equal to $1$;
• off-diagonal elements should be reasonably uniformly distributed on $(-1, 1)$. Exact distribution does not matter, but I would like to have some moderately large amount (e.g. $10\%$) of moderately large values (e.g. with absolute value of $0.5$ or higher). Basically I want to make sure that $\mathbf C$ is not almost diagonal with all off-diagonal elements $\approx 0$.

Is there a simple way to do it?

The purpose is to use such random matrices to benchmark some algorithms working with correlation (or covariance) matrices.

## Methods that do not work

Here are some ways to generate random correlation matrices that I know of, but that do not work for me here:

1. Generate random $\mathbf X$ of $s \times n$ size, center, standardize and form the correlation matrix $\mathbf C=\frac{1}{s-1}\mathbf X^\top \mathbf X$. If $s>n$, this will generally result in all off-diagonal correlations being around $0$. If $s\ll n$, some correlations will be strong, but $\mathbf C$ will not be full rank.

2. Generate random positive definite matrix $\mathbf B$ in one of the following ways:

• Generate random square $\mathbf A$ and make symmetric positive definite $\mathbf B = \mathbf A \mathbf A^\top$.

• Generate random square $\mathbf A$, make symmetric $\mathbf E = \mathbf A + \mathbf A^\top$, and make it positive definite by performing eigen-decomposition $\mathbf E = \mathbf U \mathbf S \mathbf U^\top$ and setting all negative eigenvalues to zero: $\mathbf B = \mathbf U \:\mathrm{max}\{\mathbf S, \mathbf 0\} \:\mathbf U^\top$. NB: this will result in a rank-deficient matrix.

• Generate random orthogonal $\mathbf Q$ (e.g. by generating random square $\mathbf A$ and doing its QR decomposition, or via Gram-Schmidt process) and random diagonal $\mathbf D$ with all positive elements; form $\mathbf B = \mathbf Q \mathbf D \mathbf Q^\top$.

Obtained matrix $\mathbf B$ can be easily normalized to have all ones on the diagonal: $\mathbf C = \mathbf D^{-1/2}\mathbf B \mathbf D^{-1/2}$, where $\mathbf D = \mathrm{diag}\:\mathbf B$ is the diagonal matrix with the same diagonal as $\mathbf B$. All three ways listed above to generate $\mathbf B$ result in $\mathbf C$ having off-diagonal elements close $0$.

After posting my question, I found two almost duplicates in the past:

Unfortunately, none of these threads contained a satisfactory answer (until now :)

• You could create random orthogonal matrix by QR or Gram-Schmidt processes. That will be "eigenvectors of PCA". Add scale to its columns (turn into "loadings"). Get the covariance matrix from these loadings. Something like that... – ttnphns Nov 18 '14 at 16:20
• Uhm, well.. Imagine we want to create a nXk loading matrix W, not fully random but the one we want (it will, WW'+diag(noise), define the cov matrix we seek. The only task is to correct the column-normalized W (i.e. the k "eigenvectors") to become orthogonal. Any method to de-correlate correlated variables (here variables are the eigenvectors) will probably do. (This a raw idea.) – ttnphns Nov 18 '14 at 17:08
• Ah, @whuber, now I see what you mean. Yes, you are right: if all off-diagonal elements are identical and equal to $\rho$, then the matrix is indeed full rank and positive-definite... This is of course not what I had in mind: I would like the distribution of off-diagonal elements in each matrix to be reasonably "spread", not the distribution across matrices... – amoeba Nov 18 '14 at 20:12
• You might want to look into the LKJ distribution – shadowtalker Nov 18 '14 at 22:07
• @ttnphns: I think I have finally understood that you were right all along: what you suggested is the simplest way to get to the goal. I added an update to my answer implementing essentially what you wrote above. – amoeba Nov 23 '14 at 23:53

Other answers came up with nice tricks to solve my problem in various ways. However, I found a principled approach that I think has a large advantage of being conceptually very clear and easy to adjust.

In this thread: How to efficiently generate random positive-semidefinite correlation matrices? -- I described and provided the code for two efficient algorithms of generating random correlation matrices. Both come from a paper by Lewandowski, Kurowicka, and Joe (2009), that @ssdecontrol referred to in the comments above (thanks a lot!).

Please see my answer there for a lot of figures, explanations, and matlab code. The so called "vine" method allows to generate random correlation matrices with any distribution of partial correlations and can be used to generate correlation matrices with large off-diagonal values. Here is the example figure from that thread:

The only thing that changes between subplots, is one parameter that controls how much the distribution of partial correlations is concentrated around $\pm 1$.

I copy my code to generate these matrices here as well, to show that it is not longer than the other methods suggested here. Please see my linked answer for some explanations. The values of betaparam for the figure above were ${50,20,10,5,2,1}$ (and dimensionality d was $100$).

function S = vineBeta(d, betaparam)
P = zeros(d);           %// storing partial correlations
S = eye(d);

for k = 1:d-1
for i = k+1:d
P(k,i) = betarnd(betaparam,betaparam); %// sampling from beta
P(k,i) = (P(k,i)-0.5)*2;     %// linearly shifting to [-1, 1]
p = P(k,i);
for l = (k-1):-1:1 %// converting partial correlation to raw correlation
p = p * sqrt((1-P(l,i)^2)*(1-P(l,k)^2)) + P(l,i)*P(l,k);
end
S(k,i) = p;
S(i,k) = p;
end
end

%// permuting the variables to make the distribution permutation-invariant
permutation = randperm(d);
S = S(permutation, permutation);
end


Update: eigenvalues

@psarka asks about the eigenvalues of these matrices. On the figure below I plot the eigenvalue spectra of the same six correlation matrices as above. Notice that they decrease gradually; in contrast, the method suggested by @psarka generally results in a correlation matrix with one large eigenvalue, but the rest being pretty uniform.

## Update. Really simple method: several factors

Similar to what @ttnphns wrote in the comments above and @GottfriedHelms in his answer, one very simple way to achieve my goal is to randomly generate several ($k<n$) factor loadings $\mathbf W$ (random matrix of $k \times n$ size), form the covariance matrix $\mathbf W \mathbf W^\top$ (which of course will not be full rank) and add to it a random diagonal matrix $\mathbf D$ with positive elements to make $\mathbf B = \mathbf W \mathbf W^\top + \mathbf D$ full rank. The resulting covariance matrix can be normalized to become a correlation matrix (as described in my question). This is very simple and does the trick. Here are some example correlation matrices for $k={100, 50, 20, 10, 5, 1}$:

The only downside is that the resulting matrix will have $k$ large eigenvalues and then a sudden drop, as opposed to a nice decay shown above with the vine method. Here are the corresponding spectra:

Here is the code:

d = 100;    %// number of dimensions
k = 5;      %// number of factors

W = randn(d,k);
S = W*W' + diag(rand(1,d));
S = diag(1./sqrt(diag(S))) * S * diag(1./sqrt(diag(S)));

• +1. However, here's just a reminder to your last section on "factor method". Strictly right approach calls that columns of W are orthogonal (i.e. cosines between them are 0). Simply generating random W of course doesn't provide it. If they aren't orthogonal - i.e. factors are oblique (call then W as W_) - factor theorem is not WW' but W_CW_' with C being "correlations" (cosines) between factors. Now, C=Q'Q with Q being the nonorthogonal rotation matrix of rotation W_=inv(Q)'W (and so W=W_Q'). Generate some Q - a matrix with column ss = 1 and matrix ss = size of matrix. – ttnphns Nov 24 '14 at 8:24
• ...typo: not W_=inv(Q)'W, of course W_= W inv(Q)'. – ttnphns Nov 24 '14 at 8:39
• @ttnphns: What you are saying is correct, but I don't think it matters for the purposes of generating random correlation matrices. If I generate $W$ randomly, then yes, its columns will not be exactly orthogonal, but $WW^\top+D$ (after normalization to get all ones on the diagonal) will still be full rank correlation matrix with some large off-diagonal values, which is what the question was about. Of course it is very easy to orthogonalize columns of $W$ beforehand, I just did not see why this would be necessary in this case. – amoeba Nov 24 '14 at 9:55
• Translating this into R: W = replicate(k, rnorm(d)); S = W%*%t(W) + diag(rnorm(d),nrow=d); S = diag(1/sqrt(diag(S)))%*%S%*%diag(1/sqrt(diag(S))) – Scott Worland Aug 22 '18 at 20:05
• @Mihai, good point and your suggestions are likely the simplest. You could also do S <- matrix(nearPD(S, corr = TRUE, keepDiag = TRUE)$mat@x,ncol(S),ncol(S)) – Scott Worland Apr 8 '19 at 14:07 A simple thing but maybe will work for benchmark purposes: took your 2. and injected some correlations into starting matrix. Distribution is somewhat uniform, and changing$a$you can get concentration near 1 and -1 or near 0. import numpy as np from random import choice import matplotlib.pyplot as plt n = 100 a = 2 A = np.matrix([np.random.randn(n) + np.random.randn(1)*a for i in range(n)]) A = A*np.transpose(A) D_half = np.diag(np.diag(A)**(-0.5)) C = D_half*A*D_half vals = list(np.array(C.ravel())[0]) plt.hist(vals, range=(-1,1)) plt.show() plt.imshow(C, interpolation=None) plt.show()  • (+1) Thanks! I edited your answer to add the prettify coloring scheme for python and make it python2 compatible :) Hope it's okay. I also deleted my above comments to remove clutter (you can delete yours as well). I am now trying to understand the logic of your code; do you really need crs array? What influence does$k$have? My feeling is that you can simply add a random number from$[-a, a]$to each row, no? This is similar to using matrix$X$with correlated samples (not features) and then computing its sample correlation matrix, correct? – amoeba Nov 19 '14 at 12:34 • Yes, you are completely right! (Oh boy, that was indeed silly :D). I changed the random part to randn(1)*a and now it's much better. – psarka Nov 19 '14 at 12:58 • Thanks! Now you only need to remove your mentioning of$k$above the code. Meanwhile, I found an old question that is nearly a duplicate, and posted an answer there, directly inspired by yours here: How to generate random correlation matrix that has approximately normally distributed off-diagonal entries with given standard deviation?. Seems to work nice! – amoeba Nov 19 '14 at 13:12 • Cool! I wonder if such procedure results in a known distribution. Experimenting with different$a$and$n$I can get quite a few different shapes. – psarka Nov 19 '14 at 13:31 • One downside of this method is that the resulting correlation matrix has one large eigenvalue, but the remaining ones are nearly uniform. So this procedure does not yield a "general" correlation matrix... Not that I specified it in my question. But @ssdecontrol mentioned in the comments above that apparently there are ways to sample from all correlation matrices; this looks interesting but much more complicated. – amoeba Nov 19 '14 at 13:35 Hmm, after I' done an example in my MatMate-language I see that there is already a python-answer, which might be preferable because python is widely used. But because you had still questions I show you my approach using the Matmate-matrix-language, perhaps it is more selfcommenting. Method 1 (Using MatMate): v=12 // 12 variables f=3 // subset-correlation based on 3 common factors vg = v / f // variables per subsets // generate hidden factor-matrix // randomu(rows,cols ,lowbound, ubound) gives uniform random matrix // without explicite bounds the default is: randomu(rows,cols,0,100) L = { randomu(vg,f) || randomu(vg,f)/100 || randomu(vg,f)/100 , _ randomu(vg,f)/100 || randomu(vg,f) || randomu(vg,f)/100 , _ randomu(vg,f)/100 || randomu(vg,f)/100 || randomu(vg,f) } // make sure there is itemspecific variance // by appending a diagonal-matrix with random positive entries L = L || mkdiag(randomu(v,1,10,20)) // make covariance and correlation matrix cov = L *' // L multiplied with its transpose cor = covtocorr(cov) set ccdezweite=3 ccfeldweite=8 list cor cor = 1.000, 0.321, 0.919, 0.489, 0.025, 0.019, 0.019, 0.030, 0.025, 0.017, 0.014, 0.014 0.321, 1.000, 0.540, 0.923, 0.016, 0.015, 0.012, 0.030, 0.033, 0.016, 0.012, 0.015 0.919, 0.540, 1.000, 0.679, 0.018, 0.014, 0.012, 0.029, 0.028, 0.014, 0.012, 0.012 0.489, 0.923, 0.679, 1.000, 0.025, 0.022, 0.020, 0.040, 0.031, 0.014, 0.011, 0.014 0.025, 0.016, 0.018, 0.025, 1.000, 0.815, 0.909, 0.758, 0.038, 0.012, 0.018, 0.014 0.019, 0.015, 0.014, 0.022, 0.815, 1.000, 0.943, 0.884, 0.035, 0.012, 0.014, 0.012 0.019, 0.012, 0.012, 0.020, 0.909, 0.943, 1.000, 0.831, 0.036, 0.013, 0.015, 0.010 0.030, 0.030, 0.029, 0.040, 0.758, 0.884, 0.831, 1.000, 0.041, 0.017, 0.022, 0.020 0.025, 0.033, 0.028, 0.031, 0.038, 0.035, 0.036, 0.041, 1.000, 0.831, 0.868, 0.780 0.017, 0.016, 0.014, 0.014, 0.012, 0.012, 0.013, 0.017, 0.831, 1.000, 0.876, 0.848 0.014, 0.012, 0.012, 0.011, 0.018, 0.014, 0.015, 0.022, 0.868, 0.876, 1.000, 0.904 0.014, 0.015, 0.012, 0.014, 0.014, 0.012, 0.010, 0.020, 0.780, 0.848, 0.904, 1.000  The problem here might be, that we define blocks of submatrices which have high correlations within with little correlation between and this is not programmatically but by the constant concatenation-expressions . Maybe this approach could be modeled more elegantly in python. Method 2(a) After that, there is a completely different approach, where we fill the possible remaining covariance by random amounts of 100 percent into a factor-loadings-matrix. This is done in Pari/GP: {L = matrix(8,8); \\ generate an empty factor-loadings-matrix for(r=1,8, rv=1.0; \\ remaining variance for variable is 1.0 for(c=1,8, pv=if(c<8,random(100)/100.0,1.0); \\ define randomly part of remaining variance cv= pv * rv; \\ compute current partial variance rv = rv - cv; \\ compute the now remaining variance sg = (-1)^(random(100) % 2) ; \\ also introduce randomly +- signs L[r,c] = sg*sqrt(cv) ; \\ compute factor loading as signed sqrt of cv ) );} cor = L * L~  and the produced correlation-matrix is  1.000 -0.7111 -0.08648 -0.7806 0.8394 -0.7674 0.6812 0.2765 -0.7111 1.000 0.06073 0.7485 -0.7550 0.8052 -0.8273 0.05863 -0.08648 0.06073 1.000 0.5146 -0.1614 0.1459 -0.4760 -0.01800 -0.7806 0.7485 0.5146 1.000 -0.8274 0.7644 -0.9373 -0.06388 0.8394 -0.7550 -0.1614 -0.8274 1.000 -0.5823 0.8065 -0.1929 -0.7674 0.8052 0.1459 0.7644 -0.5823 1.000 -0.7261 -0.4822 0.6812 -0.8273 -0.4760 -0.9373 0.8065 -0.7261 1.000 -0.1526 0.2765 0.05863 -0.01800 -0.06388 -0.1929 -0.4822 -0.1526 1.000  Possibly this generates a correlation-matrix with dominant principal components because of the cumulative generating-rule for the factor-loadings-matrix. Also it might be better to assure positive definiteness by making the last portion of variance a unique factor. I left it in the program to keep the focus on the general principle. A 100x100 correlation-matrix had the following frequencies of correlations (rounded to 1 dec place)  e f e: entry(rounded) f: frequency ----------------------------------------------------- -1.000, 108.000 -0.900, 460.000 -0.800, 582.000 -0.700, 604.000 -0.600, 548.000 -0.500, 540.000 -0.400, 506.000 -0.300, 482.000 -0.200, 488.000 -0.100, 464.000 0.000, 434.000 0.100, 486.000 0.200, 454.000 0.300, 468.000 0.400, 462.000 0.500, 618.000 0.600, 556.000 0.700, 586.000 0.800, 536.000 0.900, 420.000 1.000, 198.000  [update]. Hmm, the 100x100 matrix is badly conditioned; Pari/GP cannot determine the eigenvalues correctly with the polroots(charpoly())-function even with 200 digits precision. I've done a Jacobi-rotation to pca-form on the loadingsmatrix L and find mostly extremely small eigenvalues, printed them in logarithms to base 10 (which give roughly the position of the decimal point). Read from left to right and then row by row: log_10(eigenvalues): 1.684, 1.444, 1.029, 0.818, 0.455, 0.241, 0.117, -0.423, -0.664, -1.040 -1.647, -1.799, -1.959, -2.298, -2.729, -3.059, -3.497, -3.833, -4.014, -4.467 -4.992, -5.396, -5.511, -6.366, -6.615, -6.834, -7.535, -8.138, -8.263, -8.766 -9.082, -9.482, -9.940, -10.167, -10.566, -11.110, -11.434, -11.788, -12.079, -12.722 -13.122, -13.322, -13.444, -13.933, -14.390, -14.614, -15.070, -15.334, -15.904, -16.278 -16.396, -16.708, -17.022, -17.746, -18.090, -18.358, -18.617, -18.903, -19.186, -19.476 -19.661, -19.764, -20.342, -20.648, -20.805, -20.922, -21.394, -21.740, -21.991, -22.291 -22.792, -23.184, -23.680, -24.100, -24.222, -24.631, -24.979, -25.161, -25.282, -26.211 -27.181, -27.626, -27.861, -28.054, -28.266, -28.369, -29.074, -29.329, -29.539, -29.689 -30.216, -30.784, -31.269, -31.760, -32.218, -32.446, -32.785, -33.003, -33.448, -34.318  [update 2] Method 2(b) An improvement might be to increase the itemspecific variance to some non-marginal level and reduce to a reasonably smaller number of common factors (for instance integer-squareroot of itemnumber): { dimr = 100; dimc = sqrtint(dimr); \\ 10 common factors L = matrix(dimr,dimr+dimc); \\ loadings matrix \\ with dimr itemspecific and \\ dimc common factors for(r=1,dim, vr=1.0; \\ complete variance per item vu=0.05+random(100)/1000.0; \\ random variance +0.05 \\ for itemspecific variance L[r,r]=sqrt(vu); \\ itemspecific factor loading vr=vr-vu; for(c=1,dimc, cv=if(c<dimc,random(100)/100,1.0)*vr; vr=vr-cv; L[r,dimr+c]=(-1)^(random(100) % 2)*sqrt(cv) ) );} cov=L*L~ cp=charpoly(cov) \\ does not work even with 200 digits precision pr=polroots(cp) \\ spurious negative and complex eigenvalues...  The structure of the result in term of the distribution of correlations: remains similar (also the nasty non decomposability by PariGP), but the eigenvalues, when found by jacobi-rotation of the loadingsmatrix, have now a better structure, for a newly computed example I got the eigenvalues as log_10(eigenvalues): 1.677, 1.326, 1.063, 0.754, 0.415, 0.116, -0.262, -0.516, -0.587, -0.783 -0.835, -0.844, -0.851, -0.854, -0.858, -0.862, -0.862, -0.868, -0.872, -0.873 -0.878, -0.882, -0.884, -0.890, -0.895, -0.896, -0.896, -0.898, -0.902, -0.904 -0.904, -0.909, -0.911, -0.914, -0.920, -0.923, -0.925, -0.927, -0.931, -0.935 -0.939, -0.939, -0.943, -0.948, -0.951, -0.955, -0.956, -0.960, -0.967, -0.969 -0.973, -0.981, -0.986, -0.989, -0.997, -1.003, -1.005, -1.011, -1.014, -1.019 -1.022, -1.024, -1.031, -1.038, -1.040, -1.048, -1.051, -1.061, -1.064, -1.068 -1.070, -1.074, -1.092, -1.092, -1.108, -1.113, -1.120, -1.134, -1.139, -1.147 -1.150, -1.155, -1.158, -1.166, -1.171, -1.175, -1.184, -1.184, -1.192, -1.196 -1.200, -1.220, -1.237, -1.245, -1.252, -1.262, -1.269, -1.282, -1.287, -1.290  • Thanks a lot! Very interesting, but will take me some time to digest... – amoeba Nov 21 '14 at 10:11 • I still have to go carefully through your answer, but in the meantime I read a paper about sampling random correlation matrices, and one of the methods from there can be used to do exactly what I need. I posted an answer here, you might be interested to take a look! It links to a much more detailed answer I wrote in another thread. – amoeba Nov 21 '14 at 23:39 • @amoeba: happy you found something nicely working for you! It's an interesting question, I'll come back later to this myself, perhaps improve/adapt the MatMate-procedures (and make them subroutines) according to the paper that you've worked on. – Gottfried Helms Nov 24 '14 at 5:30 Interesting question (as always!). How about finding a set of example matrices that exhibit the properties you desire, and then take convex combinations thereof, since if$A$and$B$are positive definite, then so is$\lambda A + (1-\lambda)B$. As a bonus, no rescaling of the diagonals will be necessary, by the convexity of the operation. By adjusting the$\lambda$to being more concentrated towards 0 and 1 versus uniformly distributed, you could concentrate the samples on the edges of the polytope, or the interior. (You could use a beta/Dirichlet distribution to control the concentration vs uniformity). For example, you could let$A$to be component-symmetric, and$B$be toeplitz. Of course, you can always add another class$C$, and take$\lambda_A A + \lambda_B B + \lambda_C C$such that$\sum \lambda = 1$and$\lambda \geq 0$, and so on. • Thanks for the suggestion, Andrew, but of course it would be nicer to have an unbiased method that does not need to start with some predefined$A$and$B\$... In the comments to my original question @ssdecontrol referred to a paper describing algorithms to sample correlation matrices uniformly (in a certain precise sense), or biased towards identity matrix, but I cannot find a way yet to sample them biased away from identity... I have also found a couple of old threads here asking almost the same question, maybe you will be interested, see my update. – amoeba Nov 21 '14 at 0:09
• Ah, but from such an algorithm, and a suitable diversity in the "vertices" (that is, matrices) that define your polytope of positive-definite correlation matrices, you can use rejection sampling to get whatever distribution of eigenvalues, uniformity of entries, etc, that you desire. However, it's not clear to me what a good basis would be. Sounds like a question for someone who has studied abstract algebra more recently than I. – Andrew M Nov 21 '14 at 2:56
• Hi again, I read a paper about sampling random correlation matrices, and one of the methods from there can be used to do exactly what I need. I posted an answer here, you might be interested to take a look! It links to a much more detailed answer I wrote in another thread. – amoeba Nov 21 '14 at 23:38

R has a package (clusterGeneration) that implements the method in:

Example:

> (cormat10 = clusterGeneration::rcorrmatrix(10, alphad = 1/100000000000000))
[,1]   [,2]    [,3]     [,4]     [,5]   [,6]   [,7]    [,8]     [,9]   [,10]
[1,]  1.000  0.344 -0.1406 -0.65786 -0.19411  0.246  0.688 -0.6146  0.36971 -0.1052
[2,]  0.344  1.000 -0.4256 -0.35512  0.15973  0.192  0.340 -0.4907 -0.30539 -0.6104
[3,] -0.141 -0.426  1.0000  0.01775 -0.61507 -0.485 -0.273  0.3492 -0.30284  0.1647
[4,] -0.658 -0.355  0.0178  1.00000  0.00528 -0.335 -0.124  0.5256 -0.00583 -0.0737
[5,] -0.194  0.160 -0.6151  0.00528  1.00000  0.273 -0.350 -0.0785  0.08285  0.0985
[6,]  0.246  0.192 -0.4847 -0.33531  0.27342  1.000  0.278 -0.2220 -0.11010  0.0720
[7,]  0.688  0.340 -0.2734 -0.12363 -0.34972  0.278  1.000 -0.6409  0.40314 -0.2800
[8,] -0.615 -0.491  0.3492  0.52557 -0.07852 -0.222 -0.641  1.0000 -0.50796  0.1461
[9,]  0.370 -0.305 -0.3028 -0.00583  0.08285 -0.110  0.403 -0.5080  1.00000  0.3219
[10,] -0.105 -0.610  0.1647 -0.07373  0.09847  0.072 -0.280  0.1461  0.32185  1.0000
> cormat10[lower.tri(cormat10)] %>% psych::describe()
vars  n  mean   sd median trimmed mad   min  max range skew kurtosis   se
X1    1 45 -0.07 0.35  -0.08   -0.07 0.4 -0.66 0.69  1.35 0.03       -1 0.05


Unfortunately, it doesn't seem possible to simulate correlations that follow a uniform-ish distribution with this. It seems make stronger correlations when alphad is set to very small values, but even at 1/100000000000000, the range of correlations would only go up to about 1.40.

Nonetheless, I hope this might be of some use to someone.