Generieren Sie normalverteilte Zufallszahlen mit einer nicht positiv definierten Kovarianzmatrix

15

Ich habe die Sample-Kovarianzmatrix eines Samples geschätzt und eine symmetrische Matrix erhalten. Mit , würde Ich mag schaffen -variate normalverteilt rn aber deshalb muss ich die Cholesky - Zerlegung von . Was soll ich tun, wenn nicht eindeutig ist? $C$ $C$ $n$ $C$ $C$

r random-generation covariance-matrix multivariate-normal cholesky Klaus
quelle

1

Was ist der Unterschied zu dieser Frage stackoverflow.com/questions/17295627/… ?

Dickoa

1

Positive semidefinite Matrizen haben mehrere Quadratwurzeln (siehe die Erklärung am Ende von stats.stackexchange.com/a/71303/919 zum Beispiel). Sie brauchen nicht unbedingt die, die durch die Cholesky-Zerlegung entsteht. Darin liegt der Kern des Problems: Finden Sie eine Methode zur Berechnung von Quadratwurzeln, die auch dann funktioniert, wenn die Matrix singulär ist. @amoeba Der Titel deutet darauf hin, dass Ihre Interpretation korrekt ist.

Whuber

8

Die Frage ist, wie aus einer multivariaten Normalverteilung mit einer (möglicherweise) singulären Kovarianzmatrix Zufallsvariablen erzeugt werden können . Diese Antwort erklärt einen Weg, der für jede Kovarianzmatrix funktioniert . Es bietet eine Implementierung, die die Genauigkeit testet. $\mathbb{C}$ R

Algebraische Analyse der Kovarianzmatrix

Da eine Kovarianzmatrix ist, ist sie notwendigerweise symmetrisch und positiv-semidefinit. Um die Hintergrundinformationen zu vervollständigen, sei der Vektor des gewünschten Mittels. $\mathbb{C}$ $\mu$

Da symmetrisch ist, haben seine Singular Value Decomposition (SVD) und seine eigendecomposition automatisch die Form $\mathbb{C}$

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

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

für einige orthogonale Matrix und diagonale Matrix . Im Allgemeinen sind die diagonalen Elemente von nichtnegativ (was bedeutet, dass sie alle echte Quadratwurzeln haben: Wählen Sie die positiven aus, um die diagonale Matrix zu bilden ). Die Informationen, die wir über besagen, dass eines oder mehrere dieser diagonalen Elemente Null sind - dies hat jedoch keine Auswirkungen auf die nachfolgenden Operationen und verhindert auch nicht, dass die SVD berechnet wird. $\mathbb{V}$ $\mathbb{D}^2$ $\mathbb{D}^2$ $\mathbb{D}$ $\mathbb{C}$

Generierung multivariater Zufallswerte

Es sei eine Standard multivariate Normalverteilung: Jede Komponente hat einen Mittelwert von Null, Einheitsvarianz, und alle Kovarianzen gleich Null sind: seine Kovarianzmatrix die Identität . Dann wird der Zufallsvariable hat Kovarianzmatrix $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 ich 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}.$

Folglich ist die Zufallsvariable hat eine multivariate Normalverteilung mit Mittelwert und die Kovarianzmatrix . $\mu + \mathbb{Y}$ $\mu$ $\mathbb{C}$

Berechnungs- und Beispielcode

Der folgende RCode generiert eine Kovarianzmatrix mit den angegebenen Dimensionen und Rängen, analysiert sie mit der SVD (oder verwendet diese Analyse im auskommentierten Code mit einer Neukomposition), um eine bestimmte Anzahl von Realisierungen von (mit dem Mittelwertvektor ) zu generieren. und vergleicht dann die Kovarianzmatrix dieser Daten sowohl numerisch als auch grafisch mit der beabsichtigten Kovarianzmatrix. Wie gezeigt, erzeugt sie Realisierungen wo die Dimension ist , und der Rang ist . Die Ausgabe ist $Y$ $0$ $10,000$ $Y$ $100$ $C$ $50$

        rank           L2 
5.000000e+01 8.846689e-05

Das heißt, dass der Rang der Daten auch und die Kovarianzmatrix wie aus den Daten geschätzt ist innerhalb Abstand von --which nahe ist. Zur genaueren Überprüfung werden die Koeffizienten von gegen die seiner Schätzung aufgetragen. Sie alle liegen nahe an der Linie der Gleichheit: $50$ $8\times 10^{-5}$ $C$ $C$

Der Code entspricht genau der vorhergehenden Analyse und sollte daher selbsterklärend sein (auch für Nichtbenutzer R, die ihn möglicherweise in ihrer bevorzugten Anwendungsumgebung emulieren). Es zeigt sich, dass bei der Verwendung von Gleitkomma-Algorithmen Vorsicht geboten ist: Die Eingaben von können aufgrund von Ungenauigkeiten leicht negativ (aber winzig) sein. Solche Einträge müssen auf Null gesetzt werden, bevor die Quadratwurzel berechnet wird, um selbst zu finden . $\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
quelle

