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

15

Я оценил образец ковариационной матрицы образца и получил симметричную матрицу. С , я хотел бы создать -мерного нормальный распределенный гп , но поэтому мне нужно разложение Холецкого . Что мне делать, если не является положительно определенным?C n C CССNСС

Клаус
источник
1
В чем разница с этим вопросом stackoverflow.com/questions/17295627/… ?
Дикоа
1
Положительно-полуопределенные матрицы имеют несколько квадратных корней (см. , Например, объяснение в конце stats.stackexchange.com/a/71303/919 ). Вам не обязательно нужен тот, который получен в результате разложения Холецкого. В этом суть проблемы: найти метод вычисления квадратных корней, который работает, даже если матрица является единственной. @amoeba Название предполагает, что ваша интерпретация верна.
whuber

Ответы:

8

Озабоченность вопроса , как генерировать случайные из случайных величин многомерного нормального распределения с (возможно) единственным числом ковариационной матрицей . Этот ответ объясняет один способ, который будет работать для любой ковариационной матрицы. Это обеспечивает реализацию, которая проверяет его точность.CR


Алгебраический анализ ковариационной матрицы

Поскольку является ковариационной матрицей, она обязательно является симметричной и положительно-полуопределенной. Чтобы дополнить справочную информацию, пусть µ будет вектором желаемого среднего.Cμ

Поскольку симметричен, его разложение по сингулярному значению (SVD) и его собственное разложение автоматически будут иметь видC

C=VD2V

для некоторой ортогональной матрицы и диагональной матрицы D 2 . В целом, диагональные элементы D 2 неотрицательны (подразумевая, что все они имеют реальные квадратные корни: выберите положительные, чтобы сформировать диагональную матрицу D ). Информация о C, которую мы имеем, говорит о том, что один или несколько из этих диагональных элементов равны нулю, но это не повлияет ни на одну из последующих операций и не помешает вычислению SVD.VD2D2DС

Генерация многомерных случайных значений

Пусть есть стандартное многомерное нормальное распределение: каждый компонент имеет нулевое среднее, единичную дисперсии, и все ковариации равны нуль: ее ковариационная матрица является тождественным я . Тогда случайная величина Y = V D X имеет ковариационную матрицуИксяYзнак равноВDИкс

