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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.