rm(list = ls()) # Load necessary libraries library(circular) library(polynom) # Function for retrieving mean of least scattering values4 around given angle ma = function(angle_h, data) { angle_d = sqrt((data -angle_h)^2) # calculate differences of data from given angle ord = order(angle_d, data) # order data vectors of differences and original angles angle_d2 = rbind(angle_d, data)[2,ord] # ordering original angles after differences angle = mean(angle_d2[1:10]) # calculate mean angle around the given input angle return(angle) # return this meanangle } #x.s.1<-0:15 #x.s<-sample(x.s.1,1000,replace=T) #x.s<-runif(100, 1, 180) x.s <- rnorm(100, 50, 15) x.s <- 1:180 #x.s <- c(rnorm(50, 40, 12), rnorm(50, 100, 15)) #x.s<-c(0,45,90,135,180) # prepare raw data x<-x.s x <- as.numeric(x[!is.na(x)]) x <- ifelse(x < 0, x + 360, x) # convert negative angles x <- ifelse(x > 360, x - 360, x) # convert angles > 2 Pi x <- ifelse(x > 180, x - 180, x) # convert angles > Pi data<-x # calculation of cumulative angle differences xd <- seq(1:(length(x)^2)) # create target vector x_2 <- rep(x[1:length(x)], each = length(x)) x_3 <- rep(x[1:length(x)], length(x)) xd <- ifelse(abs(x_2 - x_3) < abs(180 + x_2 - x_3), abs(x_2 - x_3), abs(x_3 - x_2)) # Create histogram from cumulative angle differences hist_xd <- hist(xd, plot = F, breaks = seq(0, 180)) y_xd <- (hist_xd$counts) y_xd[1:5] <- mean (y_xd[2:15]) x_xd <- seq(1, length(y_xd)) y_xd <- y_xd / max(y_xd) # polynomial model, seven components fit <- lm(y_xd ~ x_xd + I (x_xd^2) + I (x_xd^3) + I (x_xd^4) + I (x_xd^5) + I (x_xd^6) + I (x_xd^7)) fit_pol <- polynomial(fit$coefficients) fit_fct <- as.function(fit_pol) par(mfrow=c(1,2)) rose.diag(rad(c(x, x+180)), bins = 16, prop = 2) # Calculate and rescale extreme points x_ep <- as.numeric(solve(deriv(fit_pol), 0)) x_ep_a <- x_ep for (i in 1:6) { x_ep[i] <- ifelse(abs(x_ep[i] - x_ep[i + 1]) < .0001, NA, x_ep[i]) } x_ep[6] <- x_ep_a[6] x_ep <- x_ep[!is.na(x_ep)] y_ep <- fit_fct(x_ep) y_ep <- ifelse(y_ep < 0, 0, y_ep) x_ep[1] <- ifelse(x_ep[1] < 0, NA, x_ep[1]) y_ep[1] <- ifelse(x_ep[1] < 0, NA, x_ep[1]) x_ep <- x_ep[!is.na(x_ep)] y_ep <- y_ep[!is.na(y_ep)] x_ep5 <- rep(180, 5) y_ep5 <- rep(0, 5) x_ep5[1:length(x_ep)] <- x_ep y_ep5[1:length(y_ep)] <- y_ep x_ep <- x_ep5 y_ep <- y_ep5 # Sort extreme points by y-value u_ep <- rbind(y_ep, x_ep) ii <- order(-y_ep, y_ep, x_ep) s_ep <- rbind(y_ep, x_ep)[,ii] plot(x_xd, y_xd,main="modality") lines(x_xd, fit_fct(x_xd), col = 2) abline(v = x_ep, col = 5) # Calculate modality indicators # Iuf, Index of uniformity term1 <- sum(abs(y_xd - ((180 - x_xd) / 180))) / 90 Iuf <- (1 - term1) Iuf <- ifelse(Iuf < 0, 0, Iuf) # Ibm, Index of bimodality i <- 1:7 Ibm <- (abs(y_ep[i[x_ep == s_ep[2,2]]] - ((y_ep[i[x_ep == s_ep[2,2]]-1] + y_ep[i[x_ep == s_ep[2,2]]+1]) / 2))) Ibm <- ifelse(y_ep[i[x_ep == s_ep[2,2]]-1] - s_ep[1,1] < .001, 2 * abs(y_ep[i[x_ep == s_ep[2,2]]]), Ibm) Ibm <- as.vector(Ibm) Ibm <- 1.7 * Ibm[!is.na(Ibm)] - Iuf Ibm <- ifelse(Ibm < 0, 0, Ibm) m <- (y_ep[4] - y_ep[2]) / (x_ep[4] - x_ep[2]) n <- y_ep[2] - m * x_ep[2] x_ep2_ep4 <- x_xd[x_xd >= x_ep[2] & x_xd <= x_ep[4]] y_xd_lin <- m * x_ep2_ep4 + n y_diff <- y_xd[x_xd >= x_ep[2] & x_xd <= x_ep[4]] - y_xd_lin Ibm <- 2 * max(y_diff) # Ium, Index of unimodality Ium <- (term1 - Ibm) Ium <- ifelse(Ium < 0, 0, Ium) # Rescale resulting indices Isum <- sum(Ibm, Iuf, Ium) Ibm <- Ibm / Isum Iuf <- Iuf / Isum Ium <- Ium / Isum # Infering absolute angle maxima # Group data to retrieve most common angle and calculation of absolute value for this angle H1 = hist(x, breaks = c(seq(00, 180, 10)), plot = F) # grouping of original data for finding the mode ii = order(H1$counts, H1$mids, decreasing = T) # order data groups in descending order asort1 = rbind(H1$counts, H1$mids)[2,ii] # order mid group angles after highest counting group angle_h = asort1[1] # retrieve most common angle from ordered data vectors angle1 = ma(angle_h,x) # calling function to retrieve mean angle # Calculation of the absolute value for the second-most angle angle2a = ma(angle1 + x_ep[2],x) # calling function to retrieve mean angle for most common angle plus mean angle angle2b = ma(angle1 - x_ep[2],x) # calling function to retrieve mean angle for most common angle minus mean angle H2 = hist(c(angle2a, angle2b), breaks = c(seq(0, 180, 10)), plot = F) # group angles respctively to original data ii = order(H2$counts, H2$mids, decreasing = T) # order groups by frequency iia =ii[1] # retrieve most common group, i.e. group with angle2a iib = ii[2] # retrieve second most common group, i.e. group with angle2b angle2 = ifelse(H1$counts[iia] > H1$counts[iib], angle2a, angle2b) # decide, group with which angle is suitable div_a<-abs(angle1-angle2) div_a2<-abs(div_a-x_ep[2]) par(mfrow=c(1,1)) O <- matrix(nr = 3, nc = 2) O[,1] <- c("Iuf", "Ium", "Ibm") O[,2] <- round(c(Iuf, Ium, Ibm), 2) O