¿Cómo generar una gran matriz de correlación aleatoria de rango completo con algunas correlaciones fuertes presentes?

Me gustaría generar una matriz de correlación aleatoria $\mathbf C$ de tamaño $n \times n$ modo que haya algunas correlaciones moderadamente fuertes presentes:

• matriz simétrica real cuadrada de tamaño $n \times n$ , con, por ejemplo, $n=100$ ;
• positivo-definido, es decir, con todos los valores propios reales y positivos;
• rango completo;
• todos los elementos diagonales iguales a $1$ ;
• los elementos fuera de la diagonal deben distribuirse razonablemente de manera uniforme en $(-1, 1)$ . La distribución exacta no importa, pero me gustaría tener una cantidad moderadamente grande (por ejemplo, $10\%$ ) de valores moderadamente grandes (por ejemplo, con un valor absoluto de $0.5$ o superior). Básicamente quiero para asegurarse de que $\mathbf C$ es no casi diagonal con todos los elementos fuera de la diagonal $\approx 0$ .

¿Hay una manera simple de hacerlo?

El propósito es usar tales matrices aleatorias para comparar algunos algoritmos que trabajan con matrices de correlación (o covarianza).

Métodos que no funcionan.

Aquí hay algunas formas de generar matrices de correlación aleatorias que conozco, pero que no me funcionan aquí:

1. Genere azar $\mathbf X$de tamaño $s \times n$ , centre, estandarice y forme la matriz de correlación $\mathbf C=\frac{1}{s-1}\mathbf X^\top \mathbf X$. Si$s>n$, esto generalmente dará como resultado que todas las correlaciones fuera de la diagonal estén alrededor de$0$. Si$s\ll n$, algunas correlaciones serán fuertes, pero$\mathbf C$no será rango completo.

2. Genere una matriz definida positiva aleatoria $\mathbf B$ de una de las siguientes maneras:

   • Genere el cuadrado aleatorio y haga que el simétrico positivo sea definido B = A A ⊤ .$\mathbf A$$\mathbf B = \mathbf A \mathbf A^\top$

   • Genere un cuadrado aleatorio , haga simétrico E = A + A ⊤ , y haga que sea positivo definido realizando la descomposición propia E = U S U ⊤ y estableciendo todos los valores propios negativos en cero: B = $\mathbf A$$\mathbf E = \mathbf A + \mathbf A^\top$$\mathbf E = \mathbf U \mathbf S \mathbf U^\top$ . NB: esto dará como resultado una matriz de rango deficiente.$\mathbf B = \mathbf U \:\mathrm{max}\{\mathbf S, \mathbf 0\} \:\mathbf U^\top$

   • Genere ortogonal aleatorio (por ejemplo, generando un cuadrado aleatorio A y haciendo su descomposición QR, o mediante el proceso de Gram-Schmidt) y diagonal D aleatoria con todos los elementos positivos; forma B = Q D Q ⊤ .$\mathbf Q$$\mathbf A$$\mathbf D$$\mathbf B = \mathbf Q \mathbf D \mathbf Q^\top$

   Obtenido matriz puede normalizarse fácilmente tener todos unos en la diagonal: C = D - 1 / 2 B D - 1 / 2 , donde D = delta i un $\mathbf B$$\mathbf C = \mathbf D^{-1/2}\mathbf B \mathbf D^{-1/2}$ es la matriz diagonal con la misma diagonal como B . Las tres formas mencionadas anteriormente para generar B dan como resultado que C tenga elementos fuera de la diagonal cerca de 0 .$\mathbf D = \mathrm{diag}\:\mathbf B$$\mathbf B$$\mathbf B$$\mathbf C$$0$

Actualización: hilos más antiguos

Después de publicar mi pregunta, encontré dos duplicados en el pasado:

• ¿Cómo generar una matriz de correlación aleatoria que tenga entradas fuera de diagonal distribuidas aproximadamente normalmente con una desviación estándar dada?
• ¿Cómo generar de manera eficiente matrices de correlación aleatoria positiva-semidefinida?

Desafortunadamente, ninguno de estos hilos contenía una respuesta satisfactoria (hasta ahora :)

Otras respuestas surgieron con buenos trucos para resolver mi problema de varias maneras. Sin embargo, encontré un enfoque basado en principios que creo que tiene una gran ventaja de ser conceptualmente muy claro y fácil de ajustar.

En este hilo: ¿Cómo generar eficientemente matrices de correlación semidefinidas positivas al azar? - Describí y proporcioné el código para dos algoritmos eficientes de generación de matrices de correlación aleatorias. Ambos provienen de un artículo de Lewandowski, Kurowicka y Joe (2009), al que @ssdecontrol se refirió en los comentarios anteriores (¡muchas gracias!).

Consulte mi respuesta allí para ver muchas figuras, explicaciones y códigos matlab. El llamado método "vine" permite generar matrices de correlación aleatorias con cualquier distribución de correlaciones parciales y puede usarse para generar matrices de correlación con grandes valores fuera de la diagonal. Aquí está la figura de ejemplo de ese hilo:

