Rekonstruktion und Outlier
Uwe Menzel
2023-02-25
Rekonstruktion
Rekonstruktion ist die Umkehrung der Hauptkomponentenanalyse: aus den PCA-Daten werden die ursprünglichen Daten (fast genau) wiederhergestellt. Wichtig ist das fast genau: zur Rekonstruktion werden nur solche Hauptkomponenten benutzt, die einen wesentlichen Anteil der Varianz tragen. Auf diese Weise kann Rauschen aus Datensätzen entfernt werden, außerdem können anomale Daten (Outlier) identifiziert werden.
Beispieldaten: iris
Zur Demonstration der Rekonstruktion benutzen wir wieder die iris-Daten:
data(iris)
head(iris, 3) ## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosaDiesen Datensatz hatten wir schon hier näher beschrieben. Wir extrahieren den numerischen Anteil der iris-Daten:
X = as.matrix(iris[,-5])
head(X, 3)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## [1,] 5.1 3.5 1.4 0.2
## [2,] 4.9 3.0 1.4 0.2
## [3,] 4.7 3.2 1.3 0.2Hauptkomponentenanalyse
Vor der Hauptkomponentenanalyse skalieren wir wie immer die Daten:
zscore <- function(vector) (vector - mean(vector))/sd(vector) # center/scale to zero mean and sd = 1
X = apply(X, 2, zscore)
apply(X, 2, mean) # should be 0, at least very close to it## Sepal.Length Sepal.Width Petal.Length Petal.Width
## -4.484318e-16 2.034094e-16 -2.895326e-17 -2.989362e-17apply(X, 2, sd) # should be 1## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 1 1 1Wir konstruieren die Drehmatrix aus den Eigenvektoren der Kovarianzmatrix, wie hier bereits näher beschrieben:
# cov(X) # (we see in the main diagonal that the variance of all variables is 1)
D = eigen(cov(X))$vectors # matrix of eigenvectors = rotation matrix
colnames(D) = c("PC1", "PC2", "PC3", "PC4")
rownames(D) = colnames(X)
D## PC1 PC2 PC3 PC4
## Sepal.Length 0.5210659 -0.37741762 0.7195664 0.2612863
## Sepal.Width -0.2693474 -0.92329566 -0.2443818 -0.1235096
## Petal.Length 0.5804131 -0.02449161 -0.1421264 -0.8014492
## Petal.Width 0.5648565 -0.06694199 -0.6342727 0.5235971Im Abschnitt Theorie wurde bereits beschrieben, wie die Koordinaten in der PC-Ebene (\(U\)) aus den Originaldaten berechnet (\(X\)) werden:
\[ \boldsymbol{U} = \boldsymbol{X} \cdot \boldsymbol{D} \hspace{1cm} (1) \]
Wir führen die Matritzen-Multiplikation in R aus:
# dim(X) # dimensions must be compatible
# dim(D)
U = X %*% D # rotate into PC space
head(U, 3)## PC1 PC2 PC3 PC4
## [1,] -2.257141 -0.4784238 0.1272796 0.02408751
## [2,] -2.074013 0.6718827 0.2338255 0.10266284
## [3,] -2.356335 0.3407664 -0.0440539 0.02828231dim(U) ## [1] 150 4Wir überprüfen unser Ergebnis lieber noch einmal. Wenn wir richtig gerechnet haben, dann sollte das mit dem Ergebnis von prcomp übereinstimmen:
pca_results <- prcomp(X, center = TRUE, scale = TRUE)
head(pca_results$x, 3)## PC1 PC2 PC3 PC4
## [1,] -2.257141 -0.4784238 0.1272796 0.02408751
## [2,] -2.074013 0.6718827 0.2338255 0.10266284
## [3,] -2.356335 0.3407664 -0.0440539 0.02828231Rekonstruktion mit allen Hauptkomponenten
Als Übung bzw. zur Probe führen wir die Rekonstruktion zunächst mit allen Hauptkomponenten durch. Damit sollte sich der Originaldatensatz vollständig wieder herstellen lassen.
Wir müssen im Prinzip Gleichung (1) nach \(X\) umstellen. Dazu multiplizieren wir die Gleichung von rechts mit \(D^{-1}\). Der Term \(D \cdot D^{-1}\) verschwindet (wird zur Einheitsmatrix):
\[ \boldsymbol{U} \cdot \boldsymbol{D^{-1}} = \boldsymbol{X} \hspace{1cm} (2) \]
Die inverse Matrix zu \(D\), also \(D^{-1}\), könnte z.B. mit der Funktion Inverse aus dem package matlib berechnet werden.
Jedoch brauchen wir das nicht tun. Die Matrix \(D\) wurde aus den Eigenvektoren der Kovarianzmatrix gebildet. Die Eigenvektoren zu verschiedenen Eigenwerten sind zueinander orthogonal, womit auch die Drehmatrix eine orthogonale Matrix ist. Für orthogonale Matritzen ist die inverse Matrix gleich der transponierten Matrix:
\[ D^{-1} = D^T \]
Gleichung (2) vereinfacht sich also zu
\[ \boldsymbol{X} = \boldsymbol{U} \cdot \boldsymbol{D^{T}} \hspace{1cm} (3) \]
Dies ist in der Tat eine Vereinfachung, denn die Transponierte ist natürlich viel einfacher zu bekommen als die Inverse. Gleichung (3) gibt nun also die Formel zur vollständigen Rekonstruktion der Original-Daten an.
Wir führen das gleich in R aus:
# dim(U) # 150 4
# dim(t(D)) # 4 4
X1 = U %*% t(D) # back transformed Zur Probe vergleichen wir unsere rücktransformierten Daten mit den Originaldaten, indem wir die elementweisen Differenzen beider Matritzen überprüfen:
all(abs(X - X1) < 1e-14)## [1] TRUEDie Unterschiede sind tätsächlich sehr klein (aufgrund der Rechenungenauigkeit können wir nicht verlangen, dass alle Differenzen genau Null sind). Wir können das natürlich auch grafisch darstellen, indem wir z.B. die ersten zwei Dimensionen beider Datensätze mit unterschiedlichen Farben und unterschiedlichen Symbolen plotten:
plot(X[,1:2], col = "blue", pch = 1, cex = 1.4)
points(X1[,1:2], col = "red", pch = 18)
Ebensogut könnte man die Grafik mit der 3. und 4. Komponente oder einer anderen beliebigen Kombination von 2 Komponenten machen.
Wir müssen allerdings beachten, dass \(X\) und \(X1\) nicht die iris-Daten darstellen, sondern die skalierten Daten. Wenn wir die iris-Daten rekonstruieren wollen, müssen wir die Skalierung rückgängig machen.
Die Standardisierung erfolgte für jede Spalte der Originaldaten (\(X^\prime\)) so (\(X\) sind die standardisierten Daten):
\[ X = \frac{X^\prime - mean(X^\prime)}{sd(X^\prime)} \hspace{1cm} (4) \]
Hier sind \(mean(X^\prime)\) und \(sd(X^\prime)\) die Mittelwerte und die Standardabweichungen der Spalten der Originaldaten. Die Rücktransformation lautet also:
\[ X^\prime = X \cdot sd(X^\prime) + mean(X^\prime) \hspace{1cm} (5) \]
Das führen wir sofort aus:
mean_iris = colMeans(iris[,-5]) # original mean values
sd_iris = apply(iris[,-5], 2, sd) # original standard deviations
iris2 = matrix(NA, nrow = nrow(iris), ncol = ncol(iris[,-5])) # backtransformed data
for (i in 1:4) { # loop over columns
iris2[, i] = X1[, i] * sd_iris[i] + mean_iris[i] # reverse zscore
}
colnames(iris2) = colnames(iris[,-5])
head(iris2, 3)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## [1,] 5.1 3.5 1.4 0.2
## [2,] 4.9 3.0 1.4 0.2
## [3,] 4.7 3.2 1.3 0.2Wie vergleichen die rücktransformierten Daten iris2 mit
den Originaldaten iris[,-5] mit Hilfe einer grafischen
Darstellung der ersten beiden Spalten:
plot(iris[,1:2], col = "blue", pch = 1, cex = 1.4)
points(iris2[,1:2], col = "red", pch = 18)
Das stimmt also.
Die Mittelwerte und die Standardabweichungen der Originaldaten und der rekonstruierten Daten stimmen überein:
colMeans(iris[,-5])## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 5.843333 3.057333 3.758000 1.199333colMeans(iris2)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 5.843333 3.057333 3.758000 1.199333apply(iris[,-5], 2, sd)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.8280661 0.4358663 1.7652982 0.7622377apply(iris2, 2, sd)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.8280661 0.4358663 1.7652982 0.7622377Rekonstruktion mit nur 2 Hauptkomponenten
Bei den iris-Daten wurde offensichtlich, dass ein großer
Teil der Varianz in den ersten beiden Hauptkomponenten liegt, siehe z.
B. den Scree-Plot auf dieser
Seite. Es liegt die Vermutung nahe, dass es sich bei der 3. und 4.
Hauptkomponente nur um Rauschen handelt.
Deshalb können wir probieren, zur Rekonstruktion nur die ersten
beiden Hauptkomponenten zu benutzen.
Dies realisieren wir, indem wir auch von der Drehmatrix nur 2 Spalten
benutzen. In R sieht das so aus:
# dim(U[,1:2]) # 150 2
# dim(t(D[,1:2])) # 2 4
X2 = U[,1:2] %*% t(D[,1:2]) # partielle Rekonstruktion
dim(X2) ## [1] 150 4Allerdings ist \(X2\) jetzt nicht mehr korrekt skaliert, da die Standardabweichung nicht mehr 1 ist (währenddessen die Mittelwerte noch korrekt sind):
colMeans(X2)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## -1.591320e-16 8.974303e-17 -1.753921e-16 -1.703498e-16sd_X2 = apply(X2, 2, sd)
sd_X2## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.9605200 0.9954493 0.9918316 0.9670989Das korrigieren wir, indem wir jede Spalte von \(X2\) durch ihre aktuelle Standardabweichung teilen:
for (i in 1:4) {
X2[,i] = X2[,i]/sd_X2[i] # divide by standard deviation
}
sd_X2 = apply(X2, 2, sd)
sd_X2 # must be 1 ## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 1 1 1Die Daten sind jetzt nicht genau (aber ungefähr) rekonstruiert. Wir vergleichen \(X\) und \(X2\):
head(X,3) # standardized original data## Sepal.Length Sepal.Width Petal.Length Petal.Width
## [1,] -0.8976739 1.0156020 -1.335752 -1.311052
## [2,] -1.1392005 -0.1315388 -1.335752 -1.311052
## [3,] -1.3807271 0.3273175 -1.392399 -1.311052head(X2,3) # after partial reconstruction, standardized## Sepal.Length Sepal.Width Petal.Length Petal.Width
## [1,] -1.036474 1.0544805 -1.309050 -1.285219
## [2,] -1.389120 -0.0619984 -1.230289 -1.257883
## [3,] -1.412170 0.3215078 -1.387326 -1.399860Dies sind jeweils die standardisierten Daten.
Um zu den unskalierten Daten zu gelangen, müssen wir wieder die
Skalierung rückgängig machen:
iris3 = matrix(NA, nrow = nrow(iris), ncol = ncol(iris[,-5]))
for (i in 1:4) {
iris3[, i] = X2[, i] * sd_iris[i] + mean_iris[i] # # reverse zscore
}
colnames(iris3) = colnames(iris[,-5])
head(iris3, 3)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## [1,] 4.985065 3.516946 1.447137 0.2196907
## [2,] 4.693050 3.030310 1.586173 0.2405278
## [3,] 4.673963 3.197468 1.308956 0.1323075Die Mittelwerte und die Standardabweichungen der Originaldaten und der rekonstruierten Daten sollten wieder übereinstimmen:
colMeans(iris[,-5])## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 5.843333 3.057333 3.758000 1.199333colMeans(iris3)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 5.843333 3.057333 3.758000 1.199333apply(iris[,-5], 2, sd)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.8280661 0.4358663 1.7652982 0.7622377apply(iris3, 2, sd)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.8280661 0.4358663 1.7652982 0.7622377Wie vergleichen wieder die Originaldaten und die rücktransformierten Daten. Diesmal können wir natürlich nicht erwarten, dass sie exakt übereinstimmen.
plot(iris[,1:2], col = "blue", pch = 1, cex = 1.4)
points(iris3[,1:2], col = "red", pch = 18)
Die rücktransformierten Daten scheinen kompakter zu sein, was auf das Fehlen von Rauschen zurückgeführt werden könnte.
Reconstruction Error
Der Reconstruction Error wird für jede Observation berechnet. In diesem Kontext ist er die Summe der quadratischen Differenzen zwischen den rekonstruierten und den Originaldaten:
rec_error = rowSums((iris[,-5] - iris3)^2)
names(rec_error) = 1:nrow(iris)
summary(rec_error)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.002096 0.037254 0.086025 0.144046 0.197888 0.737684plot(rec_error, col = "red", pch = 18, main = "Reconstruction Error", font.main = 1, ylab = "RE")
Einige Observationen haben einen ungewöhnlich hohen Reconstruction Error. Es kann vermutet werden, dass es sich dabei um Daten-Anomalien (Outlier) handelt, etwa durch fehlerhafte Messungen. Wir schauen uns die Observationen mit den fünf größten Reconstruction Errors an:
outlier = head(sort(rec_error, decreasing = TRUE), 5)
outlier## 107 101 115 137 149
## 0.7376841 0.7197136 0.6694112 0.5920851 0.5788424Diese Beobachtungen haben auch relativ große Varianz in der 3.
Hauptkomponente, welche wir zuvor aufgrund des
Scree-Plots
als Rauschen eingestuft hatten:
summary(U[,3])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.00204 -0.19386 -0.02468 0.00000 0.25820 0.85456ind = as.integer(names(outlier))
U[ind, 3] # close to the edge ## [1] -0.9835981 -1.0020441 -1.0005177 -0.9426957 -0.9302787Dadurch wird die Vermutung erhärtet, dass es sich bei diesen Beobachtungen um Outlier handeln könnte.
Note that all code snippets are designed for comprehensibility, not for efficiency!