## Avisos sobre implementación:
## Podría usarse %o% en vez de diag(vector)%*%...
## as.matrix es necesario para %*%
Kij <- as.matrix (read.table ("datos.txt")) # frecuencias absolutas (p307)
k   <- sum (Kij)                            # tamaño muestral
Fij <- Kij / k                              # frecuencias relativas
Fi. <- apply (Fij, 1, sum)                  # marginales de filas
F.j <- apply (Fij, 2, sum)                  # marginales de columnas
Pi  <- diag (1/Fi.) %*% Fij                 # perfiles-fila    (p309)
Pj  <- Fij %*% diag (1/F.j)                 # perfiles-columna (p310)
## distancias (p311)
d2i <- function (i1,i2) sum ((Pi[i1,]-Pi[i2,])^2 / F.j)
d2j <- function (j1,j2) sum ((Pj[j1,]-Pj[j2,])^2 / Fi.)
Di  <- Pi %*% diag (1/sqrt(F.j))            # coordenadas euclídeas (p312)
Dj  <- diag (1/sqrt(Fi.)) %*% Pj
n   <- nrow (Kij)
p   <- ncol (Kij)
J   <- matrix (1, n, p)
Ci  <- Di - J %*% diag (sqrt (F.j))         # coordenadas centradas (p313)
Cj  <- Dj - diag (sqrt (Fi.)) %*% J
Ind <- diag(Fi.) %*% J %*% diag(F.j)        # independencia
X   <- (Fij - Ind) / sqrt(Ind)              # p313f11
T   <- t(X) %*% X                           # p313f12
Xa  <- Fij / sqrt(Ind)                      # X* de p314
Ta  <- t(Xa) %*% Xa                         # T* de p314
W   <- X %*% t(X)                           # p315
Wa  <- Xa %*% t(Xa)                         # W* de p315

auto.i  <- eigen (T)
auto.iA <- eigen (Ta)
auto.j  <- eigen (W)
auto.jA <- eigen (Wa)

landa.i  <- auto.i$values                    # autovalores
landa.iA <- auto.iA$values
landa.j  <- auto.j$values                    # autovalores
landa.jA <- auto.jA$values
## desde ahora, sólo el caso «.iA» (a partir de matriz «Ta»)
U <- auto.iA$vectors
V <- auto.jA$vector                          # ¡NO! por signos arbitrarios
V <- Xa %*% U %*% diag (1 / sqrt (landa.iA)) # mejor para gráfica conjunta
U <- U[,-1]                                  # elimino el inútil
V <- V[,-1]
landa <- landa.iA[-1]                        # elimino el inútil
pv <- 100 * landa / sum (landa)              # porcentajes de varianza
cbind (landa, pv, cumsum(pv))                # tabla 11 (p320)
psi.i <- Di %*% U                            # proyecciones (p314f14)
phi.j <- t(Dj) %*% V                         # p315f15
Ca.i  <- diag(Fi.) %*% psi.i^2 %*% diag(1/landa) # contribuciones absolutas (p318)
Ca.j  <- diag(F.j) %*% phi.j^2 %*% diag(1/landa) # p319
distancias.i <- apply (psi.i, 1, function (x) sum(x^2)) # distancias al origen
distancias.j <- apply (phi.j, 1, function (x) sum(x^2)) # distancias al origen
Cr.i  <- diag (1 / distancias.i) %*% psi.i^2 # contribuciones relativas
Cr.j  <- diag (1 / distancias.j) %*% phi.j^2 
## tabla 12 (p321)
print (round (cbind (F.j,
                     phi.j[,1:3],
                     Ca.j[,1:3],
                     Cr.j[,1:3]),
              3))
## tabla 13 (p321)
print (round (cbind (Fi.,
                     psi.i[,1:3],
                     Ca.i[,1:3],
                     Cr.i[,1:3]),
              3))
## gráfica de la figura 10 (p320)
## ojo: calcúlese V a partir de U, no directamente como autovector
coor.i <- psi.i[,1:2]                   # coordenadas bidimensionales de filas
coor.j <- phi.j[,1:2]                   # ídem columnas
plot (rbind (coor.i, coor.j),
      type = "n",                       # no dibuja los puntos
      xlab = "eje 1º", ylab = "eje 2º")
text (coor.i, labels = rownames (Kij), col = "red")
text (coor.j, labels = colnames (Kij), col = "blue")

## gráfica que considera la bondad de la representación
bondad.i <- apply (Cr.i[,1:2], 1 ,sum)
bondad.j <- apply (Cr.j[,1:2], 1 ,sum)
plot (rbind (coor.i, coor.j),
      type = "n",                       # no dibuja los puntos
      xlab = "eje 1º", ylab = "eje 2º")
text (coor.i, labels = rownames (Kij), col = rgb(1,0,0,bondad.i))
text (coor.j, labels = colnames (Kij), col = rgb(0,0,1,bondad.j))