2013-05-23 12 views
5

Ich frage mich, wie man Proben in Matlab, wo ich Präzision Matrix und Mittelwert als Eingabe Argument zeichnen.Wie zeichne ich Samples von multivariaten Gaußverteilung Verteilung durch Präzision in Matlab

Ich weiß, mvnrnd ist ein typischer Weg, aber es erfordert die Kovarianzmatrix (d. H. Umgekehrt Präzision)) als Argument.

Ich habe nur Präzision Matrix, und aufgrund der Rechen Problem, ich kann nicht meine Präzision Matrix invertieren, da es zu lange dauern wird (meine Dimension etwa 2000 * 2000)

Antwort

6

Gute Frage. Beachten Sie, dass Sie anhand einer Stichprobe aus der Standardnormalverteilung anhand einer in the relevant Wikipedia article beschriebenen Prozedur Stichproben aus einer multivariaten Normalverteilung generieren können.

Im Grunde ist dies läuft darauf hinaus, nach unten, wo der Bewertung A*z + muz ein Vektor von unabhängigen Zufallsvariablen von der Standardnormalverteilung abgetastet wird, ist ein Vektor von mu Mitteln und A*A' = Sigma die Kovarianzmatrix. Da Sie die Umkehrung der letzteren Quantität haben, d. H. inv(Sigma), können Sie wahrscheinlich eine Cholesky-Zerlegung durchführen (siehe chol), um die Inverse von A zu bestimmen. Sie müssen dann A * z auswerten. Wenn Sie nur inv(A) kennen, können Sie dies immer noch tun, ohne eine Matrixinvertierung durchzuführen, indem Sie stattdessen ein lineares System lösen (z. B. über den Backslash-Operator).

Die Cholesky-Zerlegung könnte für Sie immer noch problematisch sein, aber ich hoffe, das hilft.

+0

Können Sie die Mathematik dieses Teils erweitern, "Wenn Sie nur inv (A) wissen, kann dies immer noch ohne eine Matrixinverse durchgeführt werden, indem stattdessen ein lineares System gelöst wird (zB über den Backslash-Operator)" in anderen Sprachen nachgeahmt werden? – NewNameStat

0

Wenn Sie von N (μ, Q -1) und Q ist nur zur Probe verfügbar sind, können Sie die Cholesky-Faktorisierung von Q nehmen, L, so dass LL T = Q. Als nächstes nehme man das Inverse von L T, L -T, und teste Z aus einer Standardnormalverteilung N (0, I).

Anbetracht dessen, dass L -T ist eine obere Dreiecks dxd Matrix und Z ein d-dimensionalen Spaltenvektor, μ + L -T Z als N verteilt werden (μ, Q -1) .

Wenn Sie vermeiden möchten, die Inverse von L zu verwenden, können Sie stattdessen das Dreieckssystem der Gleichungen L T v = Z durch Rücksubstitution lösen. μ + v wird dann als N verteilt (μ, Q -1).

einigen anschaulichen Matlab Code:

% einer 2x2 Kovarianzmatrix und einen Durchschnittsvektor machen

covm = [3 0,4 * (sqrt (3 * 7)); 0.4 * (sqrt (3 * 7)) 7];

mu = [100; 2];

% Holen der Präzision Matrix

Q = inv (covm);

% die Cholesky-Zerlegung von Q nehmen (Chol in Matlab liefert bereits das obere Dreieck factor)

L = Chol (Q);

% Ziehen 2000 Proben aus einer Standard bivariate Normalverteilung

Z = normrnd (0,1, [2, 2000]); Die mittlere

X = repmat (mu, 1, 2000) + L \ Z

% lösen das System und fügen;

% überprüft das Ergebnis

Mittelwert (X ')

var (X')

corrcoef (X ')

% zur Probennahme aus der Kovarianzmatrix vergleichen

Y = mvrnd (mu, covm, 2000) ';

Mittelwert (Y ')

var (Y')

corrcoef (Y ')

scatter (X (1, ​​:), X (2 :),' b‘)

Halt auf

Streuung (Y (1, :), Y (2 :), 'r')

für mehr Effizienz, ich denke, Sie som suchen Ein Paket, das dreieckige Systeme effizient löst.