set.seed (12345) unit <- function (i, n) ifelse (i == 1:n, 1, 0) x <- matrix (rnorm (10), 5, 2) x <- apply (x, 2, function (z) z - mean (z)) e <- as.matrix (dist (x) ^ 2) y <- matrix (rnorm (20), 5, 4) y <- apply (y, 2, function (z) z - mean (z)) f <- as.matrix (dist (y) ^ 2) w <- matrix (rnorm (25) ^ 2, 5, 5) w <- w + t(w) diag (w) <- 0 v <- - 2 * sqrt(w) diag (v) <- - rowSums (v) cat ("sstress of x", formatC (sum (w * (f - e) ^ 2), width = 15, digits = 10, format = "f"), "\n") ww <- ff <- ee <- array(0, c(5, 5, 5, 5)) for (i in 1:5) for (j in 1:5) for (k in 1:5) for (l in 1 :5) { ww[i, j, k, l] <- sqrt(w[i, j] * w[k, l]) if ((i == k) && (j == l)) ff[i, j, k, l] <- f[i, j] else ff[i, j, k, l] <- rchisq (1, 1) ee[i, j, k, l] = sum ((x[i, ] - x[j, ]) * (x[k, ] - x[l, ])) } cat ("sstress_star of x direct", formatC (sum (ww * (ff - ee) ^ 2), width = 15, digits = 10, format = "f"), "\n") cat ("ssq augmented distances direct", formatC (sum (ww * ee ^ 2), width = 15, digits = 10, format = "f"), "\n") cat ("ssq augmented distances with V", formatC (sum ((t(x) %*% v %*% x) ^ 2), width = 15, digits = 10, format = "f"), "\n") c <- matrix (0, 5, 5) for (i in 1:5) for (j in 1:5) for (k in 1:5) for (l in 1:5) { a <- outer (unit (i, 5) - unit (j, 5), unit (k, 5) - unit (l, 5)) c <- c + ww[i, j, k, l] * ff[i, j, k, l] * a } cat ("cross product augmented direct", formatC (sum (ww * ff * ee), width = 15, digits = 10, format = "f"), "\n") cat ("cross product augmented with C", formatC (sum (diag (t(x) %*% c %*% x)), width = 15, digits = 10, format = "f"), "\n") cat ("sstress_star of x formula", formatC (sum (ww * ff ^ 2) - 2 * sum (diag (t(x) %*% c %*% x)) + sum ((t(x) %*% v %*% x) ^ 2), width = 15, digits = 10, format = "f"), "\n") h <- - 2 * (f - e) * w diag (h) <- -rowSums (h) b <- c <- matrix (0, 5, 5) for (i in 1:5) for (j in 1:5) for (k in 1:5) for (l in 1:5) { a <- outer (unit (i, 5) - unit (j, 5), unit (k, 5) - unit (l, 5)) if ((i == k) && (j == l)) ff[i, j, k, l] <- f[i, j] else ff[i, j, k, l] <- ee[i, j, k, l] c <- c + ww[i, j, k, l] * ff[i, j, k, l] * a b <- b + ww[i, j, k, l] * ee[i, j, k, l] * a } print(c) print(h+b) print(b) u <- v %*% x print (tcrossprod (u))