Added functional learn_simple_partition

R/admix1-learn_simple_partition.R

@@ -0,0 +1,60 @@
+#' @title Learn simple partition
+#' @param ordered_points ordered_points
+#' @param fii fii
+#' @description Goes through all the ways to divide the points into two or 
+#' three groups. Chooses the partition which obtains highest logml.
+#' @export
+learn_simple_partition <- function(ordered_points, fii) {
+    npoints <- length(ordered_points)
+    
+    # One cluster:
+    val <- calculatePopLogml(ordered_points, fii)
+    bestValue <- val
+    best_type <- 'single'
+    
+    # Two clusters:
+    for (i in 1:(npoints - 1)) {
+        # The right endpoint of the first cluster.
+        val_1 <- calculatePopLogml(ordered_points[1:i], fii)
+        val_2 <- calculatePopLogml(ordered_points[(i + 1):length(ordered_points)], fii)
+        total <- val_1 + val_2
+        if (total > bestValue) {
+            bestValue <- total
+            best_type <- 'double'
+            best_i <- i
+        }
+    }
+
+    # Three clusters:
+    for (i in 1:(npoints - 2)) {
+        for (j in (i + 1):(npoints - 1)) {
+            val_1 <- calculatePopLogml(ordered_points[1:i], fii)
+            val_2 <- calculatePopLogml(ordered_points[(i + 1):j], fii)
+            val_3 <- calculatePopLogml(ordered_points[(j + 1):length(ordered_points)], fii)
+            total <- val_1 + val_2 + val_3
+            if (total > bestValue) {
+                bestValue <- total
+                best_type <- 'triple'
+                best_i <- i
+                best_j <- j
+            }
+        }
+    }
+    
+    part = matrix(0, npoints, 1)
+
+    switch(best_type,
+        'single' = {
+            part <- matrix(1, npoints, 1)
+        },
+        'double' = {
+            part[1:best_i] <- 1
+            part[(best_i + 1):length(part)] <- 2
+        },
+        'triple' = {
+            part[1:best_i] <- 1
+            part[(best_i + 1):best_j] <- 2
+            part[(best_j + 1):length(part)] <- 3
+        })
+    return(part)
+}

R/calculatePopLogml.R

@@ -0,0 +1,19 @@
+#' @title Calculate log marginal likelihood
+#' @description Calculates fuzzy (log) marginal likelihood for a population of 
+#' real values using estimate "fii" for the dispersion value, and Jeffreys prior
+#' for the mean parameter.
+#' @param points points
+#' @param fii fii
+calculatePopLogml <- function(points, fii) {
+    n <- length(points)
+    fuzzy_ones <- sum(points)
+    fuzzy_zeros <- n - fuzzy_ones
+    val = log(gamma(1)) -
+        log(gamma(1 + n / fii)) +
+        log(gamma(0.5 + fuzzy_ones / fii)) +
+        log(gamma(0.5 + fuzzy_zeros / fii)) -
+        log(gamma(0.5)) -
+        log(gamma(0.5))
+    return(val)
+}
+