library (numDeriv) library (MASS) amalgm <- function (x, w = rep (1, length (x))) { dyn.load ("pava.so") n <- length (x) a <- rep (0, n) b <- rep (0, n) y <- rep (0, n) lf <- .Fortran ( "AMALGM", n = as.integer (n), x = as.double (x), w = as.double (w), a = as.double (a), b = as.double (b), y = as.double (y), tol = as.double(1e-15), ifault = as.integer(0) ) return (lf\$y) } isotone <- function (x, y, w = rep (1, length (x)), ties = "secondary") { f <- sort(unique(x)) g <- lapply(f, function (z) which(x == z)) n <- length (x) k <- length (f) if (ties == "secondary") { w <- sapply (g, length) h <- lapply (g, function (x) y[x]) m <- sapply (h, sum) / w r <- amalgm (m, w) s <- rep (0, n) for (i in 1:k) s[g[[i]]] <- r[i] } if (ties == "primary") { h <- lapply (g, function (x) y[x]) m <- rep (0, n) for (i in 1:k) { ii <- order (h[[i]]) g[[i]] <- g[[i]][ii] h[[i]] <- h[[i]][ii] } m <- unlist (h) r <- amalgm (m, w) s <- r[order (unlist (g))] } if (ties == "tertiary") { w <- sapply (g, length) h <- lapply (g, function (x) y[x]) m <- sapply (h, sum) / w r <- amalgm (m, w) s <- rep (0, n) for (i in 1:k) s[g[[i]]] <- y[g[[i]]] + (r[i] - m[i]) } return (s) } rStress <- function (x, w, delta, a, r) { n <- length (a) s <- 0 for (i in 1:n) { xax <- sum (x * (a[[i]] %*% x)) s <- s + w[i] * (delta[i] - xax ^ r) ^ 2 } return (s) } mdsDerivatives <- function (x, w, delta, a, r, numerical = FALSE) { m <- length (x) n <- length (a) b <- c <- s <- t <- matrix (0, m, m) for (i in 1:n) { xa <- drop (a[[i]] %*% x) xax <- sum (x * xa) b <- b + w[i] * delta [i] * (xax ^ (r - 1)) * a[[i]] c <- c + w[i] * (xax ^ (2 * r - 1)) * a[[i]] s <- s + w[i] * delta[i] * (xax ^ (r - 1)) * (a[[i]] + 2 * (r - 1) * outer(xa, xa) / xax) t <- t + w[i] * (xax ^ (2 * r - 1)) * (a[[i]] + 2 * (2 * r - 1) * outer(xa, xa) / xax) } gan <- -4 * r * drop ((b - c) %*% x) han <- -4 * r * (s - t) result <- list ( b = b, c = c, s = s, t = t, gan = gan, han = han ) if (numerical) { gnu <- grad ( rStress, x, w = w, delta = delta, a = a, r = r ) hnu <- hessian ( rStress, x, w = w, delta = delta, a = a, r = r ) result <- c (result, list (gnu = gnu, hnu = hnu)) } return (result) } torgerson <- function(delta, p = 2) { doubleCenter <- function(x) { n <- dim(x)[1] m <- dim(x)[2] s <- sum(x) / (n * m) xr <- rowSums(x) / m xc <- colSums(x) / n return((x - outer(xr, xc, "+")) + s) } z <- eigen(-doubleCenter((as.matrix (delta) ^ 2) / 2), symmetric = TRUE) v <- pmax(z\$values, 0) return(z\$vectors[, 1:p] %*% diag(sqrt(v[1:p]))) } u <- function (i, n) { return (ifelse (i == 1:n, 1, 0)) } e <- function (i, j, n) { d <- u (i, n) - u (j, n) return (outer (d, d)) } directSum <- function (x) { m <- length (x) nr <- sum (sapply (x, nrow)) nc <- sum (sapply (x, ncol)) z <- matrix (0, nr, nc) kr <- 0 kc <- 0 for (i in 1:m) { ir <- nrow (x[[i]]) ic <- ncol (x[[i]]) z[kr + (1:ir), kc + (1:ic)] <- x[[i]] kr <- kr + ir kc <- kc + ic } return (z) } repList <- function(x, n) { z <- list() for (i in 1:n) z <- c(z, list(x)) return(z) } makeA <- function (n, p = 2) { m <- n * (n - 1) / 2 a <- list() for (j in 1:(n - 1)) for (i in (j + 1):n) { d <- u (i, n) - u (j, n) e <- outer (d, d) a <- c(a, list (directSum (repList (e, p)))) } return (a) } newtonMe <- function (delta, xini = NULL, w = rep (1, length (delta)), p = 2, r = .5, eps = 1e-15, itmax = 1000, linearize = FALSE, nonmetric = FALSE, ties = "secondary", verbose = TRUE) { n <- nrow (as.matrix (delta)) dhat <- delta / sqrt (sum (delta ^ 2)) if (is.null (xini)) { xold <- as.vector (torgerson (dhat, p)) } else { xold <- xini } a <- makeA (n, p) dold <- sapply (a, function (u) sum (xold * (u %*% xold))) eold <- dold ^ r sold <- sum (w * (dhat - eold) ^ 2) itel <- 1 repeat { h <- mdsDerivatives (xold, w, dhat, a, r) if (linearize) { xnew <- drop (xold - ginv (4 * r * h\$t) %*% h\$gan) } else { xnew <- drop (xold - ginv (h\$han) %*% h\$gan) } dnew <- sapply (a, function (u) sum (xnew * (u %*% xnew))) enew <- dnew ^ r if (nonmetric) { dhat <- isotone (delta, enew, ties = ties) dhat <- dhat / sqrt (sum (dhat ^ 2)) } snew <- sum (w * (dhat - enew) ^ 2) if (verbose) { cat ( formatC (itel, width = 4, format = "d"), formatC ( sold, digits = 10, width = 13, format = "f" ), formatC ( snew, digits = 10, width = 13, format = "f" ), "\n" ) } if ((itel == itmax) || (abs(sold - snew) < eps)) break itel <- itel + 1 xold <- xnew dold <- dnew sold <- snew } return (list ( x = matrix (xnew, n, p), d = dnew, dhat = dhat, rstress = snew, g = h\$gan, h = h\$han, itel = itel )) }