Cov(Y)знак равноЕ(YY')знак равноЕ(ВDИксИкс'D'В')знак равноВDЕ(ИксИкс')DВ'знак равноВDяDВ'знак равноВD2В'знак равноС,

Следовательно, случайная величина имеет многомерное нормальное распределение со средним ц и ковариационной матрицей С .μ+YμС

Расчет и пример кода

Следующий Rкод генерирует ковариационную матрицу заданных измерений и ранга, анализирует ее с помощью SVD (или, в закомментированном коде, с собственным разложением), использует этот анализ для генерации заданного числа реализаций (со средним вектором 0 ) и затем сравнивает ковариационную матрицу этих данных с предполагаемой ковариационной матрицей как в числовом, так и в графическом виде. Как показано, он генерирует 10 , 000 реализаций , где размерность Y является 100 и ранг C составляет 50 . ВыходY010,000Y100С50

        rank           L2 
5.000000e+01 8.846689e-05 

То есть, ранг данных также и ковариационная матрица по оценкам из данных находится в пределах расстояния 8 × 10 - 5 из C --which близко. В качестве более подробной проверки, коэффициенты C построены по сравнению с его оценкой. Все они лежат близко к линии равенства:508×10-5СС

фигура

Код точно соответствует предыдущему анализу и поэтому не требует пояснений (даже для не Rпользователей, которые могут эмулировать его в своей любимой прикладной среде). Одна вещь, которую он показывает, - это необходимость соблюдать осторожность при использовании алгоритмов с плавающей точкой: записи могут легко быть отрицательными (но крошечными) из-за неточности. Такие записи должны быть обнулены перед вычислением квадратного корня, чтобы найти сам D.D2D

n <- 100         # Dimension
rank <- 50
n.values <- 1e4  # Number of random vectors to generate
set.seed(17)
#
# Create an indefinite covariance matrix.
#
r <- min(rank, n)+1
X <- matrix(rnorm(r*n), r)
C <- cov(X)
#
# Analyze C preparatory to generating random values.
# `zapsmall` removes zeros that, due to floating point imprecision, might
# have been rendered as tiny negative values.
#
s <- svd(C)
V <- s$v
D <- sqrt(zapsmall(diag(s$d)))
# s <- eigen(C)
# V <- s$vectors
# D <- sqrt(zapsmall(diag(s$values)))
#
# Generate random values.
#
X <- (V %*% D) %*% matrix(rnorm(n*n.values), n)
#
# Verify their covariance has the desired rank and is close to `C`.
#
s <- svd(Sigma <- cov(t(X)))
(c(rank=sum(zapsmall(s$d) > 0), L2=sqrt(mean(Sigma - C)^2)))

plot(as.vector(C), as.vector(Sigma), col="#00000040",
     xlab="Intended Covariances",
     ylab="Estimated Covariances")
abline(c(0,1), col="Gray")
Whuber
источник
2
+1 но когда вы говорите «неопределенный» в своем первом предложении, что именно вы имеете в виду? Я проверил в Википедии, и он говорит, что положительный полуопределенный не является неопределенным, то есть неопределенный означает, что C имеет как положительные, так и отрицательные собственные значения. Это то, что вы имеете в виду там?
говорит амеба: восстанови Монику
2
@amoeba Да, это была ошибка. Спасибо, что заметили. «Неопределенный» означает, что подпись матрицы имеет как положительные, так и отрицательные знаки, тогда как «полуопределенный» означает, что подпись имеет только один знак.
whuber
6

Метод решения A :

  1. 0,5(С+СT)
  2. D+(м-мяN(еяграммеNvaLUе(D)))я

В MATLAB код будет

D = 0.5 * (C + C');
D =  D + (m - min(eig(CD)) * eye(size(D));

Метод решения B : Сформулируйте и решите выпуклую SDP (полуопределенную программу), чтобы найти ближайшую матрицу D к C в соответствии с нормой их разности Фробениуса, так что D является положительно определенным, указав минимальное собственное значение m.

Используя CVX под MATLAB, код будет:

n = size(C,1);
cvx_begin
variable D(n,n)
minimize(norm(D-C,'fro'))
D -m *eye(n) == semidefinite(n)
cvx_end

Сравнение методов решения : Помимо симметризации исходной матрицы, метод решения A корректирует (увеличивает) только диагональные элементы на некоторое общее количество и оставляет недиагональные элементы неизменными. Метод решения B находит ближайшую (к исходной матрице) положительно определенную матрицу, имеющую заданное минимальное собственное значение, в смысле минимальной фробениусовой нормы разности положительно определенной матрицы D и исходной матрицы C, которая основана на суммах квадраты разностей всех элементов D - C, включая недиагональные элементы. Таким образом, регулируя недиагональные элементы, это может уменьшить величину, на которую необходимо увеличить диагональные элементы, и элементы диагонали не обязательно все увеличиваются на одну и ту же величину.

Марк Л. Стоун
источник
2

Я бы начал с размышления о модели, которую вы оцениваете.

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

Если матрица не является положительной полуопределенной по численным причинам, то есть некоторые решения, о которых можно прочитать здесь

johneric
источник
1
Предполагается, что модель является линейной смешанной моделью. И для этого случая не имеет смысла находить правильную модель для данных, скорее данные приведены в качестве примера для некоторых расчетов. Теперь есть вероятность, что вы получите неположительную полуопределенную матрицу в качестве оценки ковариации. Так что же делать, если я хочу выяснить ковариацию по нормальному распределенному населению, откуда поступают данные. То, что выборка распределена нормально, является предположением.
Клаус
1

Одним из способов будет вычисление матрицы из разложения по собственным значениям. Теперь я признаю, что не знаю слишком много Math, стоящего за этими процессами, но из моих исследований кажется полезным посмотреть на этот файл справки:

http://stat.ethz.ch/R-manual/R-patched/library/Matrix/html/chol.html

и некоторые другие связанные команды в R.

Также проверьте 'nearPD' в пакете Matrix.

Извините, что не могу помочь, но надеюсь, что мои поиски помогут вам в правильном направлении.

Фрэнк П.
источник
Привет, спасибо за ссылки. В зависимости от разложения по собственным значениям, это разложение не помогает, потому что оттуда вы получаете сложные собственные значения для квадратно-корневой матрицы, но мне нужна матрица с переопределенными значениями.
Клаус
1

Вы можете получить результаты от функции nearPD в пакете Matrix в R. Это даст вам реальную матрицу обратно.

library(Matrix)
A <- matrix(1, 3,3); A[1,3] <- A[3,1] <- 0
n.A <- nearPD(A, corr=T, do2eigen=FALSE)
n.A$mat

# 3 x 3 Matrix of class "dpoMatrix"
#           [,1]      [,2]      [,3]
# [1,] 1.0000000 0.7606899 0.1572981
# [2,] 0.7606899 1.0000000 0.7606899
# [3,] 0.1572981 0.7606899 1.0000000
Доктор майк
источник
Для пользователей R .. в моем ответе это может быть не плохая «бедная» версия (с меньшим контролем) метода решения B.
Марк Л. Стоун
Я согласен, что это не оптимально, но иногда это помогает.
Доктор Майк