Как эффективно генерировать случайные положительно-полуопределенные корреляционные матрицы?

38

Я хотел бы иметь возможность эффективно генерировать матрицы положительно-полуопределенной (PSD) корреляции. Мой метод значительно замедляется, когда я увеличиваю размер создаваемых матриц.

  1. Не могли бы вы предложить какие-либо эффективные решения? Если вам известны какие-либо примеры в Matlab, я был бы очень благодарен.
  2. Когда вы генерируете матрицу корреляции PSD, как бы вы выбрали параметры для описания матриц, которые будут сгенерированы? Средняя корреляция, стандартное отклонение корреляций, собственные значения?
Eduardas
источник

Ответы:

16

Вы можете сделать это в обратном направлении: каждая матрица (множество всех симметричных PSD матриц) может быть разложена как p × pCR++пп×п

Сзнак равноОTDО где О - ортонормированная матрица

Чтобы получить О , сначала сгенерируйте случайный базис (v1,,,,,vп) (где vя - случайные векторы, обычно в (-1,1) ). Оттуда используйте процесс ортогонализации Грамма-Шмидта, чтобы получить (U1,,,,,,Uп)знак равноО

р имеет несколько пакетов, которые могут эффективно выполнять ортогонализацию GS на случайной основе, даже для больших измерений, например, «дальний» пакет. Несмотря на то, что вы найдете алгоритм GS в вики, вероятно, лучше не изобретать колесо и не приступить к реализации Matlab (она, безусловно, существует, я просто не могу ее рекомендовать).

Наконец, - это диагональные матрицы, все элементы которых положительны (это, опять же, легко сгенерировать: сгенерировать случайных чисел, возвести их в квадрат, отсортировать и поместить их в диагональ тождества по матрице).р р рDппп

user603
источник
3
(1) Обратите внимание, что результирующий не будет корреляционной матрицей (как запрошено OP), потому что он не будет иметь их по диагонали. Конечно , это может быть пересчитаны , чтобы те , по диагонали, установив его на , где представляет собой диагональную матрицу с той же диагонали , как . (2) Если я не ошибаюсь, это приведет к матрицам корреляции, где все недиагональные элементы сконцентрированы вокруг , поэтому нет никакой гибкости, которую искал OP (OP хотел иметь возможность установить «среднюю корреляцию»). , стандартное отклонение корреляций, собственные значения " )Е - 1 / 2 С Е - 1 / 2 Е С 0СЕ-1/2СЕ-1/2ЕС0
говорит амеба Reinstate Monica
@amoeba: я обращусь к (2), поскольку, как вы указываете, решение (1) тривиально. Одной числовой характеристикой «формы» (отношения между входящими и выходящими диагональными элементами) матрицы PSD (и, следовательно, ковариации, а также матрицы корреляции) является номер условия. И метод, описанный выше, позволяет полностью контролировать его. «Концентрация недиагональных элементов около 0» - это не особенность метода, используемого для генерации PSD-матриц, а скорее следствие ограничений, необходимых для обеспечения того, что матрица является PSD и фактом, что велико. п
user603
Вы говорите, что все большие матрицы PSD имеют недиагональные элементы, близкие к нулю? Я не согласен, это не так. Проверьте мой ответ здесь на несколько примеров: Как создать матрицу случайной корреляции, которая имеет приблизительно нормально распределенные недиагональные записи с заданным стандартным отклонением? Но сразу видно, что это не так, потому что квадратная матрица, имеющая все единицы по диагонали и фиксированное значение везде вне диагонали, является PSD, и может быть сколь угодно большим (но, конечно, ниже ). ρ 1ρρ1
говорит амеба: восстанови Монику
@amoeba: тогда я ошибался, предполагая, что по необходимости отрицательная диагональ больших корреляционных матриц, когда им разрешено быть как положительными, так и отрицательными, близка к 0. Спасибо за поучительный пример.
user603
1
Я прочитал очень хорошую статью о генерации матриц случайной корреляции и предоставил здесь свой собственный ответ (а также другой ответ в этой связанной ветке). Я думаю, вам может быть интересно.
говорит амеба, восстанови Монику
27

В статье « Генерация матриц случайной корреляции на основе лоз и расширенного метода лука», выполненной Левандовски, Куровицкой и Джо (LKJ), 2009 г., представлены единый подход и описание двух эффективных методов генерации матриц случайной корреляции. Оба метода позволяют генерировать матрицы из равномерного распределения в определенном точном смысле, определенном ниже, просты в реализации, быстры и имеют дополнительное преимущество, заключающееся в забавных именах.

Вещественная симметричная матрица размера с единицами на диагонали имеет уникальных недиагональных элемента и поэтому может быть параметризована как точка в . Каждая точка в этом пространстве соответствует симметричной матрице, но не все они являются положительно определенными (как это должно быть для корреляционных матриц). Поэтому корреляционные матрицы образуют подмножество (фактически связное выпуклое подмножество), и оба метода могут генерировать точки из равномерного распределения по этому подмножеству.d ( d - 1 ) / 2 R d ( d - 1 ) / 2 R d ( d - 1 ) / 2d×dd(d-1)/2рd(d-1)/2рd(d-1)/2

