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

Question

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

• matriz real simétrica cuadrada de $n \times n$ tamaño, 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 no diagonales deben ser razonablemente uniforme distribuido 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 asegurarme de que $\mathbf C$ es no casi diagonal con todos los elementos no diagonales $\approx 0$ .

¿Hay alguna forma sencilla de hacerlo?

El objetivo es utilizar dichas 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. Generar aleatoriamente $\mathbf X$ de $s \times n$ tamaño, centrar, normalizar y formar la matriz de correlación $\mathbf C=\frac{1}{s-1}\mathbf X^\top \mathbf X$ . Si $s>n$ Esto hará que todas las correlaciones no diagonales estén en torno a $0$ . Si $s\ll n$ algunas correlaciones serán fuertes, pero $\mathbf C$ no será de rango completo.

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

   • Generar un cuadrado aleatorio $\mathbf A$ y hacer simétrica definida positiva $\mathbf B = \mathbf A \mathbf A^\top$ .

   • Generar un cuadrado aleatorio $\mathbf A$ , hacer simétrico $\mathbf E = \mathbf A + \mathbf A^\top$ y hacerla positiva definida mediante una descomposición propia $\mathbf E = \mathbf U \mathbf S \mathbf U^\top$ y poniendo todos los valores propios negativos a cero: $\mathbf B = \mathbf U \:\mathrm{max}\{\mathbf S, \mathbf 0\} \:\mathbf U^\top$ . Nota: esto dará como resultado una matriz con rango deficiente.

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

   Matriz obtenida $\mathbf B$ puede ser fácilmente normalizado para tener todos los unos en la diagonal: $\mathbf C = \mathbf D^{-1/2}\mathbf B \mathbf D^{-1/2}$ , donde $\mathbf D = \mathrm{diag}\:\mathbf B$ es la matriz diagonal con la misma diagonal que $\mathbf B$ . Las tres formas mencionadas anteriormente para generar $\mathbf B$ resultan en $\mathbf C$ con elementos no diagonales cercanos a $0$ .

---

Actualización: Hilos antiguos

Después de publicar mi pregunta, he encontrado dos casi duplicados en el pasado:

• ¿Cómo generar una matriz de correlación aleatoria que tenga entradas fuera de la diagonal aproximadamente distribuidas normalmente con una desviación estándar dada?
• ¿Cómo generar eficazmente matrices de correlación aleatorias semidefinidas positivas?

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

Answer 1

En otras respuestas se han propuesto bonitos trucos para resolver mi problema de diversas maneras. Sin embargo, encontré un enfoque basado en principios que creo que tiene la gran ventaja de ser conceptualmente muy claro y fácil de ajustar.

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

Por favor, vea mi respuesta allí para un montón de figuras, explicaciones y código matlab. El método llamado "enredadera" permite generar matrices de correlación aleatorias con cualquier distribución de correlaciones parciales y puede utilizarse para generar matrices de correlación con grandes valores fuera de diagonal. Aquí está la figura de ejemplo de ese hilo:

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