2

+1 aber wenn du in deinem ersten Satz "unbestimmt" sagst, was genau meinst du damit? Ich habe in Wikipedia nachgesehen und es heißt, dass positives Semidefinit nicht unbestimmt ist, dh unbestimmt bedeutet, dass C sowohl positive als auch negative Eigenwerte hat. Ist es das, was du dort meinst?

Amöbe sagt Reinstate Monica

2

@amoeba Ja, das war ein Ausrutscher. Danke fürs bemerken. "Indefinit" bedeutet, dass die Signatur der Matrix sowohl positive als auch negative Vorzeichen hat, während "semidefinit" bedeutet, dass die Signatur nur ein Vorzeichen hat.

Whuber

6

Lösungsmethode A :

$0.5(C + C^T)$
$D + (m - min(eigenvalue(D)))I$

In MATLAB wäre der Code

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

Lösungsmethode B : Formulieren und lösen Sie ein konvexes SDP (Semidefinite Program), um die nächstgelegene Matrix D zu C gemäß der Frobenius-Norm ihrer Differenz zu finden, so dass D positiv definit ist und den minimalen Eigenwert m angibt.

Mit CVX unter MATLAB wäre der Code:

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

Vergleich der Lösungsmethoden: Abgesehen von der Symmetrisierung der Ausgangsmatrix werden bei der Lösungsmethode A nur die diagonalen Elemente um einen gemeinsamen Betrag angepasst (vergrößert) und die nicht diagonalen Elemente unverändert gelassen. Die Lösungsmethode B findet die nächstgelegene (der ursprünglichen Matrix) positive definite Matrix mit dem angegebenen minimalen Eigenwert im Sinne der minimalen Frobenius-Norm aus der Differenz der positiven definierten Matrix D und der ursprünglichen Matrix C, die auf den Summen von basiert Quadratische Differenzen aller Elemente von D - C, einschließlich der nicht diagonalen Elemente. Durch Anpassen von nicht diagonalen Elementen kann der Betrag verringert werden, um den diagonale Elemente vergrößert werden müssen, und diagonale Elemente werden nicht unbedingt alle um den gleichen Betrag vergrößert.

Mark L. Stone
quelle

2

Zunächst würde ich über das Modell nachdenken, das Sie schätzen.

Wenn eine Kovarianzmatrix nicht positiv und semidefinit ist, weist dies möglicherweise darauf hin, dass Ihre Variablen ein Colinearitätsproblem aufweisen, das auf ein Problem mit dem Modell hinweist und nicht unbedingt mit numerischen Methoden gelöst werden sollte.

Wenn die Matrix aus numerischen Gründen nicht positiv semidefinit ist, gibt es einige Lösungen, über die hier nachgelesen werden kann

johnerisch
quelle

1

Die Annahme ist, dass das Modell ein lineares Mischmodell ist. Und für diesen Fall ist es nicht relevant, ein korrektes Modell für die Daten zu finden, sondern die Daten werden als Beispiel für einige Berechnungen angegeben. Nun besteht die Möglichkeit, dass Sie eine nicht positive semidefinite Matrix als Schätzung für die Kovarianz erhalten. Also, was tun, wenn ich die Kovarianz aus der normalverteilten Population herausfinden möchte, aus der die Daten stammen? Dass die Probe normal verteilt ist, ist die Annahme.

Klaus

1

Eine Möglichkeit wäre, die Matrix aus einer Eigenwertzerlegung zu berechnen. Jetzt gebe ich zu, dass ich nicht allzu viel über die Mathematik hinter diesen Prozessen weiß, aber meiner Recherche nach scheint es fruchtbar, diese Hilfedatei zu lesen:

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

und einige andere verwandte Befehle in R.

Schauen Sie sich auch 'nearPD' im Matrix-Paket an.

Tut mir leid, dass ich Ihnen nicht weiterhelfen konnte, aber ich hoffe, dass meine Suche Sie in die richtige Richtung bringen kann.

Frank P.
quelle

Hallo, danke für die Links. In Bezug auf die Eigenwertzerlegung hilft diese Zerlegung nicht, da Sie von dort komplexe Eigenwerte für die Quadratwurzelmatrix erhalten, aber ich brauche eine Matrix mit reell-Werten.

Klaus

1

Sie können die Ergebnisse aus der NearPD-Funktion im Matrix-Paket in R abrufen. Dies gibt Ihnen eine wirklich wertvolle Matrix zurück.

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
quelle

Für Benutzer von R .. ist dies in meiner Antwort möglicherweise keine schlechte Version der Lösungsmethode B für "Arme" (mit weniger Kontrolle).

Mark L. Stone

Ich bin damit einverstanden, dass dies nicht optimal ist, aber manchmal ist es der Trick.

Dr. Mike

Generieren Sie normalverteilte Zufallszahlen mit einer nicht positiv definierten Kovarianzmatrix

Antworten:

Algebraische Analyse der Kovarianzmatrix

Generierung multivariater Zufallswerte

Berechnungs- und Beispielcode