Генерация многомерных нормально распределенных случайных чисел в Matlab

Если по определению вы ссылаетесь на плотность многомерного нормального распределения:

он не содержит ни разложения Холецкого, ни матричного квадратного корня из Σ, а его обратного и скалярного квадратного корня из его определителя.

Но для численного генерирования случайных чисел из этого распределения плотность бесполезна. Это даже не самое общее описание многомерного нормального распределения, поскольку формула плотности имеет смысл только для положительно определенных матриц Σ, в то время как распределение также определяется, если есть нулевые собственные значения - это просто означает, что дисперсия равна 0 в направлении соответствующего собственного вектора.

Ваш вопрос следует подходу, чтобы начать со стандартных многомерных нормально распределенных случайных чисел Z, созданных randn, а затем применить линейное преобразование. Предполагая, что mu является p-мерным вектором-строкой, нам нужна nxp-мерная случайная матрица (каждая строка - одно наблюдение, каждый столбец - одна переменная):

Z = randn(n, p);
x = mu + Z * A;

Нам нужна матрица A такая, чтобы ковариация x была Sigma. Поскольку ковариация Z является единичной матрицей, ковариация x задается A' * A. Решение этой проблемы дает разложение Холецкого, поэтому естественным выбором является

A = chol(Sigma);

где A - верхнетреугольная матрица.

Однако мы также можем найти эрмитово решение, A' = A, и тогда A' * A превратится в A^2, квадрат матрицы. Решение этой проблемы дает квадратный корень матрицы, который вычисляется путем замены каждого собственного значения Sigma квадратным корнем (или отрицательным); в общем, существует 2ⁿ возможных решений для n положительных собственных значений. Функция Matlab sqrtm возвращает квадратный корень главной матрицы, который является единственным неотрицательно-определенным решением. Следовательно,

A = sqrtm(Sigma)

тоже работает. A ^ 0.5 в принципе должен поступать так же.

Моделирование с использованием этого кода

p = 10;
n = 1000;

nr = 1000;
cp = nan(nr, 1);
sp = nan(nr, 1);
pp = nan(nr, 1);

for i = 1 : nr
    x = randn(n, p);
    Sigma = cov(x);

    cS = chol(Sigma);
    cp(i) = norm(cS' * cS - Sigma);

    sS = sqrtm(Sigma);
    sp(i) = norm(sS' * sS - Sigma);

    pS = Sigma ^ 0.5;
    pp(i) = norm(pS' * pS - Sigma);
end    

mean([cp sp pp])

yield, что chol более точен, чем два других метода, а профилирование показывает, что он также намного быстрее, как для p = 10, так и для p = 100.

Однако разложение Холецкого имеет тот недостаток, что оно определено только для положительно определенного Σ, в то время как требование матричного квадратного корня состоит просто в том, что Σ является неотрицательно-определенным (sqrtm возвращает предупреждение для единичного ввода, но возвращает действительный результат ).