Copio aquí también mi código para generar estas matrices, para mostrar que no es más largo que los otros métodos sugeridos aquí. Por favor, vea mi respuesta enlazada para algunas explicaciones. Los valores de betaparam para la figura anterior fueron ${50,20,10,5,2,1}$ (y la dimensionalidad d fue $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 por los valores propios de estas matrices. En la figura siguiente, he representado los espectros de los valores propios de las mismas seis matrices de correlación. Obsérvese que disminuyen gradualmente; en cambio, el método sugerido por @psarka suele dar como resultado una matriz de correlación con un valor propio grande, pero el resto es bastante uniforme.

---

Actualización. Método realmente sencillo: varios factores

De forma similar a lo que escribió @ttnphns en los comentarios anteriores y @GottfriedHelms en su respuesta, una forma muy sencilla de lograr mi objetivo es generar aleatoriamente varios ( $k<n$ ) cargas factoriales $\mathbf W$ (matriz aleatoria de $k \times n$ tamaño), forman la matriz de covarianza $\mathbf W \mathbf W^\top$ (que por supuesto no será de rango completo) y añadirle una matriz diagonal aleatoria $\mathbf D$ con elementos positivos para hacer $\mathbf B = \mathbf W \mathbf W^\top + \mathbf D$ rango completo. La matriz de covarianza resultante puede normalizarse para convertirse en una matriz de correlación (como se describe en mi pregunta). Esto es muy simple y hace el truco. Aquí están algunas matrices de correlación de ejemplo para $k={100, 50, 20, 10, 5, 1}$ :

El único inconveniente es que la matriz resultante tendrá $k$ grandes valores propios y luego una caída repentina, en contraposición a un bonito decaimiento mostrado anteriormente con el método de la vid. Aquí están los espectros correspondientes:

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

Answer 2

Una cosa simple, pero tal vez funcione para los propósitos de referencia: tomó su 2. e inyectó algunas correlaciones en la matriz de partida. La distribución es Algo así como uniforme, y cambiando $a$ se puede obtener una concentración cercana a 1 y -1 o cercana a 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()

Answer 3

Hmm, después de que yo' hecho un ejemplo en mi MatMate-lenguaje veo que ya hay un python-respuesta, que podría ser preferible porque python es ampliamente utilizado. Pero como todavía tienes preguntas te muestro mi enfoque usando el lenguaje Matmate-matrix, tal vez sea más autocomentario.

Método 1
(Utilizando 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 esto acercarse a podría modelarse de forma más elegante en python.

---

Método 2(a)
Después de eso, hay un enfoque completamente diferente, donde llenamos la posible covarianza restante por cantidades aleatorias del 100% en una matriz de carga de factores. Esto se hace en 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~

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

Es posible que esto genere una matriz de correlación con componentes principales dominantes debido a la regla de generación acumulativa de la matriz de cargas factoriales. También podría ser mejor asegurar la definición positiva haciendo 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 correlaciones 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

[actualización]. Hmm, la matriz 100x100 está mal condicionada; Pari/GP no puede determinar los valores propios correctamente con la función polroots(charpoly())-incluso con 200 dígitos de precisión. He hecho una rotación de Jacobi a pca-forma en la matriz de carga L y encuentro en su mayoría valores propios extremadamente pequeños, imprimiéndolos en logaritmos a base 10 (que dan aproximadamente la posición del punto decimal). Leer 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 algún 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 ítems):

{  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 la distribución de las correlaciones:

sigue siendo similar (también la desagradable no descomponibilidad por PariGP), pero los valores propios, cuando se encuentran por jacobi-rotación de la matriz de cargas, tienen ahora 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

Answer 4

Una pregunta interesante (¡como siempre!). ¿Qué tal si encuentras un conjunto de matrices de ejemplo que exhiban las propiedades que deseas, y luego tomas combinaciones convexas de las mismas, ya que si $A$ y $B$ son positivas definidas, entonces también lo son $\lambda A + (1-\lambda)B$ . Como ventaja, no será necesario reescalar las diagonales, por la convexidad de la operación. Ajustando el $\lambda$ a estar más concentrado hacia 0 y 1 frente a la distribución uniforme, podría concentrar las muestras en los bordes del politopo, o en el interior. (Podría utilizar una distribución beta/Dirichlet para controlar la concentración frente a la uniformidad).

Por ejemplo, puede dejar que $A$ para que sea simétrica en sus componentes, y $B$ sea toeplitz. Por supuesto, siempre se puede añadir 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í sucesivamente.

Answer 5

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

• Joe, H. (2006) Generación de matrices de correlación aleatorias basadas en correlaciones parciales . Journal of Multivariate Analysis, 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

Por desgracia, no parece posible simular correlaciones que sigan una distribución uniforme con esto. Parece hacer correlaciones más fuertes cuando alphad se fija en valores muy pequeños, pero incluso en 1/100000000000000 El rango de correlaciones sólo llegaría a 1,40 aproximadamente.

No obstante, espero que esto pueda ser de utilidad para alguien.