¿Cómo generar de manera eficiente matrices de correlación aleatoria positiva-semidefinida?

38

Me gustaría poder generar eficientemente matrices de correlación semidefinidas positivas (PSD). Mi método se ralentiza dramáticamente a medida que aumento el tamaño de las matrices que se generarán.

  1. ¿Podría sugerir alguna solución eficiente? Si conoce algún ejemplo en Matlab, estaría muy agradecido.
  2. Al generar una matriz de correlación PSD, ¿cómo elegiría los parámetros para describir las matrices que se generarán? ¿Una correlación promedio, desviación estándar de correlaciones, valores propios?
Eduardas
fuente

Respuestas:

16

Puede hacerlo al revés: cada matriz (el conjunto de todas las matrices PSD simétricas ) se puede descomponer comoCR++pp×p

C=OTDO donde es una matriz ortonormalO

Para obtener , primero genere una base aleatoria (donde son vectores aleatorios, generalmente en ). A partir de ahí, use el proceso de ortogonalización de Gram-Schmidt para obtener( v 1 , . . . , V p ) v i ( - 1 , 1 ) ( u 1 , . . . . , U p ) = OO(v1,...,vp)vi(1,1)(u1,....,up)=O

R tiene varios paquetes que pueden hacer la ortogonalización GS de una base aleatoria de manera eficiente, incluso para grandes dimensiones, por ejemplo, el paquete 'lejano'. Aunque encontrará el algoritmo GS en la wiki, probablemente sea mejor no reinventar la rueda y optar por una implementación de matlab (seguramente existe una, simplemente no puedo recomendar ninguna).

Finalmente, es una matriz diagonal cuyos elementos son todos positivos (esto es, nuevamente, fácil de generar: generar p números aleatorios, cuadrarlos, clasificarlos y colocarlos en la diagonal de una identidad p por p matriz).Dppp

usuario603
fuente
3
(1) Tenga en cuenta que el resultante no será una matriz de correlación (como lo solicitó el OP), porque no tendrá unos en la diagonal. Por supuesto que se puede reescalado tener unos en la diagonal, estableciéndolo en E - 1 / 2 C E - 1 / 2 , donde E es una matriz diagonal con la misma diagonal como C . (2) Si no me equivoco, esto dará como resultado matrices de correlación donde todos los elementos fuera de la diagonal se concentran alrededor de 0 , por lo que no hay flexibilidad que el OP estaba buscando (OP quería poder establecer "una correlación promedio , desviación estándar de correlaciones, valores propios "CE1/2CE1/2EC0)
ameba dice Reinstate Monica
@amoeba: abordaré (2) ya que, como usted señala, la solución a (1) es trivial. La caracterización de un número de la 'forma' (la relación entre los elementos diagonales de entrada y salida) de una matriz PSD (y, por lo tanto, una covarianza y también una matriz de correlación) es su número de condición. Y, el método descrito anteriormente permite un control total sobre él. La 'concentración de los elementos fuera de la diagonal alrededor de 0' no es una característica del método utilizado para generar matrices PSD sino, más bien, una consecuencia de las restricciones necesarias para garantizar que la matriz sea PSD y el hecho de que es grande. p
usuario603
¿Estás diciendo que todas las matrices PSD grandes tienen elementos fuera de la diagonal cercanos a cero? No estoy de acuerdo, no es así. Consulte mi respuesta aquí para ver algunos ejemplos: ¿Cómo generar una matriz de correlación aleatoria que tenga entradas fuera de la diagonal distribuidas aproximadamente normalmente con una desviación estándar dada? Pero uno puede ver directamente que no es el caso, porque una matriz cuadrada que tiene todos en diagonal y un valor fijo cualquier lugar fuera de la diagonal es PSD y ρ puede ser arbitrariamente grande (pero, por supuesto, inferior a 1 ). ρρ1
ameba dice Reinstate Monica
@amoeba: entonces me equivoqué al suponer que, por necesidad, la diagonal fuera de las matrices de correlación grandes cuando se les permite ser tanto positivas como negativas están cerca de 0. Gracias por el ejemplo esclarecedor.
usuario603
1
Leí un artículo muy bueno sobre la generación de matrices de correlación aleatorias y proporcioné mi propia respuesta aquí (así como otra respuesta en ese hilo vinculado). Creo que puede resultarle interesante.
ameba dice Reinstate Monica
27

