Oszacowałem próbkę macierzy kowariancji próbki i otrzymałem macierz symetryczną. Z C , to proszę utworzyć n -variate normalnego rozproszonego rn a zatem potrzebny jest rozkład Cholesky'iego z C . Co powinienem zrobić, jeśli C nie jest pozytywnie określony?
15
Odpowiedzi:
Pytanie dotyczy sposobu generowania przypadkowych zmiennymi z wielowymiarowym rozkładu normalnego z (ewentualnie) pojedynczej macierzy kowariancji . Ta odpowiedź wyjaśnia jeden ze sposobów, który będzie działał dla dowolnej macierzy kowariancji. Zapewnia implementację, która sprawdza jej dokładność.C
R
Analiza algebraiczna macierzy kowariancji
Ponieważ jest macierzą kowariancji, z konieczności jest symetryczna i dodatnio-pół-skończona. Aby uzupełnić informacje podstawowe, niech μ będzie wektorem pożądanych środków.C μ
Ponieważ jest symetryczny, jego rozkład wartości osobliwej (SVD) i jego składnia elektronowa będą miały automatycznie postaćC
dla niektórych macierzy ortogonalnej i macierzy diagonalnej D 2 . Zasadniczo diagonalne elementy D 2 są nieujemne (co oznacza, że wszystkie mają rzeczywiste pierwiastki kwadratowe: wybierz te dodatnie, aby utworzyć macierz diagonalną D ). Informacje, które posiadamy o C, mówią, że jeden lub więcej z tych elementów ukośnych ma wartość zero - ale nie wpłynie to na żadną z późniejszych operacji ani nie zapobiegnie obliczeniu SVD.V D2 D2 D C
Generowanie losowych wartości na wielu odmianach
Niech mają standardową wielowymiarowego rozkładu normalnego: każdy składnik ma zerową średnią jednostkową, wariancji, kowariancji i wszystkie są zerowe: jego macierz kowariancji jest tożsamość ja . Zatem zmienna losowa Y = V D X ma macierz kowariancjiX I Y=VDX
W związku z tym zmienną losową ma wielowymiarowego rozkładu normalnego o średniej ľ i macierzy kowariancji C .μ+Y μ C
Obliczenia i przykładowy kod
PoniższyY 0 10,000 Y 100 C 50
R
kod generuje macierz kowariancji danych wymiarów i rangi, analizuje ją za pomocą SVD (lub, w skomentowanym kodzie, z kompozycją elektronową), wykorzystuje tę analizę do wygenerowania określonej liczby realizacji (ze średnim wektorem 0 ) , a następnie porównuje macierz kowariancji tych danych z zamierzoną macierzą kowariancji zarówno numerycznie, jak i graficznie. Jak pokazano, generuje 10 , 000 realizacje, w których wymiar Y jest 100 i jest pozycja C to 50 . Dane wyjściowe toOznacza to, że pozycja danych jest również i macierz kowariancji oszacowana na podstawie danych w odległości 8 x 10 - 5 o C --which znajduje się w pobliżu. W ramach bardziej szczegółowej kontroli współczynniki C są wykreślane względem współczynników jego oszacowania. Wszystkie leżą blisko linii równości:50 8×10−5 C C
R
źródło
Metoda rozwiązania A :
W MATLAB byłby to kod
Metoda rozwiązania B : Sformułuj i rozwiąż wypukły SDP (program półfinałowy), aby znaleźć najbliższą macierz D do C zgodnie z normą frobeniusa ich różnicy, tak że D jest dodatnio określony, mając określoną minimalną wartość własną m.
Używając CVX pod MATLAB, kod będzie:
Porównanie metod rozwiązania : Oprócz symetryczności macierzy początkowej, metoda rozwiązania A dostosowuje (zwiększa) tylko elementy ukośne o pewną wspólną ilość i pozostawia elementy nie-ukośne bez zmian. Metoda rozwiązania B znajduje najbliższą (do pierwotnej macierzy) pozytywną określoną macierz o określonej minimalnej wartości własnej, w sensie minimalnej normy Frobeniusa różnicy dodatniej określonej macierzy D i oryginalnej macierzy C, która jest oparta na sumach kwadratowe różnice wszystkich elementów D - C, aby uwzględnić elementy nie przekątne. Tak więc, dostosowując elementy o przekątnej, może zmniejszyć ilość, o którą należy zwiększyć elementy o przekątnych, a elementy o diagoanlu niekoniecznie muszą zostać zwiększone o tę samą ilość.
źródło
Zacznę od przemyślenia modelu, który szacujesz.
Jeśli macierz kowariancji nie jest dodatnią półokreśloną, może to wskazywać, że masz problem ze współliniowością w swoich zmiennych, co wskazywałoby na problem z modelem i niekoniecznie musi być rozwiązane metodami numerycznymi.
Jeśli macierz nie jest dodatnia półfinałowa z powodów numerycznych, istnieje kilka rozwiązań, o których można przeczytać tutaj
źródło
Jednym ze sposobów byłoby obliczenie macierzy na podstawie rozkładu wartości własnych. Przyznaję, że nie znam zbyt wiele matematyki stojącej za tymi procesami, ale z moich badań wydaje się owocne spojrzenie na ten plik pomocy:
http://stat.ethz.ch/R-manual/R-pched/library/Matrix/html/chol.html
i kilka innych powiązanych poleceń w R.
Sprawdź także „nearPD” w pakiecie Matrix.
Przepraszam, że nie mogłem pomóc, ale mam nadzieję, że moje poszukiwania pomogą ci popchnąć cię we właściwym kierunku.
źródło
Możesz uzyskać wyniki z funkcji nearPD w pakiecie Matrix w R. To da ci prawdziwie cenną macierz z powrotem.
źródło