Lo único que cambia entre subtramas es un parámetro que controla cuánto se concentra la distribución de correlaciones parciales alrededor de .$\pm 1$

Copio mi código para generar estas matrices aquí también, para mostrar que no es más largo que los otros métodos sugeridos aquí. Por favor vea mi respuesta vinculada para algunas explicaciones. Los valores de betaparampara la figura anterior fueron (y la dimensionalidad fue 100 ).${50,20,10,5,2,1}$d$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

Actualización: valores propios

@psarka pregunta sobre los valores propios de estas matrices. En la figura siguiente, trazo los espectros de valores propios de las mismas seis matrices de correlación que las anteriores. Tenga en cuenta que disminuyen gradualmente; en contraste, el método sugerido por @psarka generalmente da como resultado una matriz de correlación con un valor propio grande, pero el resto es bastante uniforme.

Actualizar. Método realmente simple: varios factores

$k<n$$\mathbf W$$k \times n$$\mathbf W \mathbf W^\top$$\mathbf D$$\mathbf B = \mathbf W \mathbf W^\top + \mathbf D$$k={100, 50, 20, 10, 5, 1}$

$k$

Aquí está el código:

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)));

$a$

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()

Hmm, después de haber hecho un ejemplo en mi lenguaje MatMate, veo que ya hay una respuesta de Python, que podría ser preferible porque Python es ampliamente utilizado. Pero debido a que aún tenía preguntas, le muestro mi enfoque usando el lenguaje de matriz de Matmate, tal vez es más autocomentario.

Método 1
(usando 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

El problema aquí podría ser que definimos bloques de submatrices que tienen altas correlaciones dentro con poca correlación entre y esto no es programáticamente sino por las constantes expresiones de concatenación. Tal vez este enfoque podría modelarse de manera más elegante en Python.

{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~

y la matriz de correlación producida es

     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

Posiblemente esto genera una matriz de correlación con componentes principales dominantes debido a la regla de generación acumulativa para la matriz de carga de factores. También podría ser mejor asegurar una definición positiva al hacer que la última parte de la varianza sea un factor único. Lo dejé en el programa para mantener el enfoque en el principio general.

Una matriz de correlación de 100x100 tenía las siguientes frecuencias de correlaciones (redondeadas a 1 lugar dec.)

    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

[actualizar]. Hmm, la matriz 100x100 está mal acondicionada; Pari / GP no puede determinar los valores propios correctamente con la función polroots (charpoly ()), incluso con una precisión de 200 dígitos. Hice una rotación de Jacobi para formar pca en la matriz de carga L y encontré valores propios extremadamente pequeños, los imprimí en logaritmos a la base 10 (que dan aproximadamente la posición del punto decimal). Lea de izquierda a derecha y luego fila por fila:

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

[actualización 2]
Método 2 (b)
Una mejora podría ser aumentar la varianza específica del ítem a un nivel no marginal y reducir a un número razonablemente menor de factores comunes (por ejemplo, la raíz cuadrada entera del número de ítem):

{  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...

La estructura del resultado.

en términos de distribución de correlaciones:

sigue siendo similar (también la desagradable no descomponibilidad de PariGP), pero los valores propios, cuando se encuentran mediante la rotación jacobi de la matriz de carga, ahora tienen una mejor estructura, para un ejemplo recién calculado obtuve los valores propios como

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

Pregunta interesante (como siempre!). ¿Qué tal encontrar un conjunto de matrices de ejemplo que exhiban las propiedades que desea, y luego tomar combinaciones convexas de las mismas, ya que si$A$ y $B$ son positivos definidos, entonces también lo es $\lambda A + (1-\lambda)B$. Como beneficio adicional, no será necesario reescalar las diagonales por la convexidad de la operación. Ajustando el$\lambda$para estar más concentrado hacia 0 y 1 versus distribuido uniformemente, podría concentrar las muestras en los bordes del politopo o en el interior. (Podría usar una distribución beta / Dirichlet para controlar la concentración frente a la uniformidad).

Por ejemplo, podrías dejar $A$ ser simétrica por componentes y $B$Sé toeplitz. Por supuesto, siempre puedes agregar otra clase$C$, y tomar $\lambda_A A + \lambda_B B + \lambda_C C$ tal que $\sum \lambda = 1$ y $\lambda \geq 0$, y así.

R tiene un paquete (clusterGeneration) que implementa el método en:

• Joe, H. (2006) Generando matrices de correlación aleatorias basadas en correlaciones parciales . Revista de Análisis Multivariante, 97, 2177-2189.

Ejemplo:

> (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

Desafortunadamente, no parece posible simular correlaciones que sigan una distribución uniforme con esto. Parece hacer correlaciones más fuertes cuando alphadse establece en valores muy pequeños, pero incluso en 1/100000000000000el rango de correlaciones solo subiría a aproximadamente 1.40.

Sin embargo, espero que esto pueda ser de alguna utilidad para alguien.