R/mlsq.R

In scpropreg: Simplicially Constrained Regression Models for Proportions

## sum( ( log(y) - log(y_hat) )^2 )
# Simple Multiplicative Fixed Point
mlsq <- function(Y, x, tol = 1e-8, maxit = 50000, alpha = 0.01) {
  p <- dim(x)[2]
  tx <- t(x)
  lY <- log(Y)

  d <- dim(Y)[2]
  B <- matrix(0, p, d)
  obj <- iterations <- numeric(d)

  for ( j in 1:d ) {
    ly <- lY[, j]
    be <- rep(1/p, p)
    for (iter in 1:maxit) {
      y_hat <- drop(x %*% be)
      grad <-  -2 * tx %*% ( ( ly - log(y_hat) ) / y_hat )
      grad_centered <- grad - mean(grad)
      be_new <- be * exp(-alpha * grad_centered)
      be_new <- be_new / sum(be_new)
      if ( max( abs(be_new - be) ) < tol ) {
        be <- be_new
        break
      }
      be <- be_new
    }

    y_hat_final <- x %*% be
    obj[j] <- sum( ( ly - log(y_hat_final) )^2 )
    B[, j] <- round(be, 12)
    iterations[j] <- iter
  }

  list(coefficients = round(B, 12), value = obj, iterations = iterations)
}