R/axialnntssimulation.R

In CircNNTSRaxial: Axial Data using NNTS Models

axialnntssimulation<-function (nsim = 1, cpars = 1/pi, M = 0)
{
  size <- length(cpars)
  if (size != M + 1)
    return("Length of cpars must be equal to M+1")
  if (abs(sum(Mod(cpars)^2) - 1/(pi)) > 1e-10)
    return("Sum of the squared norms of components greater than condition")
  res <- rep(0, nsim)
  conteo <- 1
  for (k in 1:nsim) {
    U1 <- runif(1, 0, pi)
    U2 <- runif(1, 0, (M + 1)/(pi))
    while (U2 > axialnntsdensity(U1, cpars, M)) {
      U1 <- runif(1, 0, pi)
      U2 <- runif(1, 0, (M + 1)/(pi))
      conteo <- conteo + 1
    }
    res[k] <- U1
  }
  resf <- list(simulations = res, conteo = conteo)
  resf
}