Create package structure and upload code

.Rbuildignore

@@ -0,0 +1,3 @@
+^.*\.Rproj$
+^\.Rproj\.user$
+^README\.Rmd$

.gitignore

@@ -0,0 +1,6 @@
+.DS_Store
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata
+inst/doc

DESCRIPTION

@@ -0,0 +1,50 @@
+Type: Package
+Package: meteorit
+Title: Mixtures-of-ExperTs modEling for cOmplex and non-noRmal dIsTributions
+Version: 0.1.0
+Authors@R: c(person("Faicel", "Chamroukhi", role = c("aut"),
+    email = "faicel.chamroukhi@unicaen.fr",
+    comment = c(ORCID = "0000-0002-5894-3103")),
+    person("Marius", "Bartcus", role = c("aut"), 
+    email = "marius.bartcus@gmail.com"),
+    person("Florian", "Lecocq", role = c("aut", "cre"),
+    email = "florian.lecocq@outlook.com"))
+Description: MEteoritS is a toolbox containg several original and flexible mixtures-of-experts models to model, 
+    cluster and classify heteregenous data in many complex situations where the data are distributed according non-normal, 
+    possibly skewed distributions, and when they might be corrupted by atypical observations. 
+    The toolbox contains in particular sparse mixture-of-experts models for high-dimensional data.
+License: GPL (>= 3)
+Depends: 
+    R (>= 2.10)
+Imports: 
+    methods,
+    stats,
+    Rcpp
+Suggests: 
+    knitr,
+    rmarkdown
+LinkingTo: 
+    Rcpp,
+    RcppArmadillo
+Collate:
+    meteorit-package.R
+    RcppExports.R
+    utils.R
+    FData.R
+    ParamSNMoE.R
+    ParamStMoE.R
+    ParamTMoE.R
+    StatSNMoE.R
+    StatStMoE.R
+    StatTMoE.R
+    ModelSNMoE.R
+    ModelStMoE.R
+    ModelTMoE.R
+    emSNMoE.R
+    emStMoE.R
+    emTMoE.R
+VignetteBuilder: knitr
+Encoding: UTF-8
+LazyData: true
+Roxygen: list(markdown = TRUE)
+RoxygenNote: 6.1.1

NAMESPACE

@@ -0,0 +1,19 @@
+# Generated by roxygen2: do not edit by hand
+
+export(emSNMoE)
+export(emStMoE)
+export(emTMoE)
+exportClasses(FData)
+exportClasses(ModelSNMoE)
+exportClasses(ModelStMoE)
+exportClasses(ModelTMoE)
+exportClasses(ParamSNMoE)
+exportClasses(ParamStMoE)
+exportClasses(ParamTMoE)
+exportClasses(StatSNMoE)
+exportClasses(StatStMoE)
+exportClasses(StatTMoE)
+import(methods)
+importFrom(Rcpp,sourceCpp)
+importFrom(pracma,fzero)
+useDynLib(meteorit, .registration = TRUE)

R/FData.R

@@ -0,0 +1,39 @@
+#' A Reference Class which represents functional data.
+#'
+#' FData is a reference class which represents general independent and
+#' identically distributed (i.i.d.) functional objects. The data can be
+#' ordered by time (functional time series). In the last case, the field `X`
+#' represents the time.
+#'
+#' @field X Numeric vector of length \emph{m}.
+#' @field Y Matrix of size \eqn{(n, m)} representing \emph{n}
+#' functions of `X` observed at points \eqn{1,\dots,m}.
+#' @export
+FData <- setRefClass(
+  "FData",
+  fields = list(
+    X = "numeric", # Covariates
+    Y = "matrix", # Response
+    m = "numeric",
+    n = "numeric",
+    vecY = "matrix"
+  ),
+  methods = list(
+
+    initialize = function(X = numeric(1), Y = matrix(1)) {
+
+      X <<- X
+      Y <<- as.matrix(Y)
+
+      n <<- nrow(Y)
+      m <<- ncol(Y)
+
+      vecY <<- matrix(t(Y), ncol = 1)
+
+      if (n == 1) {
+        Y <<- t(Y)
+      }
+
+    }
+  )
+)

R/ModelSNMoE.R

