1. Vorbetrachtung

Wir werden jetzt die Hauptkomponentenanalyse (PCA) an einem 3-dimensionalen Datensatz durchführen, also diesen derart auf zwei Dimensionen projizieren, dass die größtmögliche Varianz sichtbar wird.

2. Parameter for 3D-normalverteilte Daten

Wir simulieren dreidimensionale normalverteilte Daten mit den Mittelwerten \(\mu_x = 1\) , \(\mu_y = 3\) und \(\mu_z = 2\). Die Kovarianzmatrix (\(\Sigma\)) muss symmetrisch und positiv definit sein. Letzteres wird mit Hilfe der Funktion is.positive.definite aus dem Package matrixcalc überprüft.

mu = c(1, 3, 2)
Sigma <- matrix(c(3,   1.5,    1, 
                  1.5,   2,  0.5, 
                  1,   0.5,    1), byrow = TRUE, nrow = 3)
isSymmetric(Sigma)
## [1] TRUE
library(matrixcalc) 
is.positive.definite(Sigma, tol=1e-8)  
## [1] TRUE


3. Simulation der Daten

Wir setzen einen seed um die Datenerzeugung wiederholbar zu machen. Zur Erzeugung der Daten wird mvrnorm aus dem Package MASS benutzt.

set.seed(77)
library(MASS) 
data3 = mvrnorm(n = 500, mu = mu, Sigma = Sigma)
colnames(data3) = c("x","y","z")
head(data3)
##                x        y         z
## [1,]  2.92817247 1.716572 3.1399360
## [2,] -0.62420281 2.018428 0.3835038
## [3,] -0.69263473 3.479778 1.0662602
## [4,] -0.08737466 1.367032 0.5041806
## [5,]  0.64412600 2.979015 1.7780631
## [6,] -0.52558127 1.558227 0.6274831


4. Zentrieren und Skalieren der Daten

Die Hauptkomponentenanalyse beruht im Wesentlichen auf einer Analyse der Varianz von Datensätzen entlang verschiedener Raumrichtungen. Die Daten müssen daher zuvor standardisiert werden, d.h. die Daten müssen so transformiert werden, dass jede Variable den Mittelwert 0 und die Varianz 1 hat. Dies ermöglicht, vereinfacht ausgedrückt, eine “Gleichbehandlung” aller Variablen. (Sonst könnte man ja die Varianz einer physikalischen Größe einfach vergrößern, indem man eine andere Einheit wählt, z.B. Meter statt Kilometer.)

Wir definieren zunächst eine Funktion, welche die Standardisierung durchführen kann:

zscore <- function(vector) {
  new_vector = (vector - mean(vector))/sd(vector)  
  return(new_vector)
}

Wir wenden die Transformation nun auf die Daten an und überprüfen das Ergebnis:

dfz = apply(data3, 2, zscore)
colMeans(dfz)  
##             x             y             z 
## -4.683753e-20  5.145797e-17  2.427225e-17
apply(dfz, 2, sd) 
## x y z 
## 1 1 1

Die Mittelwerte sind nahe genug bei Null, Rundungsfehler können nicht vermieden werden.


5. Daten plotten

Zum Erstellen interaktiver Grafik verwenden wir einige Funktionen aus dem Package rgl:

par3d(windowRect = c(819, 88, 1492, 706))  
plot3d(dfz[,1], dfz[,2], dfz[,3], aspect = c(1,1,1), col = "blue", size = 3, xlab="x", ylab = "y", zlab = "z")
rglwidget()