Willkommen bei Stats by Randolph. Hier geht es zur Hauptseite mit weiteren Informationen und Inhalten.
Autor:innen dieser Seite: An den Inhalten dieser Seite haben mitgearbeitet: Markus Janczyk und Valentin Koob. Der Inhalt dieser Seite wird in der Lehre in den Studiengängen Psychologie von der AG Forschungsmethoden und Kognitive Psychologie an der Universität Bremen verwendet, steht aber allen Interessierten zur Verfügung. Rückmeldungen/Fehler/Vorschläge können gesendet werden an randolph@uni-bremen.
Versionshistory:
Im Haupttext zu Statistik II haben wir im Zuge der multivariaten
Statistik (Teil
11), die Mahalanobis-Distanz kennengelernt. Der schwierige Teil der
Mahalanobis-Transformation ist dabei die Bestimmung der Matrix \(\mathbf{S^{-\frac{1}{2}}}\), für die gelten
soll: \[(\mathbf{S^{-\frac{1}{2}}})'\cdot
\mathbf{S^{-\frac{1}{2}}}=\mathbf{S^{-1}}\] Dabei sind wir im
Haupttext nicht auf die genaue Berechnung eingegangen, sondern haben
lediglich die Funktion sqrtm() aus dem Paket
expm benutzt. Die genaue Berechnung wollen wir nun hier
nachliefern.
Etwas vereinfacht ist zur Herleitung die Jordan-Dekomposition der Ausgangspunkt, die besagt, dass jede symmetrische Matrix \(\mathbf{A}\) geschrieben werden kann als \[\mathbf{A}=\Gamma\Lambda\Gamma'\] wobei \(\Lambda\) eine Diagonalmatrix der Eigenwerte von \(\mathbf{A}\) ist (wir schreiben dann auch \(\Lambda=\text{diag}(\lambda_i)\)) und \(\Gamma\) eine Matrix ist, deren Spalten die (standardisierten) Eigenvektoren von \(\mathbf{A}\) sind. Wenn es sich zudem bei \(\mathbf{A}\) um eine nicht-singuläre, positiv semi-definite Matrix handelt, das heißt, eine Matrix die invertierbar ist (bzw. \(\det\mathbf{A}\neq 0\)) und zudem ausschließlich nicht-negative Eigenwerte hat (Eine Kovarianzmatrix hat in der Regel nur positive Eigenwerte und ein Eigenwert von Null kann nur dann resultieren, wenn zwei der Variablen perfekt korreliert sind.), dann ist jede Potenz der Matrix \(\mathbf{A^{\frac{r}{s}}}\) (für ganze Zahl \(r\) und \(s>0\)) definiert als \[ \mathbf{A^{\frac{r}{s}}}=\Gamma\Lambda^{\frac{r}{s}}\Gamma'\qquad \text{mit }\Lambda^{\frac{r}{s}}=\text{diag}(\lambda^{\frac{r}{s}}) \] Daraus ergibt sich für \(r=-1\) und \(s=2\) die Beziehung \[ \mathbf{A^{-\frac{1}{2}}}=\Gamma\Lambda^{-\frac{1}{2}}\Gamma'\qquad \text{mit }\Lambda^{-\frac{1}{2}}=\text{diag}\left(\frac{1}{\sqrt{\lambda_i}}\right) \] Diese Berechnung “von Hand” würde mit R z.B. so aussehen:
library(expm) # benötigt für sqrtm()
# Beispiel 1:
S <- rbind( c(5, 2), # Erstellen der Kovarianzmatrix S
c(2, 5))
Gamma <- eigen(S)$vectors # Matrix der Eigenvektoren von S
Lambda <- diag(1/sqrt(eigen(S)$values)) # Diagonal-Matrix aus den Eigenwerten von S
Gamma %*% Lambda %*% t(Gamma) # Berechnung gemäß Formel## [,1] [,2]
## [1,] 0.4776574 -0.0996929
## [2,] -0.0996929 0.4776574sqrtm(solve(S)) # zum Vergleich mit sqrtm()## [,1] [,2]
## [1,] 0.4776574 -0.0996929
## [2,] -0.0996929 0.4776574# Beispiel 2:
S <- rbind( c(5, 2, 3),
c(2, 4, 4),
c(3, 4, 5))
Gamma <- eigen(S)$vectors # Matrix der Eigenvektoren von S
Lambda <- diag(1/sqrt(eigen(S)$values)) # Diagonal-Matrix aus den Eigenwerten von S
Gamma %*% Lambda %*% t(Gamma) # Berechnung gemäß Formel## [,1] [,2] [,3]
## [1,] 0.53684565 0.02738559 -0.2106659
## [2,] 0.02738559 0.99067791 -0.5925712
## [3,] -0.21066588 -0.59257121 0.9684072sqrtm(solve(S)) # zum Vergleich mit sqrtm()## [,1] [,2] [,3]
## [1,] 0.53684565 0.02738559 -0.2106659
## [2,] 0.02738559 0.99067791 -0.5925712
## [3,] -0.21066588 -0.59257121 0.9684072