@@ -0,0 +1,76 @@
+#' A Reference Class which represents a fitted SNMoE model.
+#'
+#' ModelSNMoE represents a [SNMoE][ModelSNMoE] model for which parameters have
+#' been estimated.
+#'
+#' @usage NULL
+#' @field param A [ParamSNMoE][ParamSNMoE] object. It contains the estimated values of the parameters.
+#' @field stat A [StatSNMoE][StatSNMoE] object. It contains all the statistics associated to the SNMoE model.
+#' @seealso [ParamSNMoE], [StatSNMoE]
+#' @export
+ModelSNMoE <- setRefClass(
+  "ModelSNMoE",
+  fields = list(
+    param = "ParamSNMoE",
+    stat = "StatSNMoE"
+  ),
+  methods = list(
+    plot = function(what = c("meancurve", "confregions", "clusters", "loglikelihood"), ...) {
+
+      what <- match.arg(what, several.ok = TRUE)
+
+      oldpar <- par(no.readonly = TRUE)
+      on.exit(par(oldpar), add = TRUE)
+
+      colorsvec = rainbow(param$K)
+
+      if (any(what == "meancurve")) {
+        par(mfrow = c(2, 1), mai = c(0.6, 1, 0.5, 0.5), mgp = c(2, 1, 0))
+        plot.default(param$fData$X, param$fData$Y, ylab = "y", xlab = "x", cex = 0.7, pch = 3, ...)
+        title(main = "Estimated mean and experts")
+        for (k in 1:param$K) {
+          lines(param$fData$X, stat$Ey_k[, k], col = "red", lty = "dotted", lwd = 1.5, ...)
+        }
+        lines(param$fData$X, stat$Ey, col = "red", lwd = 1.5, ...)
+
+        plot.default(param$fData$X, stat$piik[, 1], type = "l", xlab = "x", ylab = "Mixing probabilities", col = colorsvec[1], ...)
+        title(main = "Mixing probabilities")
+        for (k in 2:param$K) {
+          lines(param$fData$X, stat$piik[, k], col = colorsvec[k], ...)
+        }
+      }
+
+      if (any(what == "confregions")) {
+        # Data, Estimated mean functions and 2*sigma confidence regions
+        par(mfrow = c(1, 1), mai = c(0.6, 1, 0.5, 0.5), mgp = c(2, 1, 0))
+        plot.default(param$fData$X, param$fData$Y, ylab = "y", xlab = "x", cex = 0.7, pch = 3, ...)
+        title(main = "Estimated mean and confidence regions")
+        lines(param$fData$X, stat$Ey, col = "red", lwd = 1.5, ...)
+        lines(param$fData$X, stat$Ey - 2 * sqrt(stat$Vary), col = "red", lty = "dotted", lwd = 1.5, ...)
+        lines(param$fData$X, stat$Ey + 2 * sqrt(stat$Vary), col = "red", lty = "dotted", lwd = 1.5, ...)
+      }
+
+      if (any(what == "clusters")) {
+        # Obtained partition
+        par(mfrow = c(1, 1), mai = c(0.6, 1, 0.5, 0.5), mgp = c(2, 1, 0))
+        plot.default(param$fData$X, param$fData$Y, ylab = "y", xlab = "x", cex = 0.7, pch = 3, ...)
+        title(main = "Estimated experts and clusters")
+        for (k in 1:param$K) {
+          lines(param$fData$X, stat$Ey_k[, k], col = colorsvec[k], lty = "dotted", lwd = 1.5, ...)
+        }
+        for (k in 1:param$K) {
+          index <- stat$klas == k
+          points(param$fData$X[index], param$fData$Y[index, ], col = colorsvec[k], cex = 0.7, pch = 3, ...)
+        }
+      }
+
+      if (any(what == "loglikelihood")) {
+        # Observed data log-likelihood
+        par(mfrow = c(1, 1), mai = c(0.6, 1, 0.5, 0.5), mgp = c(2, 1, 0))
+        plot.default(unlist(stat$stored_loglik), type = "l", col = "blue", xlab = "EM iteration number", ylab = "Observed data log-likelihood", ...)
+        title(main = "Log-Likelihood")
+      }
+
+    }
+  )
+)

