library('mvtnorm')
library('dplyr')

set.seed(42)
## simulate two-class problem with two inputs
sig <- matrix(c(1.0, 0.65,
                0.65, 1.0),2,2)
x <- rbind(
  rmvnorm(100, mean=c(2,2), sigma = sig),
  rmvnorm(100, mean=c(-5,-2), sigma = sig))
y <- rep(c(1,2),c(100,100))

## how does centering affect inputs
cx <- scale(x, center=TRUE, scale=FALSE)
attributes(cx)
cols <- c('#E69F00', '#009E73')
par(mfrow=c(1,2))
plot(x[,1], x[,2],
     xlab=expression(x[1]),
     ylab=expression(x[2]),
     pch=20, col=cols[y],
     main='Original')
abline(h=0, v=0, lty=3)
plot(cx[,1], cx[,2],
     xlab=expression(x[1]),
     ylab=expression(x[2]),
     pch=20, col=cols[y],
     main='Centered')
abline(h=0, v=0, lty=3)

## compute principal components
x <- cx

## using prcomp
pc <- prcomp(x)$x
par(mfrow=c(1,2))
plot(x[,1], x[,2],
     xlab=expression(x[1]),
     ylab=expression(x[2]),
     pch=20, col=cols[y],
     main='Centered')
abline(h=0, v=0, lty=3)
plot(-pc[,1], -pc[,2],
     xlab=expression(PC[1]),
     ylab=expression(PC[2]),
     pch=20, col=cols[y],
     main='prcomp(x)$x')
abline(h=0, v=0, lty=3)

## using svd
pc <- x %*% svd(x)$v
par(mfrow=c(1,2))
plot(x[,1], x[,2],
     xlab=expression(x[1]),
     ylab=expression(x[2]),
     pch=20, col=cols[y],
     main='Centered')
abline(h=0, v=0, lty=3)
plot(-pc[,1], -pc[,2],
     xlab=expression(PC[1]),
     ylab=expression(PC[2]),
     pch=20, col=cols[y],
     main='x%*%svd(x)$v')
abline(h=0, v=0, lty=3)

## using eigen
pc <- x %*% eigen(cov(x))$vectors
par(mfrow=c(1,2))
plot(x[,1], x[,2],
     xlab=expression(x[1]),
     ylab=expression(x[2]),
     pch=20, col=cols[y],
     main='Centered')
abline(h=0, v=0, lty=3)
plot(-pc[,1], -pc[,2],
     xlab=expression(PC[1]),
     ylab=expression(PC[2]),
     pch=20, col=cols[y],
     main='x%*%eigen(cov(x))$vectors')
abline(h=0, v=0, lty=3)


## compute estimates for LDA

## estimate class means
mu_hat <- lapply(unique(y), function(cl) {
  colMeans(x[y == cl,])
})

## estimate pooled variance
sigma_hat <- lapply(unique(y), function(cl) {
  sig <- cov(x[y == cl,])
  sig * (20-1)/(60-3)
}) %>% Reduce('+', .)

## estimate pi (prior)
pi_hat <- as.list(rep(1/2,2))

## discriminant functions
delta <- function(x, k, sigmah=sigma_hat, muh=mu_hat, pih=pi_hat) {
  pik <- pih[[k]]
  muk <- muh[[k]]
  isg <- solve(sigmah)
  dmvnorm(x, mean=muk, sigma=sigmah, log=TRUE) + log(pik)
}

## plot centered x with class means
par(mfrow=c(1,2))
plot(x[,1], x[,2],
     xlab=expression(x[1]),
     ylab=expression(x[2]),
     pch=20, col=cols[y],
     main='Centered')
abline(h=0, v=0, lty=3)
points(mu_hat[[1]][1], mu_hat[[1]][2], pch=4, lwd=4)
points(mu_hat[[2]][1], mu_hat[[2]][2], pch=4, lwd=4)

## compute sphered x
svds <- svd(sigma_hat)
print(svds)
sx <- x %*% svds$v %*% diag(1/sqrt(svds$d))

## compute sphered mu
smu_hat <- lapply(mu_hat, function(mu) 
  mu %*% svds$v %*% diag(1/sqrt(svds$d)))

plot(sx[,1], sx[,2],
     xlab=expression(s[1]),
     ylab=expression(s[2]),
     pch=20, col=cols[y],
     main='Sphered')
abline(h=0, v=0, lty=3)
points(smu_hat[[1]][1], smu_hat[[1]][2], pch=4, lwd=4)
points(smu_hat[[2]][1], smu_hat[[2]][2], pch=4, lwd=4)

## plot line that connects two means
par(mfrow=c(1,1))
plot(sx[,1], sx[,2],
     xlab=expression(s[1]),
     ylab=expression(s[2]),
     pch=20, col=cols[y],
     main='Sphered')
