R/axialnntsmanifoldnewtonestimationgradientstopknownmu.R

In CircNNTSRaxial: Axial Data using NNTS Models

axialnntsmanifoldnewtonestimationgradientstopknownmu<-function(data, muknown=0, M = 0, iter = 1000, initialpoint = FALSE, cinitial,gradientstop=1e-10)
{
  data <- as.matrix(data)
  n <- length(data)
  R <- 1
  if (R != length(M))
    return("Error: M must correspond to a univariate variable distribution")
  cestmuunknown <- axialnntsmanifoldnewtonestimationgradientstop(data, M = M, iter = iter, gradientstop = gradientstop, initialpoint = FALSE)
  cestmuknown <- cestmuunknown$cestimates[, 2]
  cestmuknown[1:2] <- Mod(cestmuknown[1:2])
  muopt <- muknown
  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 <- cestmuknown/sqrt(sum(Mod(cestmuknown)^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[1:2] <- Mod(eta[1:2])
  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[1:2] <- Mod(eta[1:2])
    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
    cestmuknown <- Mod((1/sqrt(pi)) * newtonmanifoldprevious)
    muopt <- muknown
    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 )
  BIC <- -2 * loglik + (M ) * 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"
  loglikmuknown <- loglik
  loglikmuunknown <- cestmuunknown$loglik
  loglikratio <- -2 * (loglikmuknown - loglikmuunknown)
  loglikratiopvalue <- pchisq(loglikratio, df = 1, lower.tail = FALSE)
  res <- list(cestimatesmuknown = cestimatesarray, muknown = muopt,
              loglikmuknown = loglik, AICmuknown = AIC, BICmuknown = BIC, gradnormerrormuknown = gradnormerror,
              cestimatesmuunknown = cestmuunknown$cestimates, loglikmuunknown = cestmuunknown$loglik,
              AICmuunknown = cestmuunknown$AIC, BICmuunknown = cestmuunknown$BIC,
              gradnormerrormuunknown = cestmuunknown$gradnormerror, loglikratioformuknown = loglikratio,
              loglikratioformuknownpvalue = loglikratiopvalue)
  return(res)
}