R/ModelStMoE.R

@@ -0,0 +1,76 @@
+#' A Reference Class which represents a fitted StMoE model.
+#'
+#' ModelStMoE represents a [StMoE][ModelStMoE] model for which parameters have
+#' been estimated.
+#'
+#' @usage NULL
+#' @field param A [ParamStMoE][ParamStMoE] object. It contains the estimated values of the parameters.
+#' @field stat A [StatStMoE][StatStMoE] object. It contains all the statistics associated to the StMoE model.
+#' @seealso [ParamStMoE], [StatStMoE]
+#' @export
+ModelStMoE <- setRefClass(
+  "ModelStMoE",
+  fields = list(
+    param = "ParamStMoE",
+    stat = "StatStMoE"
+  ),
+  methods = list(
+    plot = function(what = c("meancurve", "confregions", "clusters", "loglikelihood"), ...) {
+
+      what <- match.arg(what, several.ok = TRUE)
+
+      oldpar <- par(no.readonly = TRUE)
+      on.exit(par(oldpar), add = TRUE)
+
+      colorsvec = rainbow(param$K)
+
+      if (any(what == "meancurve")) {
+        par(mfrow = c(2, 1), mai = c(0.6, 1, 0.5, 0.5), mgp = c(2, 1, 0))
+        plot.default(param$fData$X, param$fData$Y, ylab = "y", xlab = "x", cex = 0.7, pch = 3, ...)
+        title(main = "Estimated mean and experts")
+        for (k in 1:param$K) {
+          lines(param$fData$X, stat$Ey_k[, k], col = "red", lty = "dotted", lwd = 1.5, ...)
+        }
+        lines(param$fData$X, stat$Ey, col = "red", lwd = 1.5, ...)
+
+        plot.default(param$fData$X, stat$piik[, 1], type = "l", xlab = "x", ylab = "Mixing probabilities", col = colorsvec[1], ...)
+        title(main = "Mixing probabilities")
+        for (k in 2:param$K) {
+          lines(param$fData$X, stat$piik[, k], col = colorsvec[k], ...)
+        }
+      }
+
+      if (any(what == "confregions")) {
+        # Data, Estimated mean functions and 2*sigma confidence regions
+        par(mfrow = c(1, 1), mai = c(0.6, 1, 0.5, 0.5), mgp = c(2, 1, 0))
+        plot.default(param$fData$X, param$fData$Y, ylab = "y", xlab = "x", cex = 0.7, pch = 3, ...)
+        title(main = "Estimated mean and confidence regions")
+        lines(param$fData$X, stat$Ey, col = "red", lwd = 1.5)
+        lines(param$fData$X, stat$Ey - 2 * sqrt(stat$Vary), col = "red", lty = "dotted", lwd = 1.5, ...)
+        lines(param$fData$X, stat$Ey + 2 * sqrt(stat$Vary), col = "red", lty = "dotted", lwd = 1.5, ...)
+      }
+
+      if (any(what == "clusters")) {
+        # Obtained partition
+        par(mfrow = c(1, 1), mai = c(0.6, 1, 0.5, 0.5), mgp = c(2, 1, 0))
+        plot.default(param$fData$X, param$fData$Y, ylab = "y", xlab = "x", cex = 0.7, pch = 3, ...)
+        title(main = "Estimated experts and clusters")
+        for (k in 1:param$K) {
+          lines(param$fData$X, stat$Ey_k[, k], col = colorsvec[k], lty = "dotted", lwd = 1.5, ...)
+        }
+        for (k in 1:param$K) {
+          index <- stat$klas == k
+          points(param$fData$X[index], param$fData$Y[index, ], col = colorsvec[k], cex = 0.7, pch = 3, ...)
+        }
+      }
+
+      if (any(what == "loglikelihood")) {
+        # Observed data log-likelihood
+        par(mfrow = c(1, 1), mai = c(0.6, 1, 0.5, 0.5), mgp = c(2, 1, 0))
+        plot.default(unlist(stat$stored_loglik), type = "l", col = "blue", xlab = "EM iteration number", ylab = "Observed data log-likelihood", ...)
+        title(main = "Log-Likelihood")
+      }
+
+    }
+  )
+)

