R/axialnntsmanifoldnewtonestimationgradientstopsymmetric.R

In CircNNTSRaxial: Axial Data using NNTS Models

axialnntsmanifoldnewtonestimationgradientstopsymmetric<-function (data, M = 0, iter = 1000, gradientstop = 1e-10, pevalmu = 1000, initialpoint = FALSE, cinitial)
{
  data <- as.matrix(data)
  n <- length(data)
  R <- 1
  if (R != length(M))
    return("Error: M must correspond to a univariate variable distribution")
  cestnonsym <- axialnntsmanifoldnewtonestimationgradientstop(data, M = M, iter = iter, gradientstop = gradientstop, initialpoint = FALSE)
  cestsymmod <- Mod(cestnonsym$cestimates[, 2])
  mualt <- seq(0, pi, pi/pevalmu)
  resloglikmu <- rep(0, length(mualt))
  for (k in 1:length(mualt)) {
    resloglikmu[k] <- axialnntsloglik(data, cestsymmod * exp(seq(0,2*M,2) *
                                                               (0+1i) * (-mualt[k])), M)
  }
  muopt <- mualt[which.max(resloglikmu)]
  data <- data - muopt
  statisticsmatrix <- matrix(0, nrow = M + 1, ncol = n)
  for (k in 0:M)
    statisticsmatrix[k + 1, ] <- t(Conj(exp((0 + (0+1i)) * 2 * k * data)))
  c0 <- cestsymmod/sqrt(sum(Mod(cestsymmod)^2))
  eta <- matrix(0, nrow = M + 1, ncol = 1)
  for (k in 1:n)
    eta <- eta + as.vector(1/n) * as.vector(1/(t(Conj(c0)) %*% statisticsmatrix[, k])) * statisticsmatrix[, k]
  eta <- Mod(eta)
  eta <- eta - c0
  newtonmanifold <- (c0 + eta)
  newtonmanifold <- newtonmanifold/sqrt(sum(Mod(newtonmanifold)^2))
  newtonmanifoldprevious <- newtonmanifold
  for (j in 1:iter) {
    eta <- matrix(0, nrow = M + 1, ncol = 1)
    for (k in 1:n) {
      eta <- eta + as.vector(1/n) * as.vector(1/(t(Conj(newtonmanifold)) %*%
                                                   statisticsmatrix[, k])) * statisticsmatrix[,
                                                                                              k]
    }
    eta <- Mod(eta)
    eta <- eta - newtonmanifold
    newtonmanifold <- newtonmanifold + eta
    newtonmanifold <- newtonmanifold/sqrt(sum(Mod(newtonmanifold)^2))
    normsequence <- (sqrt(sum(Mod(newtonmanifold - newtonmanifoldprevious)^2)))
    newtonmanifoldprevious <- newtonmanifold
    data <- data + muopt
    cestsymmod <- Mod((1/sqrt(pi)) * newtonmanifoldprevious)
    mualt <- seq(0, pi, pi/pevalmu)
    resloglikmu <- rep(0, length(mualt))
    for (k in 1:length(mualt)) {
      resloglikmu[k] <- axialnntsloglik(data, cestsymmod * exp(seq(0,2*M,2) *
                                                                 (0+1i) * (-mualt[k])), M)
    }
    muopt <- mualt[which.max(resloglikmu)]
    data <- data - muopt
    if (normsequence < gradientstop)
      break
  }
  newtonmanifold <- newtonmanifold/sqrt(pi)
  data <- data + muopt
  newtonmanifold <- newtonmanifold * exp(seq(0,2*M,2) * (0+1i) * (-muopt))
  loglik <- axialnntsloglik(data, newtonmanifold, M)
  AIC <- -2 * loglik + 2 * (M + 1)
  BIC <- -2 * loglik + (M + 1) * log(n)
  gradnormerror <- normsequence
  cestimatesarray <- data.frame(cbind(0:M, newtonmanifold))
  cestimatesarray[, 1] <- as.integer(Re(as.matrix(cestimatesarray[,
                                                                  1])))
  names(cestimatesarray)[1] <- "k"
  names(cestimatesarray)[2] <- "cestimates"
  logliksym <- loglik
  logliknonsym <- cestnonsym$loglik
  loglikratio <- -2 * (logliksym - logliknonsym)
  loglikratiopvalue <- pchisq(loglikratio, df = M - 1, lower.tail = FALSE)
  res <- list(cestimatessym = cestimatesarray, mu = muopt,
              logliksym = loglik, AICsym = AIC, BICsym = BIC, gradnormerrorsym = gradnormerror,
              cestimatesnonsym = cestnonsym$cestimates, logliknonsym = cestnonsym$loglik,
              AICnonsym = cestnonsym$AIC, BICnonsym = cestnonsym$BIC,
              gradnormerrornonsym = cestnonsym$gradnormerror, loglikratioforsym = loglikratio,
              loglikratioforsympvalue = loglikratiopvalue)
  return(res)
}