Я предоставлю свою собственную реализацию MATLAB каждого метода и проиллюстрирую их с .dзнак равно100


Луковый метод

Луковый метод взят из другой статьи (ссылка № 3 в LKJ) и получил свое название от того факта, что матрицы корреляции генерируются, начиная с матрицы и увеличивая ее столбец за столбцом и строку за строкой. Результирующее распределение является равномерным. Я не очень понимаю математику метода (и все равно предпочитаю второй метод), но вот результат:1×1

Луковый метод

Здесь и далее заголовок каждого подзаговора показывает наименьшее и наибольшее собственные значения и определитель (произведение всех собственных значений). Вот код:

%// 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

Расширенный метод лука

LKJ немного изменил этот метод, чтобы иметь возможность выбирать корреляционные матрицы из распределения, пропорционального . Чем больше значение , тем больше будет определитель, означающий, что сгенерированные корреляционные матрицы будут все больше приближаться к единичной матрице. Значение соответствует равномерному распределению. На рисунке ниже матрицы генерируются с . [ д е тС η η = 1 η = 1 , 10 , 100 , 1000 , 10[dеTС]η-1ηηзнак равно1ηзнак равно1,10,100,1000,10000,100000

Расширенный метод лука

По какой-то причине, чтобы получить определитель того же порядка, что и в методе с ванильным луком, мне нужно поставить а не (как утверждает LKJ). Не уверен, где ошибка.η = 1ηзнак равно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

Винный метод

Винный метод был первоначально предложен Джо (J в LKJ) и улучшен LKJ. Мне это нравится больше, потому что это концептуально проще, а также легче модифицировать. Идея состоит в том, чтобы сгенерировать частичных корреляции (они независимы и могут иметь любые значения из без каких-либо ограничений), а затем преобразовать их в необработанные корреляции с помощью рекурсивной формулы. Удобно организовать вычисления в определенном порядке, и этот график известен как «лоза». Важно отметить, что если частичные корреляции выбираются из конкретных бета-распределений (различающихся для разных ячеек в матрице), то полученная матрица будет распределена равномерно. Здесь снова LKJ вводит дополнительный параметр для выборки из распределения, пропорционального[ - 1 , 1 ] η [ d e td(d-1)/2[-1,1]η[dеTС]η-1 . Результат идентичен расширенному луку:

Винный метод

%// 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

Винный метод с ручной выборкой частичных корреляций

Как видно выше, равномерное распределение приводит к почти диагональным корреляционным матрицам. Но можно легко изменить метод виноградной лозы, чтобы иметь более сильные корреляции (это не описано в статье LKJ, но просто): для этого следует выбрать частичные корреляции из распределения, сосредоточенного около . Ниже я приведу выборку их из бета-распределения (масштабируется с до ) с помощью . Чем меньше параметры бета-распределения, тем больше оно сосредоточено вблизи краев.[ 0 , 1 ] [ - 1 , 1 ] α = β = 50 , 20 , 10 , 5 , 2 , 1±1[0,1][-1,1]αзнак равноβзнак равно50,20,10,5,2,1

Винный метод с ручным отбором проб

Обратите внимание, что в этом случае распределение не обязательно будет инвариантным по перестановке, поэтому я дополнительно случайным образом переставляем строки и столбцы после генерации.

%// 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

Вот как выглядят гистограммы недиагональных элементов для вышеприведенных матриц (дисперсия распределения монотонно увеличивается):

Недиагональные элементы


Обновление: с использованием случайных факторов

К<dWК×dWWDС = Е - 1 / 2 Б Е - 1 / 2 Е В к = 100 , 50 , 20 , 10 , 5 , 1Взнак равноWW+DСзнак равноЕ-1/2ВЕ-1/2 , где является диагональной матрицей с той же диагонали , как . Это очень просто и делает свое дело. Вот несколько примеров корреляционных матриц для :ЕВКзнак равно100,50,20,10,5,1

матрицы случайной корреляции от случайных факторов

И код:

%// 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

Вот код переноса, использованный для генерации фигур:

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
амеба говорит восстановить монику
источник
2
Это фантастическое краткое изложение, я рад, что сказал что-то!
Shadowtalker
Когда я перевел код matlab для корреляционной матрицы на основе виноградной лозы в R и проверил его, плотность корреляций в столбце 1 всегда отличалась от более поздних столбцов. Может быть, я что-то перевел неправильно, но, возможно, эта заметка кому-то поможет.
Чарли
3
Для пользователей R функция rcorrmatrix в пакете clusterGeneration (написанная W Qui и H. Joe) реализует метод vine.
RNM
15

AATAyT(ATA)y0yYT(ATA)Yзнак равно(AY)TAYзнак равно||AY||который неотрицателен. Так что в Matlab просто попробуйте

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

В зависимости от приложения, это может не дать вам распределение собственных значений, которые вы хотите; Ответ Квака намного лучше в этом отношении. Собственные значения, Xполученные этим фрагментом кода, должны соответствовать распределению Марченко-Пастура.

Скажем, для моделирования корреляционных матриц акций может потребоваться несколько иной подход:

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
источник
1
Вы бы хотели поделиться своей функцией ascii_hist случайно?
Btown
@btown поле слишком мало, чтобы его содержать!
Шаббычеф
1
YT(ATA)Yзнак равно(AY)TAYзнак равно||AY||
8

DA=QDQTQ

JM не является статистиком
источник
М.: Хорошая справка: это наиболее эффективное решение (асимптотически).
whuber
3
@whuber: Хех, я взял это от Голуба и Ван Лоана (конечно); Я использую это все время, чтобы помочь в создании тестовых матриц для стресс-тестирования подпрограмм собственных значений / сингулярных значений. Как видно из статьи, это по сути эквивалентно QR-разложению случайной матрицы, подобной предложенной Kwak, за исключением того, что это делается более эффективно. BTL реализует это в MATLAB.
JM не является статистиком
M.:> Спасибо за реализацию Matlab. Вы случайно не узнали бы о генераторе псевдослучайных матриц Хаара в R?
user603
@kwak: Понятия не имею, но если реализации пока нет, не должно быть слишком сложно перевести код MATLAB в R (я могу попытаться поднять его, если его действительно нет); единственной предпосылкой является приличный генератор псевдослучайных нормальных переменных, который, я уверен, имеет R.
JM не является статистиком
М .: Да, я, наверное, переведу это сам. Спасибо за ссылки, Бест.
user603
4

Вы не указали распределение для матриц. Два распространенных варианта - это рассылки Wishart и обратные рассылки Wishart. Разложение Бартлетта дает Холецкого факторизации матрицы случайных Уишарта (которые также могут быть эффективно решены , чтобы получить случайную обратную Уишарт матрицу).

На самом деле, пространство Холецкого - это удобный способ генерировать другие типы случайных PSD-матриц, так как нужно только убедиться, что диагональ неотрицательна.

Саймон Бирн
источник
> Не случайно: две матрицы, сгенерированные из одного и того же Whishard, не будут независимы друг от друга. Если вы планируете менять Whishart в каждом поколении, то как вы собираетесь создавать эти Whishart в первую очередь?
user603
@kwak: Я не понимаю вашего вопроса: декомпозиция Бартлетта даст независимые ничьи из того же распределения Вишарта.
Саймон Бирн
> Позвольте мне перефразировать это, откуда вы берете масштабную матрицу вашего распределения Whishart?
user603
1
@kwak: это параметр распределения, поэтому он исправлен. Вы выбираете его в начале, основываясь на желаемых характеристиках вашего дистрибутива (таких как среднее значение).
Саймон Бирн
3

UTSU

с промежутками
источник
Если записи генерируются из нормального распределения, а не из униформы, то упомянутое вами разложение должно быть SO (n) -инвариантным (и, следовательно, равнораспределенным относительно меры Хаара).
whuber
Интересный. Можете ли вы предоставить ссылку для этого?
gappy
1
> проблема этого метода заключается в том, что вы не можете контролировать отношение наименьшего к наибольшему собственному значению (и я думаю, что, поскольку размер вашего случайно сгенерированного набора данных уходит в бесконечность, это отношение будет сходиться к 1).
user603
1

Если вы хотите лучше контролировать сгенерированную симметричную матрицу PSD, например, создать синтетический набор данных проверки, у вас есть несколько доступных параметров. Симметричная PSD-матрица соответствует гиперэллипсу в N-мерном пространстве со всеми соответствующими степенями свободы:

  1. Повороты.
  2. Длина осей.

Таким образом, для 2-мерной матрицы (т.е. 2d эллипса) у вас будет 1 поворот + 2 оси = 3 параметра.

Σзнак равноОDОTΣОD

Σ

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

U

PeriRamm
источник
0

Дешевый и веселый подход, который я использовал для тестирования, состоит в том, чтобы сгенерировать m N (0,1) n-векторов V [k], а затем использовать P = d * I + Sum {V [k] * V [k] '} как матрица nxn psd. При m <n это будет единственное число для d = 0, а для малых d будет иметь большое число условий.


источник
2
> проблема этого метода заключается в том, что вы не можете контролировать отношение наименьшего к наибольшему собственному значению (и я думаю, что, поскольку размер вашего случайно сгенерированного набора данных уходит в бесконечность, это отношение будет сходиться к 1).
user603
> кроме того, метод не очень эффективен (с вычислительной точки зрения)
user603
1
Ваша «случайная матрица» - это специально структурированная матрица, называемая «матрица диагональ плюс ранг-1» (матрица DR1), поэтому она не является хорошей репрезентативной случайной матрицей.
JM не является статистиком