Von Mises Combination Method of Moments

Von_Mises_Combination_Method_Of_Moments.R

library(circular)

## Helper function to generate a circular object in degrees, 0-360
circularDeg <- function(x) circular(x,template = "geographics",units = "degrees", modulo = "2pi")
theta <- circularDeg(seq(from = 0, to = 360,by = 1))

m1 <- circularDeg(180)
k1 <- 2.1
w1 <- 65
dvm1 <- dvonmises(theta, mu = m1, kappa = k1)

m2 <- circularDeg(175)
k2 <- 1.4
w2 <- 55
dvm2 <- dvonmises(theta, mu = m2, kappa = k2)

m3 <- circularDeg(205)
k3 <- 3.1
w3 <- 17
dvm3 <- dvonmises(theta, mu = m3, kappa = k3)


mint1 <- (Arg(circular::I.1(k1)/circular::I.0(k1) * exp(1i * toRad(m1))) %% (2*pi)) * 180 / pi * w1
mint2 <- (Arg(circular::I.1(k2)/circular::I.0(k2) * exp(1i * toRad(m2))) %% (2*pi)) * 180 / pi * w2
mint3 <- (Arg(circular::I.1(k3)/circular::I.0(k3) * exp(1i * toRad(m3))) %% (2*pi)) * 180 / pi * w3

(sum(c(mint1,mint2,mint3))/sum(w1,w2,w3))


makeJoint2 <- function(dList = NULL) {
  
  for(i in seq_len(length(dList))){
    m <- as.numeric(dList[[i]]$mu) * pi / 180
    k <- dList[[i]]$kappa
    w <- dList[[i]]$weight
    if(i == 1L){
      M <- m
      K <- k
      W <- w
      print(paste0("Sample #",i," Distribution: m = ",M / pi * 180,", k = ",K))

    } else {
      
      M <- atan2(sin(M)*W + sin(m)*w, cos(M)*W + cos(m)*w) %% (2*pi)      
      R2 <- (1/(W+w) * sum(c(rep(cos(M),W), rep(cos(m),w))))^2 + (1/(W+w) * sum(c(rep(sin(M),W), rep(sin(m),w))))^2
      R <- R2^0.5
      eq <- function(k) circular::I.1(k)/circular::I.0(k) - R + 1e-3
      
      K <- uniroot(eq,c(0,600))$root
      
      ## Kappa Bias Correction for small samples
      if (K < 2) {
        K = max(K - 2 * (i * K) ^-1, 0);
      } else {
        K = ((i - 1)^3 * K) / (i^3 + i);
      }
      
      W <- W + w
      print(paste0("Sample #",i," Distribution: m = ",M / pi * 180,", R2 = ",R2,", R = ",R2,", k = ",K))
    }
    
    
  }
  return(list(mu = circular(M / pi * 180,template = "geographics",units = "degrees", modulo = "2pi"),
              kappa = K,
              weight = W))
  
}

dList <- list(list(mu = m1,kappa = k1, weight = w1),
              list(mu = m2,kappa = k2, weight = w2),
              list(mu = m3,kappa = k3, weight = w3))
              
Estimate <- makeJoint2(dList)

# [1] "Sample #1 Distribution: m = 180, k = 2.1"
# [1] "Sample #2 Distribution: m = 177.708464726633, R2 = 0.999445331142574, R = 0.999445331142574, k = 39.1678866133919"
# [1] "Sample #3 Distribution: m = 181.010389041222, R2 = 0.981222564268574, R = 0.981222564268574, k = 12.8474416215592"