Un documento que genera matrices de correlación aleatorias basadas en vides y método de cebolla extendida por Lewandowski, Kurowicka y Joe (LKJ), 2009, proporciona un tratamiento unificado y una exposición de los dos métodos eficientes de generar matrices de correlación aleatorias. Ambos métodos permiten generar matrices a partir de una distribución uniforme en un cierto sentido preciso definido a continuación, son fáciles de implementar, rápidos y tienen la ventaja adicional de tener nombres divertidos.

Una matriz simétrica real de tamaño con unos en la diagonal tiene d ( d - 1 ) / 2 elementos únicos fuera de la diagonal y, por lo tanto, se puede parametrizar como un punto en R d ( d - 1 ) / 2 . Cada punto en este espacio corresponde a una matriz simétrica, pero no todos son positivos definidos (como deben ser las matrices de correlación). Las matrices de correlación, por lo tanto, forman un subconjunto de R d ( d - 1 ) / 2d×dd(d1)/2Rd(d1)/2Rd(d1)/2 (en realidad, un subconjunto convexo conectado), y ambos métodos pueden generar puntos a partir de una distribución uniforme sobre este subconjunto.

Proporcionaré mi propia implementación de MATLAB de cada método y los ilustraré con .d=100


Método de cebolla

El método de la cebolla proviene de otro artículo (ref. # 3 en LKJ) y su nombre se debe al hecho de que las matrices de correlación se generan comenzando con una matriz y creciendo columna por columna y fila por fila. La distribución resultante es uniforme. Realmente no entiendo las matemáticas detrás del método (y prefiero el segundo método de todos modos), pero aquí está el resultado:1×1

Método de cebolla

Aquí y debajo, el título de cada subtrama muestra los valores propios más pequeños y más grandes, y el determinante (producto de todos los valores propios). Aquí está el código:

%// ONION METHOD to generate random correlation matrices distributed randomly
function S = onion(d)
    S = 1;
    for k = 2:d
        y = betarnd((k-1)/2, (d-k)/2); %// sampling from beta distribution
        r = sqrt(y);
        theta = randn(k-1,1);
        theta = theta/norm(theta);
        w = r*theta;
        [U,E] = eig(S);
        R = U*E.^(1/2)*U';             %// R is a square root of S
        q = R*w;
        S = [S q; q' 1];               %// increasing the matrix size
    end
end

Método de cebolla extendida

LKJ modifica ligeramente este método para poder muestrear matrices de correlación de una distribución proporcional a [ d e tC . Cuanto mayor sea η , mayor será el determinante, lo que significa que las matrices de correlación generadas se acercarán cada vez más a la matriz de identidad. El valor η = 1 corresponde a una distribución uniforme. En la figura siguiente, las matrices se generan con η = 1 , 10 , 100 , 1000 , 10[detC]η1ηη=1 .η=1,10,100,1000,10000,100000

Método de cebolla extendida

Por alguna razón para obtener el determinante del mismo orden de magnitud que en el método de la cebolla de vainilla, necesito poner y no η = 1 (como afirma LKJ). No estoy seguro de dónde está el error.η=0η=1

%// EXTENDED ONION METHOD to generate random correlation matrices
%// distributed ~ det(S)^eta [or maybe det(S)^(eta-1), not sure]
function S = extendedOnion(d, eta)
    beta = eta + (d-2)/2;
    u = betarnd(beta, beta);
    r12 = 2*u - 1;
    S = [1 r12; r12 1];  

    for k = 3:d
        beta = beta - 1/2;
        y = betarnd((k-1)/2, beta);
        r = sqrt(y);
        theta = randn(k-1,1);
        theta = theta/norm(theta);
        w = r*theta;
        [U,E] = eig(S);
        R = U*E.^(1/2)*U';
        q = R*w;
        S = [S q; q' 1];
    end
end

Método de la vid

El método Vine fue originalmente sugerido por Joe (J en LKJ) y mejorado por LKJ. Me gusta más, porque es conceptualmente más fácil y también más fácil de modificar. La idea es generar correlaciones parciales (son independientes y pueden tener cualquier valor de [ - 1 , 1 ]d(d1)/2[1,1]sin restricciones) y luego conviértalos en correlaciones crudas a través de una fórmula recursiva. Es conveniente organizar el cálculo en un cierto orden, y este gráfico se conoce como "vid". Es importante destacar que si se muestrean correlaciones parciales de distribuciones beta particulares (diferentes para diferentes células en la matriz), la matriz resultante se distribuirá uniformemente. Aquí nuevamente, LKJ introduce un parámetro adicional para muestrear desde una distribución proporcional a [ d e tη . El resultado es idéntico a la cebolla extendida:[detC]η1

Método de la vid

%// VINE METHOD to generate random correlation matrices
%// distributed ~ det(S)^eta [or maybe det(S)^(eta-1), not sure]
function S = vine(d, eta)
    beta = eta + (d-1)/2;   
    P = zeros(d);           %// storing partial correlations
    S = eye(d);

    for k = 1:d-1
        beta = beta - 1/2;
        for i = k+1:d
            P(k,i) = betarnd(beta,beta); %// 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
end

Método de vid con muestreo manual de correlaciones parciales.

±1[0,1][1,1]α=β=50,20,10,5,2,1. Cuanto más pequeños son los parámetros de la distribución beta, más se concentra cerca de los bordes.

Método de viñedo con muestreo manual.

Tenga en cuenta que en este caso no se garantiza que la distribución sea invariante de permutación, por lo que adicionalmente permuto aleatoriamente filas y columnas después de la generación.

%// VINE METHOD to generate random correlation matrices
%// with all partial correlations distributed ~ beta(betaparam,betaparam)
%// rescaled to [-1, 1]
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

Así es como los histogramas de los elementos fuera de la diagonal buscan las matrices anteriores (la variación de la distribución aumenta monotónicamente):

Elementos fuera de diagonal


Actualización: utilizando factores aleatorios

k<dWk×dWWDB=WW+DC=E1/2BE1/2EBk=100,50,20,10,5,1

matrices de correlación aleatoria de factores aleatorios

Y el codigo:

%// FACTOR method
function S = factor(d,k)
    W = randn(d,k);
    S = W*W' + diag(rand(1,d));
    S = diag(1./sqrt(diag(S))) * S * diag(1./sqrt(diag(S)));
end

Aquí está el código de ajuste utilizado para generar las figuras:

d = 100; %// size of the correlation matrix

figure('Position', [100 100 1100 600])
for repetition = 1:6
    S = onion(d);

    %// etas = [1 10 100 1000 1e+4 1e+5];
    %// S = extendedOnion(d, etas(repetition));

    %// S = vine(d, etas(repetition));

    %// betaparams = [50 20 10 5 2 1];
    %// S = vineBeta(d, betaparams(repetition));

    subplot(2,3,repetition)

    %// use this to plot colormaps of S
    imagesc(S, [-1 1])
    axis square
    title(['Eigs: ' num2str(min(eig(S)),2) '...' num2str(max(eig(S)),2) ', det=' num2str(det(S),2)])

    %// use this to plot histograms of the off-diagonal elements
    %// offd = S(logical(ones(size(S))-eye(size(S))));
    %// hist(offd)
    %// xlim([-1 1])
end
ameba dice Reinstate Monica
fuente
2
Este es un resumen fantástico, ¡me alegro de haber dicho algo!
shadowtalker
Cuando traduje el código matlab para la matriz de correlación basada en vid a R y lo probé, la densidad de correlaciones en la columna 1 siempre fue diferente a las columnas posteriores. Puede ser que haya traducido algo incorrectamente, pero tal vez esta nota ayude a alguien.
Charlie
3
Para los usuarios de R, la función rcorrmatrix en el paquete clusterGeneration (escrito por W Qui y H. Joe) implementa el método vine.
RNM
15

AATAyT(ATA)y0yyT(ATA)y=(Ay)TAy=||Ay||que no es negativo Entonces, en Matlab, simplemente intente

A = randn(m,n);   %here n is the desired size of the final matrix, and m > n
X = A' * A;

Dependiendo de la aplicación, esto puede no darle la distribución de los valores propios que desea; La respuesta de Kwak es mucho mejor en ese sentido. Los valores propios Xproducidos por este fragmento de código deben seguir la distribución Marchenko-Pastur.

Para simular las matrices de correlación de acciones, por ejemplo, es posible que desee un enfoque ligeramente diferente:

k = 7;      % # of latent dimensions;
n = 100;    % # of stocks;
A = 0.01 * randn(k,n);  % 'hedgeable risk'
D = diag(0.001 * randn(n,1));   % 'idiosyncratic risk'
X = A'*A + D;
ascii_hist(eig(X));    % this is my own function, you do a hist(eig(X));
-Inf <= x <  -0.001 : **************** (17)
-0.001 <= x <   0.001 : ************************************************** (53)
 0.001 <= x <   0.002 : ******************** (21)
 0.002 <= x <   0.004 : ** (2)
 0.004 <= x <   0.005 :  (0)
 0.005 <= x <   0.007 : * (1)
 0.007 <= x <   0.008 : * (1)
 0.008 <= x <   0.009 : *** (3)
 0.009 <= x <   0.011 : * (1)
 0.011 <= x <     Inf : * (1)
shabbychef
fuente
1
¿estarías dispuesto a compartir tu función ascii_hist por casualidad?
btown
@btown, ¡el margen es demasiado pequeño para contenerlo!
shabbychef
1
yT(ATA)y=(Ay)TAy=||Ay||
8

DA=QDQTQ

JM no es un estadístico
fuente
M.: Buena referencia: esta parece ser la solución más eficiente (asintóticamente).
whuber
3
@whuber: Je, lo recogí de Golub y Van Loan (por supuesto); Lo uso todo el tiempo para ayudar a generar matrices de prueba para rutinas de valor propio / valor singular de prueba de esfuerzo. Como se puede ver en el documento, es esencialmente equivalente a descomponer QR una matriz aleatoria como lo sugerido por kwak, excepto que se hace de manera más eficiente. Hay una implementación de MATLAB de esto en Higham's Text Matrix Toolbox, por cierto.
JM no es un estadístico
M .:> Gracias por la implementación de matlab. ¿Por casualidad conocerías un generador de matriz pseudoaleatoria de Haar en R?
user603
@kwak: No tengo idea, pero si todavía no hay una implementación, no debería ser demasiado difícil traducir el código MATLAB a R (puedo intentar crear uno si realmente no lo hay); El único requisito previo es un generador decente para las variantes normales pseudoaleatorias, que estoy seguro de que R tiene.
JM no es un estadístico
M .:> sí, probablemente lo traduciré yo mismo. Gracias por los enlaces, mejor.
user603
4

No ha especificado una distribución para las matrices. Dos comunes son las distribuciones Wishart y Wishart inversa. La descomposición de Bartlett proporciona una factorización de Cholesky de una matriz aleatoria de Wishart (que también se puede resolver de manera eficiente para obtener una matriz inversa aleatoria de Wishart).

De hecho, el espacio Cholesky es una forma conveniente de generar otros tipos de matrices PSD aleatorias, ya que solo debe asegurarse de que la diagonal no sea negativa.

Simon Byrne
fuente
> No aleatorio: dos matrices generadas a partir del mismo Whishard no serán independientes entre sí. Si planea cambiar el Whishart en cada generación, ¿cómo piensa generar esos Whishart en primer lugar?
user603
@kwak: No entiendo su pregunta: la descomposición de Bartlett dará sorteos independientes de la misma distribución de Wishart.
Simon Byrne
> Permítanme reformular esto, ¿de dónde obtienen la matriz de escala de su distribución de whishart?
user603
1
@kwak: es un parámetro de la distribución, por lo que está arreglado. Lo selecciona al principio, en función de las características deseadas de su distribución (como la media).
Simon Byrne
3

UTSU

alegre
fuente
Si las entradas se generan a partir de una distribución Normal en lugar de una uniforme, la descomposición que mencione debe ser SO (n) invariante (y, por lo tanto, equidistribuida en relación con la medida de Haar).
whuber
Interesante. ¿Puede proporcionar una referencia para esto?
Gappy
1
> el problema con este método es que no puede controlar la relación del valor propio más pequeño a más grande (y creo que a medida que el tamaño de su conjunto de datos generado aleatoriamente es infinito, esta relación convergerá a 1).
user603
1

Si desea tener más control sobre su matriz PSD simétrica generada, por ejemplo, generar un conjunto de datos de validación sintética, tiene una serie de parámetros disponibles. Una matriz PSD simétrica corresponde a una hiper-elipse en el espacio N-dimensional, con todos los grados de libertad relacionados:

  1. Rotaciones
  2. Longitudes de ejes.

Entonces, para una matriz bidimensional (es decir, elipse 2d), tendrá 1 rotación + 2 ejes = 3 parámetros.

Si las rotaciones traen a la mente matrices ortogonales, es un tren correcto, ya que la construcción es nuevamente Σ=OreOT, con Σ siendo la matriz Sym.PSD producida, O la matriz de rotación (que es ortogonal) y rela matriz diagonal, cuyos elementos diagonales controlarán la longitud de los ejes de la elipse.

El siguiente código de Matlab traza 16 conjuntos de datos distribuidos en Gauss bidimensionales basados ​​en Σ, con un ángulo creciente. El código para la generación aleatoria de parámetros está en los comentarios.

figure;
mu = [0,0];
for i=1:16
    subplot(4,4,i)
    theta = (i/16)*2*pi;   % theta = rand*2*pi;
    U=[cos(theta), -sin(theta); sin(theta) cos(theta)];
    % The diagonal's elements control the lengths of the axes
    D = [10, 0; 0, 1]; % D = diag(rand(2,1));    
    sigma = U*D*U';
    data = mvnrnd(mu,sigma,1000);
    plot(data(:,1),data(:,2),'+'); axis([-6 6 -6 6]); hold on;
end

Para más dimensiones, la matriz Diagonal es directa (como arriba), y el U debe derivar de la multiplicación de las matrices de rotación.

PeriRamm
fuente
0

Un enfoque barato y alegre que he usado para probar es generar m N (0,1) n-vectores V [k] y luego usar P = d * I + Suma {V [k] * V [k] '} como una matriz nxn psd. Con m <n esto será singular para d = 0, y para d pequeño tendrá un número de condición alto.


fuente
2
> el problema con este método es que no puede controlar la relación del valor propio más pequeño a más grande (y creo que a medida que el tamaño de su conjunto de datos generado aleatoriamente es infinito, esta relación convergerá a 1).
user603
> además, el método no es muy eficiente (desde un punto de vista computacional)
user603
1
Su "matriz aleatoria" es una estructura especialmente estructurada llamada "matriz diagonal más rango-1" (matriz DR1), por lo que no es realmente una buena matriz aleatoria representativa.
JM no es un estadístico