Generuj normalnie rozmieszczone liczby losowe z nieokreśloną dodatnio macierzą kowariancji

15

Oszacowałem próbkę macierzy kowariancji próbki i otrzymałem macierz symetryczną. Z , to proszę utworzyć -variate normalnego rozproszonego rn a zatem potrzebny jest rozkład Cholesky'iego z . Co powinienem zrobić, jeśli nie jest pozytywnie określony? $C$ $C$ $n$ $C$ $C$

— Klaus
źródło

1

Jaka jest różnica z tym pytaniem stackoverflow.com/questions/17295627/... ?

— dickoa

1

Macierze dodatnio-pół-skończone mają wiele pierwiastków kwadratowych (na przykład wyjaśnienie na końcu stats.stackexchange.com/a/71303/919 ). Niekoniecznie potrzebujesz tego, który powstał w wyniku rozkładu Cholesky'ego. Na tym polega sedno problemu: znaleźć metodę obliczania pierwiastków kwadratowych, która działa nawet wtedy, gdy macierz jest pojedyncza. @amoeba Tytuł sugeruje, że twoja interpretacja jest poprawna.

— whuber

8

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ść. $\mathbb{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. $\mathbb{C}$ $\mu$

Ponieważ jest symetryczny, jego rozkład wartości osobliwej (SVD) i jego składnia elektronowa będą miały automatycznie postać $\mathbb{C}$

C = V D^{2} V^{'}

$\mathbb{C} = \mathbb{V\, D^2\, V^\prime}$

dla niektórych macierzy ortogonalnej i macierzy diagonalnej . Zasadniczo diagonalne elementy są nieujemne (co oznacza, że wszystkie mają rzeczywiste pierwiastki kwadratowe: wybierz te dodatnie, aby utworzyć macierz diagonalną ). Informacje, które posiadamy o 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. $\mathbb{V}$ $\mathbb{D}^2$ $\mathbb{D}^2$ $\mathbb{D}$ $\mathbb{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ść . Zatem zmienna losowa ma macierz kowariancji $X$ $\mathbb{I}$ $Y=\mathbb{VD}X$

Cov (Y) = E (Y Y^{'}) = E (V D X X^{'} D^{'} V^{'}) = V D E (X X^{'}) D V^{'} = V D I D V^{'} = V D^{2} V^{'} = C .

$\operatorname{Cov}(Y) = \mathbb{E}(Y Y^\prime) = \mathbb{E}(\mathbb{V D}X\, X^\prime \mathbb{D^\prime V^\prime}) = \mathbb{V D}\mathbb{E}(X X^\prime)\mathbb{D V^\prime} = \mathbb{V D I D V^\prime} = \mathbb{V D^2 V^\prime} = \mathbb{C}.$

W związku z tym zmienną losową ma wielowymiarowego rozkładu normalnego o średniej i macierzy kowariancji . $\mu + \mathbb{Y}$ $\mu$ $\mathbb{C}$

Obliczenia i przykładowy kod

Poniższy Rkod 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 ) , a następnie porównuje macierz kowariancji tych danych z zamierzoną macierzą kowariancji zarówno numerycznie, jak i graficznie. Jak pokazano, generuje realizacje, w których wymiar jest i jest pozycja to . Dane wyjściowe to $Y$ $0$ $10,000$ $Y$ $100$ $C$ $50$

        rank           L2 
5.000000e+01 8.846689e-05

Oznacza to, że pozycja danych jest również i macierz kowariancji oszacowana na podstawie danych w odległości o --which znajduje się w pobliżu. W ramach bardziej szczegółowej kontroli współczynniki są wykreślane względem współczynników jego oszacowania. Wszystkie leżą blisko linii równości: $50$ $8\times 10^{-5}$ $C$ $C$

R $\mathbb{D}^2$ $\mathbb{D}$

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
źródło

2

+1, ale kiedy mówisz „nieokreślony” w pierwszym zdaniu, co dokładnie masz na myśli? Sprawdziłem w Wikipedii i jest napisane, że dodatnia półfinał nie jest nieokreślony, tzn. Nieokreślony oznacza, że C ma zarówno dodatnie, jak i ujemne wartości własne. Czy to masz na myśli?

— ameba mówi Przywróć Monikę

2

@amoeba Tak, to był poślizg. Dzięki za zauważenie. „Nieokreślony” oznacza, że podpis macierzy ma zarówno znaki dodatnie, jak i ujemne, podczas gdy „półfinał” oznacza, że podpis ma tylko jeden znak.

— whuber

6

Metoda rozwiązania A :

$0.5(C + C^T)$
Dodaj wielokrotność macierzy Tożsamości do symetrycznego C wystarczającą, aby uzyskać wynik dodatni z dowolnym pożądanym marginesem, m, tj. Tak, aby najmniejsza wartość własna nowej macierzy miała minimalną wartość własną = m. W szczególności D <- $D + (m - min(eigenvalue(D)))I$

W MATLAB byłby to kod

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

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:

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

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ść.

— Mark L. Stone
źródło

2

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

— johneric
źródło

1

Zakłada się, że model jest liniowym modelem mieszanym. W tym przypadku znalezienie właściwego modelu danych nie jest istotne, raczej dane podano jako przykład niektórych obliczeń. Teraz istnieje możliwość otrzymania nie dodatniej macierzy półfinału jako oszacowania dla kowariancji. Co więc zrobić, jeśli chcę dowiedzieć się o kowariancji z normalnej populacji rozproszonej, z której pochodzą dane. To, że próbka jest normalnie rozłożona, jest przypuszczeniem.

— Klaus,

1

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.

— Frank P.
źródło

Cześć, dzięki za linki. W odniesieniu do rozkładu wartości własnych rozkład ten nie pomaga, ponieważ stamtąd otrzymujesz złożone wartości własne dla macierzy pierwiastkowej, ale potrzebuję macierzy wartościowanej ponownie.

— Klaus,

1

Możesz uzyskać wyniki z funkcji nearPD w pakiecie Matrix w R. To da ci prawdziwie cenną macierz z powrotem.

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

— Dr Mike
źródło

Dla użytkowników R .. może to nie być zła wersja „biednego człowieka” (z mniejszą kontrolą) Rozwiązania Metoda B w mojej odpowiedzi.

— Mark L. Stone,

Zgadzam się, że to nie jest optymalne, ale czasami to działa.

— Dr Mike