-
Notifications
You must be signed in to change notification settings - Fork 17
/
em-pca.Rmd
187 lines (132 loc) · 4.07 KB
/
em-pca.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
## PCA
The following is an EM algorithm for principal components analysis. See Murphy,
2012 Probabilistic Machine Learning 12.2.5. Some of the constructed object is
based on output from <span class="func" style = "">pca</span> function used below.
### Data Setup
The `state.x77` is from base R, which includes various state demographics. We
will first standardize the data.
```{r em-pca-setup}
library(tidyverse)
X = scale(state.x77)
```
### Function
The estimating function. Note that it uses <span class="func" style =
"">orth</span> from <span class="pack" style = "">pracma</span>, but I show the
core of the underlying code if you don't want to install it.
```{r em-pca-func}
# orth <- function(M) {
# svdM = svd(M)
# U = svdM$u
# s = svdM$d
# tol = max(dim(M)) * max(s) * .Machine$double.eps
# r = sum(s > tol)
#
# U[, 1:r, drop = FALSE]
# }
em_pca <- function(
X,
n_comp = 2,
tol = .00001,
maxits = 100,
showits = TRUE
) {
# Arguments
# X: numeric data
# n_comp: number of components
# tol = tolerance level
# maxits: maximum iterations
# showits: show iterations
# starting points and other initializations
N = nrow(X)
D = ncol(X)
L = n_comp
Xt = t(X)
Z = t(replicate(L, rnorm(N))) # latent variables
W = replicate(L, rnorm(D)) # loadings
it = 0
converged = FALSE
if (showits)
cat(paste("Iterations of EM:", "\n"))
# while no convergence and we haven't reached our max iterations do this stuff
while ((!converged) & (it < maxits)) {
Z_old = Z # create 'old' values for comparison
Z = solve(t(W)%*%W) %*% crossprod(W, Xt) # E
W = Xt%*%t(Z) %*% solve(tcrossprod(Z)) # M
it = it + 1
# if showits, show first and every 5th iteration
if (showits & (it == 1 | it%%5 == 0))
cat(paste(format(it), "...", "\n", sep = ""))
converged = max(abs(Z_old-Z)) <= tol
}
# calculate reconstruction error
Xrecon_em = W %*% Z
reconerr = sum((Xrecon_em - t(X))^2)
# orthogonalize
W = pracma::orth(W) # for orthonormal basis of W; pcaMethods package has also
evs = eigen(cov(X %*% W))
evals = evs$values
evecs = evs$vectors
W = W %*% evecs
Z = X %*% W
if (showits) # Show last iteration
cat(paste0(format(it), "...", "\n"))
list(
scores = Z,
loadings = W,
reconerr = reconerr,
Xrecon_em = t(Xrecon_em)
)
}
```
### Estimation
```{r em-pca-est}
fit_em = em_pca(
X = X,
n_comp = 2,
tol = 1e-12,
maxit = 1000
)
str(fit_em) # examine results
```
### Comparison
Extract reconstructed values and loadings for comparison.
```{r em-pca-extract}
Xrecon_em = fit_em$Xrecon_em
loadings_em = fit_em$loadings
scores_em = fit_em$scores
```
Compare results to output from <span class="pack" style = "">pcaMethods</span>, which also has probabilistic PCA (demonstrated next). Note that the signs for loadings/scores may be different in sign, but otherwise should be comparable.
```{r em-pca-compare}
library(pcaMethods) # install via BiocManager::install("pcaMethods")
fit_pcam = pca(
X,
nPcs = 2,
method = 'svd',
scale = 'none',
center = FALSE
)
loadings_pcam = loadings(fit_pcam)
scores_pcam = scores(fit_pcam)
```
Compare loadings and scores.
```{r em-pca-compare-loadings}
sum((abs(loadings_pcam) - abs(loadings_em))^2)
cbind(scores_pcam, data.frame(EM = scores_em)) %>%
head()
```
Calculate mean squared reconstruction error and compare.
```{r em-pca-compare-recon}
Xrecon_pcam = scores_pcam %*% t(loadings_pcam)
mean((Xrecon_em - X)^2)
mean((Xrecon_pcam - X)^2)
mean(abs(Xrecon_pcam - Xrecon_em))
```
### Visualize
```{r em-pca-vis}
# qplot(Xrecon_pcam[,1], X[,1])
# qplot(Xrecon_pcam[,2], X[,2])
qplot(Xrecon_em[,1], Xrecon_pcam[,1])
```
### Source
Original code available at
https://github.com/m-clark/Miscellaneous-R-Code/blob/master/ModelFitting/EM%20Examples/EM%20for%20pca.R