abline(h=0, v=0, lty=3)
points(smu_hat[[1]][1], smu_hat[[1]][2], pch=4, lwd=4)
points(smu_hat[[2]][1], smu_hat[[2]][2], pch=4, lwd=4)

smu_hat_mx <- do.call(rbind, smu_hat)
b_plot <- diff(smu_hat_mx[,2])/diff(smu_hat_mx[,1]) 
a_plot <- smu_hat_mx[1,2] - b_plot * smu_hat_mx[1,1]
abline(b=b_plot, a=a_plot)

## compute PCA on sphered centers
svdb <- svd(cov(smu_hat_mx))
svdb$v <- svdb$v
print(svdb)
print(eigen(cov(smu_hat_mx)))

## compute canonical variables and centers
cx <- sx %*% svdb$v[,1]
cmu_hat_mx <- smu_hat_mx %*% svdb$v

## plot how first canonical variable looks in sphered space
scx <- cx[,1,drop=F] %*% t(svdb$v[,1,drop=F])
par(mfrow=c(1,1))
plot(sx[,1], sx[,2],
     xlab=expression(s[1]),
     ylab=expression(s[2]),
     main='Sphered', 
     type='n')
points(sx[,1], sx[,2],
       pch=1, col=cols[y])
points(scx[,1], scx[,2],
       pch=20, col=cols[y])
abline(h=0, v=0, lty=3)
abline(b=b_plot, a=a_plot)
points(smu_hat[[1]][1], smu_hat[[1]][2], pch=4, lwd=4)
points(smu_hat[[2]][1], smu_hat[[2]][2], pch=4, lwd=4)

## plot canonical variables
par(mfrow=c(1,1))
plot(sx[,1], sx[,2],
     xlab=expression(c[1]),
     ylab=expression(c[2]),
     main='Canonical', 
     type='n')
points(cx[,1], cx[,2], #rep(0,nrow(cx)),
       pch=20, col=cols[y])
abline(h=0, v=0, lty=3)
points(cmu_hat_mx[1,1], 0, pch=4, lwd=4)
points(cmu_hat_mx[2,1], 0, pch=4, lwd=4)

## plot just first canonical variable
par(mfrow=c(1,1))
plot(sx[,1], sx[,2],
     xlab=expression(c[1]),
     ylab=expression(c[2]),
     main='Canonical', 
     type='n')
points(cx[,1], rep(0,nrow(cx)),
       pch=20, col=cols[y])
abline(h=0, v=0, lty=3)
points(cmu_hat_mx[1,1], 0, pch=4, lwd=4)
points(cmu_hat_mx[2,1], 0, pch=4, lwd=4)


## transform a new inputs into to canonical variable space
x0 <- matrix(c(-1.6, 1.7,
                1.5, 2.1),
             2,2,byrow=TRUE)

## centered
m0 <- t(t(x0) - attr(x, 'scaled:center'))

## plot centered x with class means
par(mfrow=c(1,3))
plot(x[,1], x[,2],
     xlab=expression(x[1]),
     ylab=expression(x[2]),
     pch=20, col=cols[y],
     main='Centered')
abline(h=0, v=0, lty=3)
points(mu_hat[[1]][1], mu_hat[[1]][2], pch=4, lwd=4)
points(mu_hat[[2]][1], mu_hat[[2]][2], pch=4, lwd=4)
points(m0[,1], m0[,2], pch=1:2, lwd=2, col='#AA4499')

## sphered
s0 <- m0 %*% svds$v %*% diag(1/sqrt(svds$d))

plot(sx[,1], sx[,2],
     xlab=expression(s[1]),
     ylab=expression(s[2]),
     pch=20, col=cols[y],
     main='Sphered')
abline(h=0, v=0, lty=3)
points(smu_hat[[1]][1], smu_hat[[1]][2], pch=4, lwd=4)
points(smu_hat[[2]][1], smu_hat[[2]][2], pch=4, lwd=4)
points(s0[,1], s0[,2], pch=1:2, lwd=2, col='#AA4499')


## canonical
c0 <- s0 %*% svdb$v

## plot just first canonical variable
plot(sx[,1], sx[,2],
     xlab=expression(c[1]),
     ylab=expression(c[2]),
     main='Canonical', 
     type='n')
points(cx[,1], rep(0,nrow(cx)),
       pch=20, col=cols[y])
abline(h=0, v=0, lty=3)
points(cmu_hat_mx[1,1], 0, pch=4, lwd=4)
points(cmu_hat_mx[2,1], 0, pch=4, lwd=4)
points(c0[,1], rep(0,nrow(c0)), pch=1:2, lwd=2, col='#AA4499')