ODA: Ocean Data Analysis package in R

  est_mov <- function(ang, np = 5, n.values = 360, nruns = 150,
                      main = "von-Mises Mixture Distribution", h.class = 36,
                      cex.size = 1.2){



    if(require(pacman)==FALSE){install.packages("pacman")}
    pacman::p_load(movMF,circular,NPCirc,Directional,dplyr,latex2exp)



    bic.mixvmf2 <- function (x, A = 3){

      x <- as.matrix(x)
      n <- nrow(x)
      p <- ncol(x)
      x <- x/sqrt(rowSums(x^2))
      BIC <- 1:A
      mod <- vmf(x)
      BIC[1] <- -2 * mod$loglik + p * log(n)

      for (vim in 2:A) {

        a <- mix.vmf(x, vim)
        BIC[vim] <- -2 * a$loglik + ((vim - 1) + vim * p) * log(n)

      }

      names(BIC) <- 1:A
      ina <- rep(1, A)
      ina[which.min(BIC)] <- 2
      BIC

    }



    ang <- na.omit(ang)
    rd <- rad(ang)
    rad.comp <- cbind( cos(rd), sin(rd) )
    vm.bic <- bic.mixvmf2( rad.comp, np )
    onp <- which.min( vm.bic )
    pe <- movMF(rad.comp, onp, nruns=nruns)

    ope.mu <- atan2(pe$theta[,2], pe$theta[,1])
    ope.kappa <- sqrt(rowSums(pe$theta^2))
    ope.alpha <- pe$alpha

    mu <- matrix(ope.mu, nrow=1, ncol=length(ope.mu))
    k <- matrix(ope.kappa, nrow=1, ncol=length(ope.kappa))
    alpha <- matrix(ope.alpha, nrow=1 , ncol=length(ope.alpha))


    theta<- matrix(seq(0, 2*pi, length.out = n.values))
    a1<- matrix(rep(1,length(mu)), ncol=length(mu), nrow=1)
    a2<- matrix(rep(1,length(theta)), nrow=length(theta), ncol=1)
    e.c<- exp((cos(theta%*%a1-a2%*%mu)*a2%*%k))
    f.theta <- e.c%*%t(alpha/(2*pi*I.0(k)))


    c.ang <- circular(rad(ang), units="radians", template="geographics")
    k.ang <- kern.den.circ(na.omit(c.ang), bw=bw.pi(na.omit(c.ang)), len=n.values)



    #dev.new()
    par(family="serif", mar = c(5,5,4,2))
    plot(k.ang, plot.type = "l", axes = FALSE, ylim = c(0, 1.1), lwd = 2,
         xlab = expression(bold("Direction"(N, Clockwise, degree))),
         ylab = expression(bold("Probability")),
         main = main, font.lab = 2, font.main = 2,
         cex.main = 2, cex.lab = cex.size)

    grid()
    box()

    axis(1, seq(0, 2*pi, pi/2), labels=seq(0, 360, 90), cex.axis=cex.size, font = 2)
    axis(2, las=2, cex.axis=cex.size, font = 2)
    lines(theta, f.theta, type="l", col=2, lwd=2)
    hist(rad(ang), prob = TRUE, nclass = h.class, add = TRUE)

    legend("topright", legend=c("Kernel", "MoV"), cex = cex.size,
           col=1:2, lty=1, lwd=3, horiz=TRUE, text.font = 2)

    mtext(paste("Optimal Number of Parameters =",onp), adj = 0.02,
          padj = 2.2, font = 2, cex = 1.1)



    mov.value <- list(Bic=vm.bic, n.OP=as.numeric(onp),
                      parameters = data.frame(Mu=ifelse(mu[1,]<0,(mu[1,]+2*pi)*180/pi, mu[1,]*180/pi),
                                              Kappa=k[1,], a=alpha[1,]),
                      values=list(MoV.radians=as.numeric(f.theta),
                                  Kernel.radians=c(rev(k.ang$y)[(0.75*n.values+1):n.values],
                                                   rev(k.ang$y)[1:(0.75*n.values)])))



    points(x = mov.value$parameters$Mu/(180/pi), y = f.theta[mov.value$parameters$Mu*(n.values/360), 1],
           cex=cex.size, pch=16)

    label <- onp %>% seq %>% paste0("$\\mu_{", ., "}$") %>% latex2exp
    text(x = mov.value$parameters$Mu/(180/pi), y = f.theta[mov.value$parameters$Mu*(n.values/360), 1],
         labels = label, pos = 3, offset = 1.3, cex = cex.size, font = 2)



    class(mov.value) <- "fmov"
    print.fmov <<- function(x) print(x[c(1:3)])
    return(mov.value)



  }