Я оценил образец ковариационной матрицы образца и получил симметричную матрицу. С , я хотел бы создать -мерного нормальный распределенный гп , но поэтому мне нужно разложение Холецкого . Что мне делать, если не является положительно определенным?C n C C
15
Ответы:
Озабоченность вопроса , как генерировать случайные из случайных величин многомерного нормального распределения с (возможно) единственным числом ковариационной матрицей . Этот ответ объясняет один способ, который будет работать для любой ковариационной матрицы. Это обеспечивает реализацию, которая проверяет его точность.C
R
Алгебраический анализ ковариационной матрицы
Поскольку является ковариационной матрицей, она обязательно является симметричной и положительно-полуопределенной. Чтобы дополнить справочную информацию, пусть µ будет вектором желаемого среднего.C μ
Поскольку симметричен, его разложение по сингулярному значению (SVD) и его собственное разложение автоматически будут иметь видC
для некоторой ортогональной матрицы и диагональной матрицы D 2 . В целом, диагональные элементы D 2 неотрицательны (подразумевая, что все они имеют реальные квадратные корни: выберите положительные, чтобы сформировать диагональную матрицу D ). Информация о C, которую мы имеем, говорит о том, что один или несколько из этих диагональных элементов равны нулю, но это не повлияет ни на одну из последующих операций и не помешает вычислению SVD.V D2 D2 D C
Генерация многомерных случайных значений
Пусть есть стандартное многомерное нормальное распределение: каждый компонент имеет нулевое среднее, единичную дисперсии, и все ковариации равны нуль: ее ковариационная матрица является тождественным я . Тогда случайная величина Y = V D X имеет ковариационную матрицуИкс я Y= V D X
Следовательно, случайная величина имеет многомерное нормальное распределение со средним ц и ковариационной матрицей С .μ + Y μ С
Расчет и пример кода
СледующийY 0 10 , 000 Y 100 С 50
R
код генерирует ковариационную матрицу заданных измерений и ранга, анализирует ее с помощью SVD (или, в закомментированном коде, с собственным разложением), использует этот анализ для генерации заданного числа реализаций (со средним вектором 0 ) и затем сравнивает ковариационную матрицу этих данных с предполагаемой ковариационной матрицей как в числовом, так и в графическом виде. Как показано, он генерирует 10 , 000 реализаций , где размерность Y является 100 и ранг C составляет 50 . ВыходТо есть, ранг данных также и ковариационная матрица по оценкам из данных находится в пределах расстояния 8 × 10 - 5 из C --which близко. В качестве более подробной проверки, коэффициенты C построены по сравнению с его оценкой. Все они лежат близко к линии равенства:50 8 × 10- 5 С С
Код точно соответствует предыдущему анализу и поэтому не требует пояснений (даже для неD2 D
R
пользователей, которые могут эмулировать его в своей любимой прикладной среде). Одна вещь, которую он показывает, - это необходимость соблюдать осторожность при использовании алгоритмов с плавающей точкой: записи могут легко быть отрицательными (но крошечными) из-за неточности. Такие записи должны быть обнулены перед вычислением квадратного корня, чтобы найти сам D.источник
Метод решения A :
В MATLAB код будет
Метод решения B : Сформулируйте и решите выпуклую SDP (полуопределенную программу), чтобы найти ближайшую матрицу D к C в соответствии с нормой их разности Фробениуса, так что D является положительно определенным, указав минимальное собственное значение m.
Используя CVX под MATLAB, код будет:
Сравнение методов решения : Помимо симметризации исходной матрицы, метод решения A корректирует (увеличивает) только диагональные элементы на некоторое общее количество и оставляет недиагональные элементы неизменными. Метод решения B находит ближайшую (к исходной матрице) положительно определенную матрицу, имеющую заданное минимальное собственное значение, в смысле минимальной фробениусовой нормы разности положительно определенной матрицы D и исходной матрицы C, которая основана на суммах квадраты разностей всех элементов D - C, включая недиагональные элементы. Таким образом, регулируя недиагональные элементы, это может уменьшить величину, на которую необходимо увеличить диагональные элементы, и элементы диагонали не обязательно все увеличиваются на одну и ту же величину.
источник
Я бы начал с размышления о модели, которую вы оцениваете.
Если ковариационная матрица не является положительно-полуопределенной, это может означать, что у вас есть проблема коллинеарности в ваших переменных, которая может указывать на проблему с моделью и не обязательно должна решаться численными методами.
Если матрица не является положительной полуопределенной по численным причинам, то есть некоторые решения, о которых можно прочитать здесь
источник
Одним из способов будет вычисление матрицы из разложения по собственным значениям. Теперь я признаю, что не знаю слишком много Math, стоящего за этими процессами, но из моих исследований кажется полезным посмотреть на этот файл справки:
http://stat.ethz.ch/R-manual/R-patched/library/Matrix/html/chol.html
и некоторые другие связанные команды в R.
Также проверьте 'nearPD' в пакете Matrix.
Извините, что не могу помочь, но надеюсь, что мои поиски помогут вам в правильном направлении.
источник
Вы можете получить результаты от функции nearPD в пакете Matrix в R. Это даст вам реальную матрицу обратно.
источник