R/ModelTMoE.R

@@ -0,0 +1,76 @@
+#' A Reference Class which represents a fitted TMoE model.
+#'
+#' ModelMRHLP represents a [TMoE][ModelTMoE] model for which parameters have
+#' been estimated.
+#'
+#' @usage NULL
+#' @field param A [ParamTMoE][ParamTMoE] object. It contains the estimated values of the parameters.
+#' @field stat A [StatTMoE][StatTMoE] object. It contains all the statistics associated to the TMoE model.
+#' @seealso [ParamTMoE], [StatTMoE]
+#' @export
+ModelTMoE <- setRefClass(
+  "ModelTMoE",
+  fields = list(
+    param = "ParamTMoE",
+    stat = "StatTMoE"
+  ),
+  methods = list(
+    plot = function(what = c("meancurve", "confregions", "clusters", "loglikelihood"), ...) {
+
+      what <- match.arg(what, several.ok = TRUE)
+
+      oldpar <- par(no.readonly = TRUE)
+      on.exit(par(oldpar), add = TRUE)
+
+      colorsvec = rainbow(param$K)
+
+      if (any(what == "meancurve")) {
+        par(mfrow = c(2, 1), mai = c(0.6, 1, 0.5, 0.5), mgp = c(2, 1, 0))
+        plot.default(param$fData$X, param$fData$Y, ylab = "y", xlab = "x", cex = 0.7, pch = 3, ...)
+        title(main = "Estimated mean and experts")
+        for (k in 1:param$K) {
+          lines(param$fData$X, stat$Ey_k[, k], col = "red", lty = "dotted", lwd = 1.5, ...)
+        }
+        lines(param$fData$X, stat$Ey, col = "red", lwd = 1.5, ...)
+
+        plot.default(param$fData$X, stat$piik[, 1], type = "l", xlab = "x", ylab = "Mixing probabilities", col = colorsvec[1], ...)
+        title(main = "Mixing probabilities")
+        for (k in 2:param$K) {
+          lines(param$fData$X, stat$piik[, k], col = colorsvec[k], ...)
+        }
+      }
+
+      if (any(what == "confregions")) {
+        # Data, Estimated mean functions and 2*sigma confidence regions
+        par(mfrow = c(1, 1), mai = c(0.6, 1, 0.5, 0.5), mgp = c(2, 1, 0))
+        plot.default(param$fData$X, param$fData$Y, ylab = "y", xlab = "x", cex = 0.7, pch = 3, ...)
+        title(main = "Estimated mean and confidence regions")
+        lines(param$fData$X, stat$Ey, col = "red", lwd = 1.5)
+        lines(param$fData$X, stat$Ey - 2 * sqrt(stat$Vary), col = "red", lty = "dotted", lwd = 1.5, ...)
+        lines(param$fData$X, stat$Ey + 2 * sqrt(stat$Vary), col = "red", lty = "dotted", lwd = 1.5, ...)
+      }
+
+      if (any(what == "clusters")) {
+        # Obtained partition
+        par(mfrow = c(1, 1), mai = c(0.6, 1, 0.5, 0.5), mgp = c(2, 1, 0))
+        plot.default(param$fData$X, param$fData$Y, ylab = "y", xlab = "x", cex = 0.7, pch = 3, ...)
+        title(main = "Estimated experts and clusters")
+        for (k in 1:param$K) {
+          lines(param$fData$X, stat$Ey_k[, k], col = colorsvec[k], lty = "dotted", lwd = 1.5, ...)
+        }
+        for (k in 1:param$K) {
+          index <- stat$klas == k
+          points(param$fData$X[index], param$fData$Y[index, ], col = colorsvec[k], cex = 0.7, pch = 3, ...)
+        }
+      }
+
+      if (any(what == "loglikelihood")) {
+        # Observed data log-likelihood
+        par(mfrow = c(1, 1), mai = c(0.6, 1, 0.5, 0.5), mgp = c(2, 1, 0))
+        plot.default(unlist(stat$stored_loglik), type = "l", col = "blue", xlab = "EM iteration number", ylab = "Observed data log-likelihood", ...)
+        title(main = "Log-Likelihood")
+      }
+
+    }
+  )
+)