diff --git a/.Rbuildignore b/.Rbuildignore
new file mode 100644
index 0000000..931662b
--- /dev/null
+++ b/.Rbuildignore
@@ -0,0 +1,3 @@
+^.*\.Rproj$
+^\.Rproj\.user$
+^README\.Rmd$
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..53b5c81
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,6 @@
+.DS_Store
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata
+inst/doc
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..e7d949f
--- /dev/null
+++ b/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
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..54722a6
--- /dev/null
+++ b/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)
diff --git a/R/FData.R b/R/FData.R
new file mode 100644
index 0000000..0c46ebf
--- /dev/null
+++ b/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)
+      }
+
+    }
+  )
+)
diff --git a/R/ModelSNMoE.R b/R/ModelSNMoE.R
new file mode 100644
index 0000000..4314940
--- /dev/null
+++ b/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")
+      }
+
+    }
+  )
+)
diff --git a/R/ModelStMoE.R b/R/ModelStMoE.R
new file mode 100644
index 0000000..116c249
--- /dev/null
+++ b/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")
+      }
+
+    }
+  )
+)
diff --git a/R/ModelTMoE.R b/R/ModelTMoE.R
new file mode 100644
index 0000000..4934fe5
--- /dev/null
+++ b/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")
+      }
+
+    }
+  )
+)
diff --git a/R/ParamSNMoE.R b/R/ParamSNMoE.R
new file mode 100755
index 0000000..3083180
--- /dev/null
+++ b/R/ParamSNMoE.R
@@ -0,0 +1,147 @@
+#' A Reference Class which contains parameters of a SNMoE model.
+#'
+#' ParamSNMoE contains all the parameters of a SNMoE model.
+#'
+#' @field fData [FData][FData] object representing the sample.
+#' @field K The number of mixture components.
+#' @field p The order of the polynomial regression.
+#' @field q The dimension of the logistic regression. For the purpose of
+#' segmentation, it must be set to 1.
+#' @field nu degree of freedom
+#' @field alpha is the parameter vector of the logistic model with \eqn{alpha_K} being the null vector.
+#' @field beta is the vector of regression coefficients of component k,
+#' the updates for each of the expert component parameters consist in analytically solving a weighted
+#' Gaussian linear regression problem.
+#' @field sigma The variances for the \emph{K} mixture component.
+#' @field lambda skewness parameter
+#' @field delta the skewness parameter lambda (by equivalence delta)
+#' @seealso [FData]
+#' @export
+ParamSNMoE <- setRefClass(
+  "ParamSNMoE",
+  fields = list(
+    fData = "FData",
+    phiBeta = "list",
+    phiAlpha = "list",
+
+    K = "numeric",
+    # number of regimes
+    p = "numeric",
+    # dimension of beta (order of polynomial regression)
+    q = "numeric",
+    # dimension of w (order of logistic regression)
+    nu = "numeric", # degree of freedom
+
+    alpha = "matrix",
+    beta = "matrix",
+    sigma = "matrix",
+    lambda = "matrix",
+    delta = "matrix"
+  ),
+  methods = list(
+    initialize = function(fData = FData(numeric(1), matrix(1)), K = 1, p = 3, q = 1) {
+
+      fData <<- fData
+
+      phiBeta <<- designmatrix(x = fData$X, p = p)
+      phiAlpha <<- designmatrix(x = fData$X, p = q)
+
+      nu <<- (p + q + 3) * K - (q + 1)
+      K <<- K
+      p <<- p
+      q <<- q
+
+      alpha <<- matrix(0, q + 1, K - 1)
+      beta <<- matrix(NA, p + 1, K)
+      sigma <<- matrix(NA, 1, K)
+      lambda <<- matrix(NA, K)
+      delta <<- matrix(NA, K)
+    },
+
+    initParam = function(try_EM, segmental = FALSE) {
+      alpha <<- matrix(runif((q + 1) * (K - 1)), nrow = q + 1, ncol = K - 1) #initialisation al??atoire du vercteur param???tre du IRLS
+
+      #Initialise the regression parameters (coeffecients and variances):
+      if (segmental == FALSE) {
+        Zik <- zeros(fData$n, K)
+
+        klas <- floor(K * matrix(runif(fData$n), fData$n)) + 1
+
+        Zik[klas %*% ones(1, K) == ones(fData$n, 1) %*% seq(K)] <- 1
+
+        Tauik <- Zik
+
+
+        #beta <<- matrix(0, modelRHLP$p + 1, modelRHLP$K)
+        #sigma <<- matrix(0, modelRHLP$K)
+
+        for (k in 1:K) {
+          Xk <- phiBeta$XBeta * (sqrt(Tauik[, k] %*% ones(1, p + 1)))
+          yk <- fData$Y * sqrt(Tauik[, k])
+
+          beta[, k] <<- solve(t(Xk) %*% Xk) %*% t(Xk) %*% yk
+
+          sigma[k] <<- sum(Tauik[, k] * ((fData$Y - phiBeta$XBeta %*% beta[, k]) ^ 2)) / sum(Tauik[, k])
+        }
+      }
+      else{
+        #segmental : segment uniformly the data and estimate the parameters
+        nk <- round(fData$n / K) - 1
+
+        for (k in 1:K) {
+          i <- (k - 1) * nk + 1
+          j <- (k * nk)
+          yk <- matrix(fData$Y[i:j])
+          Xk <- phiBeta$XBeta[i:j,]
+
+          beta[, k] <<- solve(t(Xk) %*% Xk) %*% (t(Xk) %*% yk)
+
+          muk <- Xk %*% beta[, k]
+
+          sigma[k] <<- t(yk - muk) %*% (yk - muk) / length(yk)
+        }
+      }
+
+      if (try_EM == 1) {
+        alpha <<- zeros(q + 1, K - 1)
+      }
+
+      # Initialize the skewness parameter Lambdak (by equivalence delta)
+      delta <<- -1 + 2 * rand(1, K)
+
+      lambda <<- delta / sqrt(1 - delta ^ 2)
+    },
+
+    MStep = function(statSNMoE, verbose_IRLS) {
+      # M-Step
+      res_irls <- IRLS(phiAlpha$XBeta, statSNMoE$tik, ones(nrow(statSNMoE$tik), 1), alpha, verbose_IRLS)
+      statSNMoE$piik <- res_irls$piik
+      reg_irls <- res_irls$reg_irls
+
+      alpha <<- res_irls$W
+
+      for (k in 1:K) {
+        #update the regression coefficients
+
+        tauik_Xbeta <- (statSNMoE$tik[, k] %*% ones(1, p + 1)) * phiBeta$XBeta
+        beta[, k] <<- solve((t(tauik_Xbeta) %*% phiAlpha$XBeta)) %*% (t(tauik_Xbeta) %*% (fData$Y - delta[k] * statSNMoE$E1ik[, k]))
+
+        # update the variances sigma2k
+
+        sigma[k] <<- sum(statSNMoE$tik[, k] * ((fData$Y - phiBeta$XBeta %*% beta[, k]) ^ 2 - 2 * delta[k] * statSNMoE$E1ik[, k] * (fData$Y - phiBeta$XBeta %*% beta[, k]) + statSNMoE$E2ik[, k])) / (2 * (1 - delta[k] ^ 2) * sum(statSNMoE$tik[, k]))
+
+        # update the deltak (the skewness parameter)
+        delta[k] <<- uniroot(f <- function(dlt) {
+          sigma[k] * dlt * (1 - dlt ^ 2) * sum(statSNMoE$tik[, k]) + (1 + dlt ^ 2) * sum(statSNMoE$tik[, k] * (fData$Y - phiBeta$XBeta %*% beta[, k]) * statSNMoE$E1ik[, k])
+          - dlt * sum(statSNMoE$tik[, k] * (statSNMoE$E2ik[, k] + (fData$Y - phiBeta$XBeta %*% beta[, k]) ^ 2))
+        }, c(-1, 1))$root
+
+
+        lambda[k] <<- delta[k] / sqrt(1 - delta[k] ^ 2)
+
+      }
+
+      return(reg_irls)
+    }
+  )
+)
diff --git a/R/ParamStMoE.R b/R/ParamStMoE.R
new file mode 100755
index 0000000..29b886f
--- /dev/null
+++ b/R/ParamStMoE.R
@@ -0,0 +1,165 @@
+#' A Reference Class which contains parameters of a MRHLP model.
+#'
+#' ParamMRHLP contains all the parameters of a MRHLP model.
+#'
+#' @field fData [FData][FData] object representing the sample.
+#' @field K The number of mixture components.
+#' @field p The order of the polynomial regression.
+#' @field q The dimension of the logistic regression. For the purpose of
+#' segmentation, it must be set to 1.
+#' @field nu degree of freedom
+#' @field alpha is the parameter vector of the logistic model with \eqn{alpha_K} being the null vector.
+#' @field beta is the vector of regression coefficients of component k,
+#' the updates for each of the expert component parameters consist in analytically solving a weighted
+#' Gaussian linear regression problem.
+#' @field sigma The variances for the \emph{K} mixture component.
+#' @field lambda skewness parameter
+#' @field delta the skewness parameter lambda (by equivalence delta)
+#' @field nuk degrees of freedom
+#' @seealso [FData]
+#' @export
+ParamStMoE <- setRefClass(
+  "ParamStMoE",
+  fields = list(
+    fData = "FData",
+    phiBeta = "list",
+    phiAlpha = "list",
+
+    K = "numeric",
+    # number of components
+    p = "numeric",
+    # dimension of beta (order of polynomial regression)
+    q = "numeric",
+    # dimension of w (order of logistic regression)
+    nu = "numeric", # degree of freedom
+
+    alpha = "matrix",
+    beta = "matrix",
+    sigma = "matrix",
+    lambda = "matrix",
+    delta = "matrix",
+    nuk = "matrix"
+  ),
+  methods = list(
+    initialize = function(fData = FData(numeric(1), matrix(1)), K = 1, p = 3, q = 1) {
+
+      fData <<- fData
+
+      phiBeta <<- designmatrix(x = fData$X, p = p)
+      phiAlpha <<- designmatrix(x = fData$X, p = q)
+
+      nu <<- (p + q + 3) * K - (q + 1)
+      K <<- K
+      p <<- p
+      q <<- q
+
+      alpha <<- matrix(0, q + 1, K - 1)
+      beta <<- matrix(NA, p + 1, K)
+      sigma <<- matrix(NA, 1, K)
+      lambda <<- matrix(NA, K)
+      delta <<- matrix(NA, K)
+    },
+
+    initParam = function(try_EM, segmental = FALSE) {
+      "Method to initialize parameters \\code{alpha}, \\code{beta} and
+      \\code{sigma}."
+
+      alpha <<- matrix(runif((q + 1) * (K - 1)), nrow = q + 1, ncol = K - 1) #initialisation al??atoire du vercteur param???tre du IRLS
+
+      #Initialise the regression parameters (coeffecients and variances):
+      if (segmental == FALSE) {
+        Zik <- zeros(fData$n, K)
+
+        klas <- floor(K * matrix(runif(fData$n), fData$n)) + 1
+
+        Zik[klas %*% ones(1, K) == ones(fData$n, 1) %*% seq(K)] <- 1
+
+        Tauik <- Zik
+
+
+        #beta <<- matrix(0, modelRHLP$p + 1, modelRHLP$K)
+        #sigma <<- matrix(0, modelRHLP$K)
+
+        for (k in 1:K) {
+          Xk <- phiBeta$XBeta * (sqrt(Tauik[, k] %*% ones(1, p + 1)))
+          yk <- fData$Y * sqrt(Tauik[, k])
+
+          beta[, k] <<- solve(t(Xk) %*% Xk) %*% t(Xk) %*% yk
+
+          sigma[k] <<- sum(Tauik[, k] * ((fData$Y - phiBeta$XBeta %*% beta[, k]) ^ 2)) / sum(Tauik[, k])
+        }
+      }
+      else{
+        #segmental : segment uniformly the data and estimate the parameters
+        nk <- round(fData$n / K) - 1
+
+        for (k in 1:K) {
+          i <- (k - 1) * nk + 1
+          j <- (k * nk)
+          yk <- matrix(fData$Y[i:j])
+          Xk <- phiBeta$XBeta[i:j,]
+
+          beta[, k] <<- solve(t(Xk) %*% Xk) %*% (t(Xk) %*% yk)
+
+          muk <- Xk %*% beta[, k]
+
+          sigma[k] <<- t(yk - muk) %*% (yk - muk) / length(yk)
+        }
+      }
+
+      if (try_EM == 1) {
+        alpha <<- rand(q + 1, K - 1)
+      }
+
+      # Initialize the skewness parameter Lambdak (by equivalence delta)
+      delta <<- -0.9 + 1.8 * rand(1, K)
+
+      lambda <<- delta / sqrt(1 - delta ^ 2)
+
+      # Intitialization of the degrees of freedm
+      nuk <<- 1 + 5 * rand(1, K)
+    },
+
+    MStep = function(statStMoE, verbose_IRLS) {
+      "Method used in the EM algorithm to learn the parameters of the StMoE model
+      based on statistics provided by \\code{statStMoE}."
+      # M-Step
+      res_irls <- IRLS(phiAlpha$XBeta, statStMoE$tik, ones(nrow(statStMoE$tik), 1), alpha, verbose_IRLS)
+      # statStMoE$piik <- res_irls$piik
+      reg_irls <- res_irls$reg_irls
+
+      alpha <<- res_irls$W
+
+      for (k in 1:K) {
+        #update the regression coefficients
+        TauikWik <- (statStMoE$tik[,k] * statStMoE$wik[,k]) %*% ones(1, p+1)
+        TauikX <- phiBeta$XBeta * (statStMoE$tik[,k] %*% ones(1, p+1))
+        betak <- solve((t(TauikWik * phiBeta$XBeta) %*% phiAlpha$XBeta)) %*% (t(TauikX) %*% ( (statStMoE$wik[,k] * fData$Y) - (delta[k] * statStMoE$E1ik[ ,k]) ))
+
+        beta[,k] <<- betak;
+        # update the variances sigma2k
+
+        sigma[k] <<- sum(statStMoE$tik[, k]*(statStMoE$wik[,k] * ((fData$Y-phiBeta$XBeta%*%betak)^2) - 2 * delta[k] * statStMoE$E1ik[,k] * (fData$Y - phiBeta$XBeta %*% betak) + statStMoE$E2ik[,k]))/(2*(1-delta[k]^2) * sum(statStMoE$tik[,k]))
+
+        sigmak <- sqrt(sigma[k])
+
+        # update the deltak (the skewness parameter)
+        delta[k] <<- uniroot(f <- function(dlt) {
+          return(dlt*(1-dlt^2)*sum(statStMoE$tik[, k])
+          + (1+ dlt^2)*sum(statStMoE$tik[, k] * statStMoE$dik[,k]*statStMoE$E1ik[,k]/sigmak)
+          - dlt * sum(statStMoE$tik[, k] * (statStMoE$wik[,k] * (statStMoE$dik[,k]^2) + statStMoE$E2ik[,k]/(sigmak^2))))
+        }, c(-1, 1))$root
+
+
+        lambda[k] <<- delta[k] / sqrt(1 - delta[k] ^ 2)
+
+
+        nuk[k] <<- uniroot(f <- function(nnu) {
+          return(- psigamma((nnu)/2) + log((nnu)/2) + 1 + sum(statStMoE$tik[,k] * (statStMoE$E3ik[,k] - statStMoE$wik[,k]))/sum(statStMoE$tik[,k]))
+        }, c(0.1, 200))$root
+      }
+
+      return(reg_irls)
+    }
+  )
+)
diff --git a/R/ParamTMoE.R b/R/ParamTMoE.R
new file mode 100755
index 0000000..cf3fb4a
--- /dev/null
+++ b/R/ParamTMoE.R
@@ -0,0 +1,149 @@
+#' A Reference Class which contains parameters of a TMoE model.
+#'
+#' ParamTMoE contains all the parameters of a TMoE model.
+#'
+#' @field fData [FData][FData] object representing the sample.
+#' @field K The number of mixture components.
+#' @field p The order of the polynomial regression.
+#' @field q The dimension of the logistic regression. For the purpose of
+#' segmentation, it must be set to 1.
+#' @field nu degree of freedom
+#' @field alpha is the parameter vector of the logistic model with \eqn{alpha_K} being the null vector.
+#' @field beta is the vector of regression coefficients of component k,
+#' the updates for each of the expert component parameters consist in analytically solving a weighted
+#' Gaussian linear regression problem.
+#' @field sigma The variances for the \emph{K} mixture component.
+#' @field delta the skewness parameter lambda (by equivalence delta)
+#' @seealso [FData]
+#' @importFrom pracma fzero
+#' @export
+ParamTMoE <- setRefClass(
+  "ParamTMoE",
+  fields = list(
+    fData = "FData",
+    phiBeta = "list",
+    phiAlpha = "list",
+
+    K = "numeric",
+    # number of regimes
+    p = "numeric",
+    # dimension of beta (order of polynomial regression)
+    q = "numeric",
+    # dimension of w (order of logistic regression)
+    nu = "numeric", # degree of freedom
+
+    alpha = "matrix",
+    beta = "matrix",
+    sigma = "matrix",
+    delta = "matrix"
+  ),
+  methods = list(
+    initialize = function(fData = FData(numeric(1), matrix(1)), K = 1, p = 3, q = 1) {
+
+      fData <<- fData
+
+      phiBeta <<- designmatrix(x = fData$X, p = p)
+      phiAlpha <<- designmatrix(x = fData$X, p = q)
+
+      K <<- K
+      p <<- p
+      q <<- q
+
+      nu <<- (p + q + 3) * K - (q + 1)
+
+      alpha <<- matrix(0, q + 1, K - 1)
+      beta <<- matrix(NA, p + 1, K)
+      sigma <<- matrix(NA, 1, K)
+      delta <<- matrix(NA, K)
+    },
+
+    initParam = function(try_EM, segmental = FALSE) {
+      alpha <<- matrix(runif((q + 1) * (K - 1)), nrow = q + 1, ncol = K - 1) #initialisation aleatoire du vercteur parametre du IRLS
+
+      #Initialise the regression parameters (coeffecients and variances):
+      if (segmental == FALSE) {
+        Zik <- zeros(n, K)
+
+        klas <- floor(K * matrix(runif(fData$n), fData$n)) + 1
+
+        Zik[klas %*% ones(1, K) == ones(fData$n, 1) %*% seq(K)] <- 1
+
+        Tauik <- Zik
+
+        #beta <<- matrix(0, modelRHLP$p + 1, modelRHLP$K)
+        #sigma <<- matrix(0, modelRHLP$K)
+
+        for (k in 1:K) {
+          Xk <- phiBeta$XBeta * (sqrt(Tauik[, k] %*% ones(1, p + 1)))
+          yk <- fData$Y * sqrt(Tauik[, k])
+
+          beta[, k] <<- solve(t(Xk) %*% Xk) %*% t(Xk) %*% yk
+
+          sigma[k] <<- sum(Tauik[, k] * ((fData$Y - phiBeta$XBeta %*% beta[, k]) ^ 2)) / sum(Tauik[, k])
+        }
+      }
+      else{
+        #segmental : segment uniformly the data and estimate the parameters
+        nk <- round(fData$n / K) - 1
+
+        for (k in 1:K) {
+          i <- (k - 1) * nk + 1
+          j <- (k * nk)
+          yk <- matrix(fData$Y[i:j])
+          Xk <- phiBeta$XBeta[i:j,]
+
+          beta[, k] <<- solve(t(Xk) %*% Xk) %*% (t(Xk) %*% yk)
+
+          muk <- Xk %*% beta[, k]
+
+          sigma[k] <<- t(yk - muk) %*% (yk - muk) / length(yk)
+        }
+      }
+
+      if (try_EM == 1) {
+        alpha <<- zeros(q + 1, K - 1)
+      }
+
+      # Initialize the skewness parameter Lambdak (by equivalence delta)
+      delta <<- 50 * rand(1, K)
+
+    },
+
+    MStep = function(statTMoE, verbose_IRLS) {
+      # M-Step
+      res_irls <- IRLS(phiAlpha$XBeta, statTMoE$tik, ones(nrow(statTMoE$tik), 1), alpha, verbose_IRLS)
+      statTMoE$piik <- res_irls$piik
+      reg_irls <- res_irls$reg_irls
+
+      alpha <<- res_irls$W
+
+      for (k in 1:K) {
+        #update the regression coefficients
+
+        Xbeta <- phiBeta$XBeta * (matrix( sqrt(statTMoE$tik[,k] * statTMoE$Wik[,k] )) %*% ones(1, p + 1))
+        yk <- fData$Y * sqrt(statTMoE$tik[,k] * statTMoE$Wik[,k])
+
+        #update the regression coefficients
+        beta[, k] <<- solve((t(Xbeta) %*% Xbeta)) %*% (t(Xbeta) %*% yk)
+
+        # update the variances sigma2k
+        sigma[k] <<- sum(statTMoE$tik[, k] * statTMoE$Wik[,k] * ((fData$Y - phiBeta$XBeta %*% beta[, k])^2)) / sum(statTMoE$tik[,k])
+
+        # if ECM (use an additional E-Step with the updatated betak and sigma2k
+        dik <- (fData$Y - phiBeta$XBeta %*% beta[, k]) / sqrt(sigma[k])
+
+
+        # update the deltak (the skewness parameter)
+
+
+        delta[k] <<- fzero(f <- function(dlt) {
+          return(-psigamma(dlt/2) + log(dlt/2) + 1 + (1 / sum(statTMoE$tik[, k])) * sum(statTMoE$tik[, k] * (log(statTMoE$Wik[,k]) - statTMoE$Wik[,k]))
+                 + psigamma((delta[k] + 1)/2) - log((delta[k] + 1)/2))
+        }, delta[k])$x
+
+      }
+
+      return(reg_irls)
+    }
+  )
+)
diff --git a/R/RcppExports.R b/R/RcppExports.R
new file mode 100644
index 0000000..b5a8233
--- /dev/null
+++ b/R/RcppExports.R
@@ -0,0 +1,11 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+IRLS <- function(X, Tau, Gamma, Winit, verbose = FALSE) {
+    .Call(`_meteorit_IRLS`, X, Tau, Gamma, Winit, verbose)
+}
+
+multinomialLogit <- function(W, X, Y, Gamma) {
+    .Call(`_meteorit_multinomialLogit`, W, X, Y, Gamma)
+}
+
diff --git a/R/StatSNMoE.R b/R/StatSNMoE.R
new file mode 100755
index 0000000..fda2100
--- /dev/null
+++ b/R/StatSNMoE.R
@@ -0,0 +1,193 @@
+#' A Reference Class which contains statistics of a SNMoE model.
+#'
+#' StatMRHLP contains all the parameters of a [SNMoE][ParamSNMoE] model.
+#'
+#' @field piik Matrix of size \eqn{(n, K)} representing the probabilities
+#' \eqn{P(zi = k; W) = P(z_{ik} = 1; W)}{P(zi = k; W) = P(z_ik = 1; W)} of the
+#' latent variable \eqn{zi,\ i = 1,\dots,m}{zi, i = 1,\dots,n}.
+#' @field z_ik Hard segmentation logical matrix of dimension \eqn{(n, K)}
+#' obtained by the Maximum a posteriori (MAP) rule:
+#' \eqn{z_{ik} = 1 \ \textrm{if} \ z_{ik} = \textrm{arg} \ \textrm{max}_{k} \
+#' P(z_i = k | Y, W, \beta);\ 0 \ \textrm{otherwise}}{z_ik = 1 if z_ik =
+#' arg max_k P(z_i = k | Y, W, \beta); 0 otherwise}, \eqn{k = 1,\dots,K}.
+#' @field klas Column matrix of the labels issued from `z_ik`. Its elements are
+#' \eqn{klas(i) = k}, \eqn{k = 1,\dots,K}.
+#' @field Ey_k Matrix of dimension \emph{(n,K)}.
+#' @field Ey Column matrix of dimension \emph{n}.
+#' @field Var_yk Column matrix of dimension \emph{K}.
+#' @field Vary Column matrix of dimension \emph{n}.
+#' @field log\_lik Numeric. Log-likelihood of the SNMoE model.
+#' @field com_loglik Numeric. Complete log-likelihood of the SNMoE model.
+#' @field stored_loglik List. Stored values of the log-likelihood at each EM
+#' iteration.
+#' @field BIC Numeric. Value of the BIC (Bayesian Information Criterion)
+#' criterion. The formula is \eqn{BIC = log\_lik - nu \times
+#' \textrm{log}(n) / 2}{BIC = log\_lik - nu x log(n) / 2} with \emph{nu} the
+#' degree of freedom of the SNMoE model.
+#' @field ICL Numeric. Value of the ICL (Integrated Completed Likelihood)
+#' criterion. The formula is \eqn{ICL = com\_loglik - nu \times
+#' \textrm{log}(n) / 2}{ICL = com_loglik - nu x log(n) / 2} with \emph{nu} the
+#' degree of freedom of the SNMoE model.
+#' @field AIC Numeric. Value of the AIC (Akaike Information Criterion)
+#' criterion. The formula is \eqn{AIC = log\_lik - nu}{AIC = log\_lik - nu}.
+#' @field log_piik_fik Matrix of size \eqn{(n, K)} giving the values of the
+#' logarithm of the joint probability
+#' \eqn{P(Y_{i}, \ zi = k)}{P(Yi, zi = k)}, \eqn{i = 1,\dots,n}.
+#' @field log_sum_piik_fik Column matrix of size \emph{n} giving the values of
+#' \eqn{\sum_{k = 1}^{K} \textrm{log} P(Y_{i}, \ zi = k)}{\sum_{k = 1}^{K} log
+#' P(Yi, zi = k)}, \eqn{i = 1,\dots,n}.
+#' @field tik Matrix of size \eqn{(n, K)} giving the posterior probability that
+#' \eqn{Y_{i}}{Yi} originates from the \eqn{k}-th regression model
+#' \eqn{P(zi = k | Y, W, \beta)}.
+#' @field E1ik To define.
+#' @field E2ik To define.
+#' @seealso [ParamSNMoE], [FData]
+#' @export
+StatSNMoE <- setRefClass(
+  "StatSNMoE",
+  fields = list(
+    piik = "matrix",
+    z_ik = "matrix",
+    klas = "matrix",
+    # Ex = "matrix",
+    Ey_k = "matrix",
+    Ey = "matrix",
+    Var_yk = "matrix",
+    Vary = "matrix",
+    log_lik = "numeric",
+    com_loglik = "numeric",
+    stored_loglik = "list",
+    BIC = "numeric",
+    ICL = "numeric",
+    AIC = "numeric",
+    log_piik_fik = "matrix",
+    log_sum_piik_fik = "matrix",
+    tik = "matrix",
+    E1ik = "matrix",
+    E2ik = "matrix"
+  ),
+  methods = list(
+    initialize = function(paramSNMoE = ParamSNMoE()) {
+      piik <<- matrix(NA, paramSNMoE$fData$n, paramSNMoE$K)
+      z_ik <<- matrix(NA, paramSNMoE$fData$n, paramSNMoE$K)
+      klas <<- matrix(NA, paramSNMoE$fData$n, 1)
+      Ey_k <<- matrix(NA, paramSNMoE$fData$n, paramSNMoE$K)
+      Ey <<- matrix(NA, paramSNMoE$fData$n, 1)
+      Var_yk <<- matrix(NA, 1, paramSNMoE$K)
+      Vary <<- matrix(NA, paramSNMoE$fData$n, 1)
+      log_lik <<- -Inf
+      com_loglik <<- -Inf
+      stored_loglik <<- list()
+      BIC <<- -Inf
+      ICL <<- -Inf
+      AIC <<- -Inf
+      log_piik_fik <<- matrix(0, paramSNMoE$fData$n, paramSNMoE$K)
+      log_sum_piik_fik <<- matrix(NA, paramSNMoE$fData$n, 1)
+      tik <<- matrix(0, paramSNMoE$fData$n, paramSNMoE$K)
+      E1ik <<- matrix(0, paramSNMoE$fData$m * paramSNMoE$fData$n, paramSNMoE$K)
+      E2ik <<- matrix(0, paramSNMoE$fData$m * paramSNMoE$fData$n, paramSNMoE$K)
+    },
+
+    MAP = function() {
+      "
+      calcule une partition d'un echantillon par la regle du Maximum A Posteriori ?? partir des probabilites a posteriori
+      Entrees : post_probas , Matrice de dimensions [n x K] des probabibiltes a posteriori (matrice de la partition floue)
+      n : taille de l'echantillon
+      K : nombres de classes
+      klas(i) = arg   max (post_probas(i,k)) , for all i=1,...,n
+      1<=k<=K
+      = arg   max  p(zi=k|xi;theta)
+      1<=k<=K
+      = arg   max  p(zi=k;theta)p(xi|zi=k;theta)/sum{l=1}^{K}p(zi=l;theta) p(xi|zi=l;theta)
+      1<=k<=K
+      Sorties : classes : vecteur collones contenant les classe (1:K)
+      Z : Matrice de dimension [nxK] de la partition dure : ses elements sont zik, avec zik=1 si xi
+      appartient ?? la classe k (au sens du MAP) et zero sinon.
+      "
+      N <- nrow(piik)
+      K <- ncol(piik)
+      ikmax <- max.col(piik)
+      ikmax <- matrix(ikmax, ncol = 1)
+      z_ik <<- ikmax %*% ones(1, K) == ones(N, 1) %*% (1:K) # partition_MAP
+      klas <<- ones(N, 1)
+      for (k in 1:K) {
+        klas[z_ik[, k] == 1] <<- k
+      }
+    },
+    #######
+    # compute loglikelihood
+    #######
+    computeLikelihood = function(reg_irls) {
+      log_lik <<- sum(log_sum_piik_fik) + reg_irls
+
+    },
+    #######
+    #
+    #######
+    #######
+    # compute the final solution stats
+    #######
+    computeStats = function(paramSNMoE) {
+
+      # E[yi|zi=k]
+      Ey_k <<- paramSNMoE$phiBeta$XBeta[1:paramSNMoE$fData$n, ] %*% paramSNMoE$beta + ones(paramSNMoE$fData$n, 1) %*% (sqrt(2 / pi) * paramSNMoE$delta * paramSNMoE$sigma)
+
+      # E[yi]
+      Ey <<- matrix(apply(piik * Ey_k, 1, sum))
+
+      # Var[yi|zi=k]
+      Var_yk <<- (1 - (2 / pi) * (paramSNMoE$delta ^ 2)) * (paramSNMoE$sigma ^ 2)
+
+      # Var[yi]
+      Vary <<- apply(piik * (Ey_k ^ 2 + ones(paramSNMoE$fData$n, 1) %*% Var_yk), 1, sum) - Ey ^ 2
+
+
+      ### BIC AIC et ICL
+
+      BIC <<- log_lik - (paramSNMoE$nu * log(paramSNMoE$fData$n * paramSNMoE$fData$m) / 2)
+      AIC <<- log_lik - paramSNMoE$nu
+      ## CL(theta) : complete-data loglikelihood
+      zik_log_piik_fk <- (repmat(z_ik, paramSNMoE$fData$m, 1)) * log_piik_fik
+      sum_zik_log_fik <- apply(zik_log_piik_fk, 1, sum)
+      com_loglik <<- sum(sum_zik_log_fik)
+
+      ICL <<- com_loglik - (paramSNMoE$nu * log(paramSNMoE$fData$n * paramSNMoE$fData$m) / 2)
+      # solution.XBeta = XBeta(1:m,:);
+      # solution.XAlpha = XAlpha(1:m,:);
+    },
+    #######
+    # EStep
+    #######
+    EStep = function(paramSNMoE) {
+      "Method used in the EM algorithm to update statistics based on parameters
+      provided by \\code{paramSNMoE} (prior and posterior probabilities)."
+      piik <<- multinomialLogit(paramSNMoE$alpha, paramSNMoE$phiAlpha$XBeta, ones(paramSNMoE$fData$n, paramSNMoE$K), ones(paramSNMoE$fData$n, 1))$piik
+
+      piik_fik <- zeros(paramSNMoE$fData$m * paramSNMoE$fData$n, paramSNMoE$K)
+
+      for (k in (1:paramSNMoE$K)) {
+        muk <- paramSNMoE$phiBeta$XBeta %*% paramSNMoE$beta[, k]
+
+        sigma2k <- paramSNMoE$sigma[k]
+        sigmak <- sqrt(sigma2k)
+        dik <- (paramSNMoE$fData$Y - muk) / sigmak
+
+        mu_uk <- (paramSNMoE$delta[k] * abs(paramSNMoE$fData$Y - muk))
+        sigma2_uk <- (1 - paramSNMoE$delta[k] ^ 2) * paramSNMoE$sigma[k]
+        sigma_uk <- sqrt(sigma2_uk)
+
+        E1ik[, k] <<- mu_uk + sigma_uk * dnorm(paramSNMoE$lambda[k] * dik, 0, 1) / pnorm(paramSNMoE$lambda[k] * dik, 0, 1)
+        E2ik[, k] <<- mu_uk ^ 2 + sigma_uk ^ 2 + sigma_uk * mu_uk * dnorm(paramSNMoE$lambda[k] * dik, 0, 1) / pnorm(paramSNMoE$lambda[k] * dik, 0, 1)
+
+        # weighted skew normal linear expert likelihood
+        piik_fik[, k] <- piik[, k] * (2 / sigmak) * dnorm(dik, 0, 1) * pnorm(paramSNMoE$lambda[k] * dik)
+      }
+
+      log_piik_fik <<- log(piik_fik)
+
+      log_sum_piik_fik <<- matrix(log(rowSums(piik_fik)))
+
+      tik <<- piik_fik / (rowSums(piik_fik) %*% ones(1, paramSNMoE$K))
+    }
+  )
+)
diff --git a/R/StatStMoE.R b/R/StatStMoE.R
new file mode 100755
index 0000000..98c18ad
--- /dev/null
+++ b/R/StatStMoE.R
@@ -0,0 +1,244 @@
+#' A Reference Class which contains statistics of a StMoE model.
+#'
+#' StatMRHLP contains all the parameters of a [StMoE][ParamStMoE] model.
+#'
+#' @field piik Matrix of size \eqn{(n, K)} representing the probabilities
+#' \eqn{P(zi = k; W) = P(z_{ik} = 1; W)}{P(zi = k; W) = P(z_ik = 1; W)} of the
+#' latent variable \eqn{zi,\ i = 1,\dots,m}{zi, i = 1,\dots,n}.
+#' @field z_ik Hard segmentation logical matrix of dimension \eqn{(n, K)}
+#' obtained by the Maximum a posteriori (MAP) rule:
+#' \eqn{z_{ik} = 1 \ \textrm{if} \ z_{ik} = \textrm{arg} \ \textrm{max}_{k} \
+#' P(z_i = k | Y, W, \beta);\ 0 \ \textrm{otherwise}}{z_ik = 1 if z_ik =
+#' arg max_k P(z_i = k | Y, W, \beta); 0 otherwise}, \eqn{k = 1,\dots,K}.
+#' @field klas Column matrix of the labels issued from `z_ik`. Its elements are
+#' \eqn{klas(i) = k}, \eqn{k = 1,\dots,K}.
+#' @field Ey_k Matrix of dimension \emph{(n,K)}.
+#' @field Ey Column matrix of dimension \emph{n}.
+#' @field Var_yk Column matrix of dimension \emph{K}.
+#' @field Var_y Column matrix of dimension \emph{n}.
+#' @field log\_lik Numeric. Log-likelihood of the StMoE model.
+#' @field com_loglik Numeric. Complete log-likelihood of the StMoE model.
+#' @field stored_loglik List. Stored values of the log-likelihood at each EM
+#' iteration.
+#' @field BIC Numeric. Value of the BIC (Bayesian Information Criterion)
+#' criterion. The formula is \eqn{BIC = log\_lik - nu \times
+#' \textrm{log}(n) / 2}{BIC = log\_lik - nu x log(n) / 2} with \emph{nu} the
+#' degree of freedom of the StMoE model.
+#' @field ICL Numeric. Value of the ICL (Integrated Completed Likelihood)
+#' criterion. The formula is \eqn{ICL = com\_loglik - nu \times
+#' \textrm{log}(n) / 2}{ICL = com_loglik - nu x log(n) / 2} with \emph{nu} the
+#' degree of freedom of the StMoE model.
+#' @field AIC Numeric. Value of the AIC (Akaike Information Criterion)
+#' criterion. The formula is \eqn{AIC = log\_lik - nu}{AIC = log\_lik - nu}.
+#' @field log_piik_fik Matrix of size \eqn{(n, K)} giving the values of the
+#' logarithm of the joint probability
+#' \eqn{P(Y_{i}, \ zi = k)}{P(Yi, zi = k)}, \eqn{i = 1,\dots,n}.
+#' @field log_sum_piik_fik Column matrix of size \emph{n} giving the values of
+#' \eqn{\sum_{k = 1}^{K} \textrm{log} P(Y_{i}, \ zi = k)}{\sum_{k = 1}^{K} log
+#' P(Yi, zi = k)}, \eqn{i = 1,\dots,n}.
+#' @field tik Matrix of size \eqn{(n, K)} giving the posterior probability that
+#' \eqn{Y_{i}}{Yi} originates from the \eqn{k}-th regression model
+#' \eqn{P(zi = k | Y, W, \beta)}.
+#' @field wik To define.
+#' @field dik To define.
+#' @field stme_pdf skew-t mixture of experts density
+#' @field E1ik To define.
+#' @field E2ik To define.
+#' @field E3ik To define.
+#' @seealso [ParamStMoE], [FData]
+#' @export
+StatStMoE <- setRefClass(
+  "StatStMoE",
+  fields = list(
+    piik = "matrix",
+    z_ik = "matrix",
+    klas = "matrix",
+    # Ex = "matrix",
+    Ey_k = "matrix",
+    Ey = "matrix",
+    Var_yk = "matrix",
+    Vary = "matrix",
+    log_lik = "numeric",
+    com_loglik = "numeric",
+    stored_loglik = "list",
+    BIC = "numeric",
+    ICL = "numeric",
+    AIC = "numeric",
+    log_piik_fik = "matrix",
+    log_sum_piik_fik = "matrix",
+    tik = "matrix",
+    wik = "matrix",
+    dik = "matrix",
+    stme_pdf = "matrix",
+    E1ik = "matrix",
+    E2ik = "matrix",
+    E3ik = "matrix"
+  ),
+  methods = list(
+    initialize = function(paramStMoE = ParamStMoE()) {
+      piik <<- matrix(NA, paramStMoE$fData$n, paramStMoE$K)
+      z_ik <<- matrix(NA, paramStMoE$fData$n, paramStMoE$K)
+      klas <<- matrix(NA, paramStMoE$fData$n, 1)
+      Ey_k <<- matrix(NA, paramStMoE$fData$n, paramStMoE$K)
+      Ey <<- matrix(NA, paramStMoE$fData$n, 1)
+      Var_yk <<- matrix(NA, 1, paramStMoE$K)
+      Vary <<- matrix(NA, paramStMoE$fData$n, 1)
+      log_lik <<- -Inf
+      com_loglik <<- -Inf
+      stored_loglik <<- list()
+      BIC <<- -Inf
+      ICL <<- -Inf
+      AIC <<- -Inf
+      log_piik_fik <<- matrix(0, paramStMoE$fData$n, paramStMoE$K)
+      log_sum_piik_fik <<- matrix(NA, paramStMoE$fData$n, 1)
+      tik <<- matrix(0, paramStMoE$fData$n, paramStMoE$K)
+      wik <<- matrix(0, paramStMoE$fData$m * paramStMoE$fData$n, paramStMoE$K)
+      dik <<- matrix(0, paramStMoE$fData$m * paramStMoE$fData$n, paramStMoE$K)
+      stme_pdf <<- matrix(0, paramStMoE$fData$n, 1)
+      E1ik <<- matrix(0, paramStMoE$fData$m * paramStMoE$fData$n, paramStMoE$K)
+      E2ik <<- matrix(0, paramStMoE$fData$m * paramStMoE$fData$n, paramStMoE$K)
+      E3ik <<- matrix(0, paramStMoE$fData$m * paramStMoE$fData$n, paramStMoE$K)
+    },
+
+    MAP = function() {
+      "
+      calcule une partition d'un echantillon par la regle du Maximum A Posteriori ?? partir des probabilites a posteriori
+      Entrees : post_probas , Matrice de dimensions [n x K] des probabibiltes a posteriori (matrice de la partition floue)
+      n : taille de l'echantillon
+      K : nombres de classes
+      klas(i) = arg   max (post_probas(i,k)) , for all i=1,...,n
+      1<=k<=K
+      = arg   max  p(zi=k|xi;theta)
+      1<=k<=K
+      = arg   max  p(zi=k;theta)p(xi|zi=k;theta)/sum{l=1}^{K}p(zi=l;theta) p(xi|zi=l;theta)
+      1<=k<=K
+      Sorties : classes : vecteur collones contenant les classe (1:K)
+      Z : Matrice de dimension [nxK] de la partition dure : ses elements sont zik, avec zik=1 si xi
+      appartient et la classe k (au sens du MAP) et zero sinon.
+      "
+      N <- nrow(piik)
+      K <- ncol(piik)
+      ikmax <- max.col(piik)
+      ikmax <- matrix(ikmax, ncol = 1)
+      z_ik <<- ikmax %*% ones(1, K) == ones(N, 1) %*% (1:K) # partition_MAP
+      klas <<- ones(N, 1)
+      for (k in 1:K) {
+        klas[z_ik[, k] == 1] <<- k
+      }
+    },
+    #######
+    # compute loglikelihood
+    #######
+    computeLikelihood = function(reg_irls) {
+      log_lik <<- sum(log_sum_piik_fik) + reg_irls
+
+    },
+    #######
+    #
+    #######
+    #######
+    # compute the final solution stats
+    #######
+    computeStats = function(paramStMoE) {
+
+      Xi_nuk = sqrt(paramStMoE$nuk/pi) * (gamma(paramStMoE$nuk/2 - 1/2)) / (gamma(paramStMoE$nuk/2));
+
+      # E[yi|zi=k]
+      Ey_k <<- paramStMoE$phiBeta$XBeta[1:paramStMoE$fData$n, ] %*% paramStMoE$beta + ones(paramStMoE$fData$n, 1) %*% (paramStMoE$delta * sqrt(paramStMoE$sigma) * Xi_nuk)
+
+      # E[yi]
+      Ey <<- matrix(apply(piik * Ey_k, 1, sum))
+
+      # Var[yi|zi=k]
+      Var_yk <<- (paramStMoE$nuk/(paramStMoE$nuk-2) - (paramStMoE$delta^2) * (Xi_nuk^2)) * paramStMoE$sigma
+
+      # Var[yi]
+      Vary <<- apply(piik * (Ey_k ^ 2 + ones(paramStMoE$fData$n, 1) %*% Var_yk), 1, sum) - Ey ^2
+
+
+      ### BIC AIC et ICL
+
+      BIC <<- log_lik - (paramStMoE$nu * log(paramStMoE$fData$n * paramStMoE$fData$m) / 2)
+      AIC <<- log_lik - paramStMoE$nu
+      ## CL(theta) : complete-data loglikelihood
+      zik_log_piik_fk <- (repmat(z_ik, paramStMoE$fData$m, 1)) * log_piik_fik
+      sum_zik_log_fik <- apply(zik_log_piik_fk, 1, sum)
+      com_loglik <<- sum(sum_zik_log_fik)
+
+      ICL <<- com_loglik - (paramStMoE$nu * log(paramStMoE$fData$n * paramStMoE$fData$m) / 2)
+      # solution.XBeta = XBeta(1:m,:);
+      # solution.XAlpha = XAlpha(1:m,:);
+    },
+
+    #######
+    # Intial value for STMoE density
+    #######
+    univStMoEpdf = function(paramStMoE){
+      piik <<- multinomialLogit(paramStMoE$alpha, paramStMoE$phiAlpha$XBeta, ones(paramStMoE$fData$n, paramStMoE$K), ones(paramStMoE$fData$n, 1))$piik
+
+      piik_fik <- zeros(paramStMoE$fData$n, paramStMoE$K)
+      dik <<- zeros(paramStMoE$fData$n, paramStMoE$K)
+      mik <- zeros(paramStMoE$fData$n, paramStMoE$K)
+
+      for (k in (1:paramStMoE$K)){
+        dik[,k] <<- (paramStMoE$fData$Y - paramStMoE$phiBeta$XBeta %*% paramStMoE$beta[, k])/ paramStMoE$sigma[k]
+        mik[,k] <- paramStMoE$lambda[k] %*% dik[,k] * sqrt(paramStMoE$nuk[k]+1)/(paramStMoE$nuk[k] + dik[,k]^2)
+
+        piik_fik[,k] <- piik[,k]*(2/paramStMoE$sigma[k])*dt(dik[,k], paramStMoE$nuk[k])*pt(mik[,k], paramStMoE$nuk[k]+1)
+      }
+
+      stme_pdf <<- matrix(rowSums(piik_fik)) #  skew-t mixture of experts density
+    },
+
+    #######
+    # EStep
+    #######
+    EStep = function(paramStMoE) {
+      "Method used in the EM algorithm to update statistics based on parameters
+      provided by \\code{paramStMoE} (prior and posterior probabilities)."
+      piik <<- multinomialLogit(paramStMoE$alpha, paramStMoE$phiAlpha$XBeta, ones(paramStMoE$fData$n, paramStMoE$K), ones(paramStMoE$fData$n, 1))$piik
+
+      piik_fik <- zeros(paramStMoE$fData$m * paramStMoE$fData$n, paramStMoE$K)
+
+      dik <<- zeros(paramStMoE$fData$m * paramStMoE$fData$n, paramStMoE$K)
+      mik <- zeros(paramStMoE$fData$m * paramStMoE$fData$n, paramStMoE$K)
+      wik <<- zeros(paramStMoE$fData$m * paramStMoE$fData$n, paramStMoE$K)
+
+
+      for (k in (1:paramStMoE$K)) {
+        muk <- paramStMoE$phiBeta$XBeta %*% paramStMoE$beta[, k]
+
+        sigma2k <- paramStMoE$sigma[k]
+        sigmak <- sqrt(sigma2k)
+        dik[,k] <<- (paramStMoE$fData$Y - muk) / sigmak
+
+        mik[,k] <- paramStMoE$lambda[k] %*% dik[,k] * sqrt((paramStMoE$nuk[k] +1)/(paramStMoE$nuk[k] + dik[,k]^2))
+
+        # E[Wi|yi,zik=1]
+        wik[,k] <<- ((paramStMoE$nuk[k] + 1)/(paramStMoE$nuk[k] + dik[,k]^2)) * pt(mik[,k]*sqrt((paramStMoE$nuk[k] + 3)/(paramStMoE$nuk[k] + 1)), paramStMoE$nuk[k] + 3)/pt(mik[,k], paramStMoE$nuk[k] + 1)
+
+        # E[Wi Ui |yi,zik=1]
+        deltak <- paramStMoE$delta[k]
+
+        E1ik[, k] <<- deltak * abs(paramStMoE$fData$Y - muk) * wik[,k] + (sqrt(1 - deltak^2)/(pi * stme_pdf)) * ((dik[,k]^2/(paramStMoE$nuk[k]*(1 - deltak^2)) + 1)^(-(paramStMoE$nuk[k]/2 + 1)))
+        E2ik[, k] <<- deltak^2 * ((paramStMoE$fData$Y - muk)^2) * wik[,k] + (1 - deltak^2) * sigmak^2 + ((deltak * (paramStMoE$fData$Y - muk) * sqrt(1 - deltak^2))/(pi * stme_pdf)) * (((dik[,k]^2)/(paramStMoE$nuk[k] * (1 - deltak^2)) + 1)^(-(paramStMoE$nuk[k]/2 + 1)))
+
+        Integgtx <- 0
+        E3ik[, k] <<- wik[,k] - log((paramStMoE$nuk[k] + dik[,k]^2)/2) -(paramStMoE$nuk[k] + 1)/(paramStMoE$nuk[k] + dik[,k]^2) + psigamma((paramStMoE$nuk[k] + 1)/2) + ((paramStMoE$lambda[k] * dik[,k] * (dik[,k]^2 - 1)) / sqrt((paramStMoE$nuk[k] + 1)*((paramStMoE$nuk[k] + dik[,k]^2)^3))) * dt(mik[,k], paramStMoE$nuk[k]+1) / pt(mik[,k], paramStMoE$nuk[k] + 1) + (1./pt(mik[,k], paramStMoE$nuk[k] + 1)) * Integgtx;
+
+        # weighted skew normal linear expert likelihood
+        piik_fik[, k] <- piik[, k] * (2 / sigmak) * dt(dik[, k], paramStMoE$nuk[k]) * pt(mik[,k], paramStMoE$nuk[k] + 1);
+      }
+
+
+      stme_pdf <<- matrix(rowSums(piik_fik)) # skew-t mixture of experts density
+
+      log_piik_fik <<- log(piik_fik)
+
+      log_sum_piik_fik <<- matrix(log(rowSums(piik_fik)))
+
+      #E[Zi=k|yi]
+      tik <<- piik_fik / (stme_pdf %*% ones(1,paramStMoE$K))
+    }
+  )
+)
diff --git a/R/StatTMoE.R b/R/StatTMoE.R
new file mode 100755
index 0000000..5797a58
--- /dev/null
+++ b/R/StatTMoE.R
@@ -0,0 +1,186 @@
+#' A Reference Class which contains statistics of a TMoE model.
+#'
+#' StatTMoE contains all the parameters of a [TMoE][ParamTMoE] model.
+#'
+#' @field piik Matrix of size \eqn{(n, K)} representing the probabilities
+#' \eqn{P(zi = k; W) = P(z_{ik} = 1; W)}{P(zi = k; W) = P(z_ik = 1; W)} of the
+#' latent variable \eqn{zi,\ i = 1,\dots,m}{zi, i = 1,\dots,n}.
+#' @field z_ik Hard segmentation logical matrix of dimension \eqn{(n, K)}
+#' obtained by the Maximum a posteriori (MAP) rule:
+#' \eqn{z_{ik} = 1 \ \textrm{if} \ z_{ik} = \textrm{arg} \ \textrm{max}_{k} \
+#' P(z_i = k | Y, W, \beta);\ 0 \ \textrm{otherwise}}{z_ik = 1 if z_ik =
+#' arg max_k P(z_i = k | Y, W, \beta); 0 otherwise}, \eqn{k = 1,\dots,K}.
+#' @field klas Column matrix of the labels issued from `z_ik`. Its elements are
+#' \eqn{klas(i) = k}, \eqn{k = 1,\dots,K}.
+#' @field Wik Matrix of dimension \emph{(nm,K)}.
+#' @field Ey_k Matrix of dimension \emph{(n,K)}.
+#' @field Ey Column matrix of dimension \emph{n}.
+#' @field Var_yk Column matrix of dimension \emph{K}.
+#' @field Vary Column matrix of dimension \emph{n}.
+#' @field log\_lik Numeric. Log-likelihood of the StMoE model.
+#' @field com_loglik Numeric. Complete log-likelihood of the StMoE model.
+#' @field stored_loglik List. Stored values of the log-likelihood at each EM
+#' iteration.
+#' @field BIC Numeric. Value of the BIC (Bayesian Information Criterion)
+#' criterion. The formula is \eqn{BIC = log\_lik - nu \times
+#' \textrm{log}(n) / 2}{BIC = log\_lik - nu x log(n) / 2} with \emph{nu} the
+#' degree of freedom of the StMoE model.
+#' @field ICL Numeric. Value of the ICL (Integrated Completed Likelihood)
+#' criterion. The formula is \eqn{ICL = com\_loglik - nu \times
+#' \textrm{log}(n) / 2}{ICL = com_loglik - nu x log(n) / 2} with \emph{nu} the
+#' degree of freedom of the StMoE model.
+#' @field AIC Numeric. Value of the AIC (Akaike Information Criterion)
+#' criterion. The formula is \eqn{AIC = log\_lik - nu}{AIC = log\_lik - nu}.
+#' @field log_piik_fik Matrix of size \eqn{(n, K)} giving the values of the
+#' logarithm of the joint probability
+#' \eqn{P(Y_{i}, \ zi = k)}{P(Yi, zi = k)}, \eqn{i = 1,\dots,n}.
+#' @field log_sum_piik_fik Column matrix of size \emph{n} giving the values of
+#' \eqn{\sum_{k = 1}^{K} \textrm{log} P(Y_{i}, \ zi = k)}{\sum_{k = 1}^{K} log
+#' P(Yi, zi = k)}, \eqn{i = 1,\dots,n}.
+#' @field tik Matrix of size \eqn{(n, K)} giving the posterior probability that
+#' \eqn{Y_{i}}{Yi} originates from the \eqn{k}-th regression model
+#' \eqn{P(zi = k | Y, W, \beta)}.
+#' @seealso [ParamStMoE], [FData]
+#' @export
+StatTMoE <- setRefClass(
+  "StatTMoE",
+  fields = list(
+    piik = "matrix",
+    z_ik = "matrix",
+    klas = "matrix",
+    Wik = "matrix",
+    # Ex = "matrix",
+    Ey_k = "matrix",
+    Ey = "matrix",
+    Var_yk = "matrix",
+    Vary = "matrix",
+    log_lik = "numeric",
+    com_loglik = "numeric",
+    stored_loglik = "list",
+    BIC = "numeric",
+    ICL = "numeric",
+    AIC = "numeric",
+    log_piik_fik = "matrix",
+    log_sum_piik_fik = "matrix",
+    tik = "matrix"
+  ),
+  methods = list(
+    initialize = function(paramTMoE = ParamTMoE()) {
+      piik <<- matrix(NA, paramTMoE$fData$n, paramTMoE$K)
+      z_ik <<- matrix(NA, paramTMoE$fData$n, paramTMoE$K)
+      klas <<- matrix(NA, paramTMoE$fData$n, 1)
+      Wik <<- matrix(0, paramTMoE$fData$n * paramTMoE$fData$m, paramTMoE$K)
+      Ey_k <<- matrix(NA, paramTMoE$fData$n, paramTMoE$K)
+      Ey <<- matrix(NA, paramTMoE$fData$n, 1)
+      Var_yk <<- matrix(NA, 1, paramTMoE$K)
+      Vary <<- matrix(NA, paramTMoE$fData$n, 1)
+      log_lik <<- -Inf
+      com_loglik <<- -Inf
+      stored_loglik <<- list()
+      BIC <<- -Inf
+      ICL <<- -Inf
+      AIC <<- -Inf
+      log_piik_fik <<- matrix(0, paramTMoE$fData$n, paramTMoE$K)
+      log_sum_piik_fik <<- matrix(NA, paramTMoE$fData$n, 1)
+      tik <<- matrix(0, paramTMoE$fData$n, paramTMoE$K)
+    },
+
+    MAP = function() {
+      "
+      calcule une partition d'un echantillon par la regle du Maximum A Posteriori ?? partir des probabilites a posteriori
+      Entrees : post_probas , Matrice de dimensions [n x K] des probabibiltes a posteriori (matrice de la partition floue)
+      n : taille de l'echantillon
+      K : nombres de classes
+      klas(i) = arg   max (post_probas(i,k)) , for all i=1,...,n
+      1<=k<=K
+      = arg   max  p(zi=k|xi;theta)
+      1<=k<=K
+      = arg   max  p(zi=k;theta)p(xi|zi=k;theta)/sum{l=1}^{K}p(zi=l;theta) p(xi|zi=l;theta)
+      1<=k<=K
+      Sorties : classes : vecteur collones contenant les classe (1:K)
+      Z : Matrice de dimension [nxK] de la partition dure : ses elements sont zik, avec zik=1 si xi
+      appartient ?? la classe k (au sens du MAP) et zero sinon.
+      "
+      N <- nrow(piik)
+      K <- ncol(piik)
+      ikmax <- max.col(piik)
+      ikmax <- matrix(ikmax, ncol = 1)
+      z_ik <<- ikmax %*% ones(1, K) == ones(N, 1) %*% (1:K) # partition_MAP
+      klas <<- ones(N, 1)
+      for (k in 1:K) {
+        klas[z_ik[, k] == 1] <<- k
+      }
+    },
+    #######
+    # compute loglikelihood
+    #######
+    computeLikelihood = function(reg_irls) {
+      log_lik <<- sum(log_sum_piik_fik) + reg_irls
+
+    },
+    #######
+    #
+    #######
+    #######
+    # compute the final solution stats
+    #######
+    computeStats = function(paramTMoE) {
+
+      # E[yi|zi=k]
+      Ey_k <<- paramTMoE$phiBeta$XBeta[1:paramTMoE$fData$n, ] %*% paramTMoE$beta
+
+      # E[yi]
+      Ey <<- matrix(apply(piik * Ey_k, 1, sum))
+
+      # Var[yi|zi=k]
+      Var_yk <<- paramTMoE$delta/(paramTMoE$delta - 2) * paramTMoE$sigma
+
+      # Var[yi]
+      Vary <<- apply(piik * (Ey_k ^ 2 + ones(paramTMoE$fData$n, 1) %*% Var_yk), 1, sum) - Ey ^ 2
+
+
+      ### BIC AIC et ICL
+
+      BIC <<- log_lik - (paramTMoE$nu * log(paramTMoE$fData$n * paramTMoE$fData$m) / 2)
+      AIC <<- log_lik - paramTMoE$nu
+      ## CL(theta) : complete-data loglikelihood
+      zik_log_piik_fk <- (repmat(z_ik, paramTMoE$fData$m, 1)) * log_piik_fik
+      sum_zik_log_fik <- apply(zik_log_piik_fk, 1, sum)
+      com_loglik <<- sum(sum_zik_log_fik)
+
+      ICL <<- com_loglik - (paramTMoE$nu * log(paramTMoE$fData$n * paramTMoE$fData$m) / 2)
+      # solution.XBeta = XBeta(1:m,:);
+      # solution.XAlpha = XAlpha(1:m,:);
+    },
+    #######
+    # EStep
+    #######
+    EStep = function(paramTMoE) {
+      piik <<- multinomialLogit(paramTMoE$alpha, paramTMoE$phiAlpha$XBeta, ones(paramTMoE$fData$n, paramTMoE$K), ones(paramTMoE$fData$n, 1))$piik
+      piik_fik <- zeros(paramTMoE$fData$m * paramTMoE$fData$n, paramTMoE$K)
+
+      for (k in (1:paramTMoE$K)) {
+        muk <- paramTMoE$phiBeta$XBeta %*% paramTMoE$beta[, k]
+
+        sigma2k <- paramTMoE$sigma[k]
+        sigmak <- sqrt(sigma2k)
+        dik <- (paramTMoE$fData$Y - muk) / sigmak
+
+
+        nuk <- paramTMoE$delta[k]
+        Wik[,k] <<- (nuk + 1)/(nuk + dik^2)
+
+
+        # weighted t linear expert likelihood
+
+        piik_fik[, k] <- piik[, k] *  (1/sigmak * dt((paramTMoE$fData$Y - muk)/sigmak, nuk)) #pdf('tlocationscale', y, muk, sigmak, nuk);
+      }
+
+      log_piik_fik <<- log(piik_fik)
+
+      log_sum_piik_fik <<- matrix(log(rowSums(piik_fik)))
+
+      tik <<- piik_fik / (rowSums(piik_fik) %*% ones(1, paramTMoE$K))
+    }
+  )
+)
diff --git a/R/emSNMoE.R b/R/emSNMoE.R
new file mode 100755
index 0000000..ef064c4
--- /dev/null
+++ b/R/emSNMoE.R
@@ -0,0 +1,118 @@
+#' emSNMoE is used to fit a SNMoE model.
+#'
+#' emSNMoE is used to fit a SNMoE model. The estimation method is performed by
+#' the Expectation-Maximization algorithm.
+#'
+#' @details emSNMoE function is based on the EM algorithm. This functions starts
+#' with an initialization of the parameters done by the method `initParam` of
+#' the class [ParamSNMoE][ParamSNMoE], then it alternates between a E-Step
+#' (method of the class [StatSNMoE][StatSNMoE]) and a M-Step (method of the class
+#' [ParamSNMoE][ParamSNMoE]) until convergence (until the absolute difference of
+#' log-likelihood between two steps of the EM algorithm is less than the
+#' `threshold` parameter).
+#'
+#' @param X Numeric vector of length \emph{m} representing the covariates.
+#' @param Y Matrix of size \eqn{(n, m)} representing \emph{n} functions of `X`
+#' observed at points \eqn{1,\dots,m}.
+#' @param K The number of mixture components.
+#' @param p The order of the polynomial regression.
+#' @param q The dimension of the logistic regression. For the purpose of
+#' segmentation, it must be set to 1.
+#' @param n_tries Number of times EM algorithm will be launched.
+#' The solution providing the highest log-likelihood will be returned.
+#'
+#' If `n_tries` > 1, then for the first pass, parameters are initialized
+#' by uniformly segmenting the data into K segments, and for the next passes,
+#' parameters are initialized by randomly segmenting the data into K contiguous
+#'  segments.
+#' @param max_iter The maximum number of iterations for the EM algorithm.
+#' @param threshold A numeric value specifying the threshold for the relative
+#'  difference of log-likelihood between two steps  of the EM as stopping
+#'  criteria.
+#' @param verbose A logical value indicating whether values of the
+#' log-likelihood should be printed during EM iterations.
+#' @param verbose_IRLS A logical value indicating whether values of the
+#' criterion optimized by IRLS should be printed at each step of the EM
+#' algorithm.
+#' @return EM returns an object of class [ModelSNMoE][ModelSNMoE].
+#' @seealso [ModelSNMoE], [ParamSNMoE], [StatSNMoE]
+#' @export
+emSNMoE <- function(X, Y, K, p = 3, q = 1, n_tries = 1, max_iter = 1500, threshold = 1e-6, verbose = FALSE, verbose_IRLS = FALSE) {
+
+    fData <- FData(X, Y)
+
+    top <- 0
+    try_EM <- 0
+    best_loglik <- -Inf
+
+    while (try_EM < n_tries) {
+      try_EM <- try_EM + 1
+      message("EM try nr ", try_EM)
+
+      # Initializations
+      param <- ParamSNMoE$new(fData = fData, K = K, p = p, q = q)
+      param$initParam(try_EM, segmental = TRUE)
+
+
+
+      iter <- 0
+      converge <- FALSE
+      prev_loglik <- -Inf
+
+      stat <- StatSNMoE(paramSNMoE = param)
+
+      while (!converge && (iter <= max_iter)) {
+        stat$EStep(param)
+
+        reg_irls <- param$MStep(stat, verbose_IRLS)
+
+        stat$computeLikelihood(reg_irls)
+        # FIN EM
+
+        iter <- iter + 1
+        if (verbose) {
+          message("EM : Iteration : ", iter," log-likelihood : "  , stat$log_lik)
+        }
+        if (prev_loglik - stat$log_lik > 1e-5) {
+          message("!!!!! EM log-likelihood is decreasing from ", prev_loglik, "to ", stat$log_lik)
+          top <- top + 1
+          if (top > 20)
+            break
+        }
+
+        # TEST OF CONVERGENCE
+        converge <- abs((stat$log_lik - prev_loglik) / prev_loglik) <= threshold
+        if (is.na(converge)) {
+          converge <- FALSE
+        } # Basically for the first iteration when prev_loglik is Inf
+
+        prev_loglik <- stat$log_lik
+        stat$stored_loglik[iter] <- stat$log_lik
+      }# FIN EM LOOP
+
+      # at this point we have computed param and stat that contains all the information
+
+      if (stat$log_lik > best_loglik) {
+        statSolution <- stat$copy()
+        paramSolution <- param$copy()
+
+        best_loglik <- stat$log_lik
+      }
+      if (n_tries > 1) {
+        message("max value: ", stat$log_lik)
+      }
+    }
+
+    # Computation of c_ig the hard partition of the curves and klas
+    statSolution$MAP()
+
+    if (n_tries > 1) {
+      message("max value: ", statSolution$log_lik)
+    }
+
+
+    # FINISH computation of statSolution
+    statSolution$computeStats(paramSolution)
+
+    return(ModelSNMoE(param = paramSolution, stat = statSolution))
+  }
diff --git a/R/emStMoE.R b/R/emStMoE.R
new file mode 100755
index 0000000..54bc80e
--- /dev/null
+++ b/R/emStMoE.R
@@ -0,0 +1,122 @@
+#' emStMoE is used to fit a StMoE model.
+#'
+#' emStMoE is used to fit a StMoE model. The estimation method is performed by
+#' the Expectation-Maximization algorithm.
+#'
+#' @details emStMoE function is based on the EM algorithm. This functions starts
+#' with an initialization of the parameters done by the method `initParam` of
+#' the class [ParamStMoE][ParamStMoE], then it alternates between a E-Step
+#' (method of the class [StatStMoE][StatStMoE]) and a M-Step (method of the class
+#' [ParamStMoE][ParamStMoE]) until convergence (until the absolute difference of
+#' log-likelihood between two steps of the EM algorithm is less than the
+#' `threshold` parameter).
+#'
+#' @param X Numeric vector of length \emph{m} representing the covariates.
+#' @param Y Matrix of size \eqn{(n, m)} representing \emph{n} functions of `X`
+#' observed at points \eqn{1,\dots,m}.
+#' @param K The number of mixture components.
+#' @param p The order of the polynomial regression.
+#' @param q The dimension of the logistic regression. For the purpose of
+#' segmentation, it must be set to 1.
+#' @param variance_type Numeric indicating if the model is homoskedastic
+#' (`variance_type` = 1) or heteroskedastic (`variance_type` = 2).
+#' @param n_tries Number of times EM algorithm will be launched.
+#' The solution providing the highest log-likelihood will be returned.
+#'
+#' If `n_tries` > 1, then for the first pass, parameters are initialized
+#' by uniformly segmenting the data into K segments, and for the next passes,
+#' parameters are initialized by randomly segmenting the data into K contiguous
+#'  segments.
+#' @param max_iter The maximum number of iterations for the EM algorithm.
+#' @param threshold A numeric value specifying the threshold for the relative
+#'  difference of log-likelihood between two steps  of the EM as stopping
+#'  criteria.
+#' @param verbose A logical value indicating whether values of the
+#' log-likelihood should be printed during EM iterations.
+#' @param verbose_IRLS A logical value indicating whether values of the
+#' criterion optimized by IRLS should be printed at each step of the EM
+#' algorithm.
+#' @return EM returns an object of class [ModelStMoE][ModelStMoE].
+#' @seealso [ModelStMoE], [ParamStMoE], [StatStMoE]
+#' @export
+emStMoE <- function(X, Y, K, p = 3, q = 1, n_tries = 1, max_iter = 1500, threshold = 1e-6, verbose = FALSE, verbose_IRLS = FALSE) {
+
+    fData <- FData(X, Y)
+
+    top <- 0
+    try_EM <- 0
+    best_loglik <- -Inf
+
+    while (try_EM < n_tries) {
+      try_EM <- try_EM + 1
+      message("EM try nr ", try_EM)
+
+      # Initializations
+      param <- ParamStMoE$new(fData = fData, K = K, p = p, q = q)
+      param$initParam(try_EM, segmental = FALSE)
+
+
+
+      iter <- 0
+      converge <- FALSE
+      prev_loglik <- -Inf
+
+      stat <- StatStMoE(paramStMoE = param)
+      stat$univStMoEpdf(param)
+
+      while (!converge && (iter <= max_iter)) {
+        stat$EStep(param)
+
+        reg_irls <- param$MStep(stat, verbose_IRLS)
+
+        stat$computeLikelihood(reg_irls)
+        # FIN EM
+
+        iter <- iter + 1
+        if (verbose) {
+          message("EM : Iteration : ", iter," log-likelihood : "  , stat$log_lik)
+        }
+        if (prev_loglik - stat$log_lik > 1e-5) {
+          message("!!!!! EM log-likelihood is decreasing from ", prev_loglik, "to ", stat$log_lik)
+          top <- top + 1
+          if (top > 20)
+            break
+        }
+
+        # TEST OF CONVERGENCE
+        converge <- abs((stat$log_lik - prev_loglik) / prev_loglik) <= threshold
+        if (is.na(converge)) {
+          converge <- FALSE
+        } # Basically for the first iteration when prev_loglik is -Inf
+
+        prev_loglik <- stat$log_lik
+        stat$stored_loglik[iter] <- stat$log_lik
+      }# FIN EM LOOP
+
+      # end of computation of all estimates (param and stat)
+
+      if (stat$log_lik > best_loglik) {
+        statSolution <- stat$copy()
+        paramSolution <- param$copy()
+
+        best_loglik <- stat$log_lik
+      }
+      if (n_tries > 1) {
+        message("max value: ", stat$log_lik)
+      }
+    }
+
+    # Computation of c_ig the hard partition of the curves and klas
+    statSolution$MAP()
+
+    if (n_tries > 1) {
+      message("max value: ", statSolution$log_lik)
+    }
+
+
+    # FINISH computation of statSolution
+    statSolution$computeStats(paramSolution)
+
+    return(ModelStMoE(param = paramSolution, stat = statSolution))
+
+  }
diff --git a/R/emTMoE.R b/R/emTMoE.R
new file mode 100755
index 0000000..297f1dd
--- /dev/null
+++ b/R/emTMoE.R
@@ -0,0 +1,121 @@
+#' emTMoE is used to fit a TMoE model.
+#'
+#' emTMoE is used to fit a TMoE model. The estimation method is performed by
+#' the Expectation-Maximization algorithm.
+#'
+#' @details emStMoE function is based on the EM algorithm. This functions starts
+#' with an initialization of the parameters done by the method `initParam` of
+#' the class [ParamTMoE][ParamTMoE], then it alternates between a E-Step
+#' (method of the class [StatTMoE][StatTMoE]) and a M-Step (method of the class
+#' [ParamTMoE][ParamTMoE]) until convergence (until the absolute difference of
+#' log-likelihood between two steps of the EM algorithm is less than the
+#' `threshold` parameter).
+#'
+#' @param X Numeric vector of length \emph{m} representing the covariates.
+#' @param Y Matrix of size \eqn{(n, m)} representing \emph{n} functions of `X`
+#' observed at points \eqn{1,\dots,m}.
+#' @param K The number of mixture components.
+#' @param p The order of the polynomial regression.
+#' @param q The dimension of the logistic regression. For the purpose of
+#' segmentation, it must be set to 1.
+#' @param variance_type Numeric indicating if the model is homoskedastic
+#' (`variance_type` = 1) or heteroskedastic (`variance_type` = 2).
+#' @param n_tries Number of times EM algorithm will be launched.
+#' The solution providing the highest log-likelihood will be returned.
+#'
+#' If `n_tries` > 1, then for the first pass, parameters are initialized
+#' by uniformly segmenting the data into K segments, and for the next passes,
+#' parameters are initialized by randomly segmenting the data into K contiguous
+#'  segments.
+#' @param max_iter The maximum number of iterations for the EM algorithm.
+#' @param threshold A numeric value specifying the threshold for the relative
+#'  difference of log-likelihood between two steps  of the EM as stopping
+#'  criteria.
+#' @param verbose A logical value indicating whether values of the
+#' log-likelihood should be printed during EM iterations.
+#' @param verbose_IRLS A logical value indicating whether values of the
+#' criterion optimized by IRLS should be printed at each step of the EM
+#' algorithm.
+#' @return EM returns an object of class [ModelTMoE][ModelTMoE].
+#' @seealso [ModelTMoE], [ParamTMoE], [StatTMoE]
+#' @export
+emTMoE <- function(X, Y, K, p = 3, q = 1, n_tries = 1, max_iter = 1500, threshold = 1e-6, verbose = FALSE, verbose_IRLS = FALSE) {
+
+    fData <- FData(X, Y)
+
+    top <- 0
+    try_EM <- 0
+    best_loglik <- -Inf
+
+    while (try_EM < n_tries) {
+      try_EM <- try_EM + 1
+      message("EM try nr ", try_EM)
+
+      # Initializations
+      param <- ParamTMoE$new(fData = fData, K = K, p = p, q = q)
+      param$initParam(try_EM, segmental = TRUE)
+
+
+
+      iter <- 0
+      converge <- FALSE
+      prev_loglik <- -Inf
+
+      stat <- StatTMoE(paramTMoE = param)
+
+      while (!converge && (iter <= max_iter)) {
+        stat$EStep(param)
+
+        reg_irls <- param$MStep(stat, verbose_IRLS)
+
+        stat$computeLikelihood(reg_irls)
+        # FIN EM
+
+        iter <- iter + 1
+        if (verbose) {
+          message("EM : Iteration : ", iter," log-likelihood : "  , stat$log_lik)
+        }
+        if (prev_loglik - stat$log_lik > 1e-5) {
+          message("!!!!! EM log-likelihood is decreasing from ", prev_loglik, "to ", stat$log_lik)
+          top <- top + 1
+          if (top > 20)
+            break
+        }
+
+        # TEST OF CONVERGENCE
+        converge <- abs((stat$log_lik - prev_loglik) / prev_loglik) <= threshold
+        if (is.na(converge)) {
+          converge <- FALSE
+        } # Basically for the first iteration when prev_loglik is Inf
+
+        prev_loglik <- stat$log_lik
+        stat$stored_loglik[iter] <- stat$log_lik
+      }# FIN EM LOOP
+
+      # at this point we have computed param and stat that contains all the information
+
+      if (stat$log_lik > best_loglik) {
+        statSolution <- stat$copy()
+        paramSolution <- param$copy()
+
+        best_loglik <- stat$log_lik
+      }
+      if (n_tries > 1) {
+        message("max value: ", stat$log_lik)
+      }
+    }
+
+    # Computation of c_ig the hard partition of the curves and klas
+    statSolution$MAP()
+
+    if (n_tries > 1) {
+      message("max value: ", statSolution$log_lik)
+    }
+
+
+    # FINISH computation of statSolution
+    statSolution$computeStats(paramSolution)
+
+    return(ModelTMoE(param = paramSolution, stat = statSolution))
+
+  }
diff --git a/R/meteorit-package.R b/R/meteorit-package.R
new file mode 100644
index 0000000..beae311
--- /dev/null
+++ b/R/meteorit-package.R
@@ -0,0 +1,9 @@
+#' @import methods
+## usethis namespace: start
+#' @useDynLib meteorit, .registration = TRUE
+## usethis namespace: end
+NULL
+## usethis namespace: start
+#' @importFrom Rcpp sourceCpp
+## usethis namespace: end
+NULL
diff --git a/R/utils.R b/R/utils.R
new file mode 100644
index 0000000..19c146c
--- /dev/null
+++ b/R/utils.R
@@ -0,0 +1,109 @@
+designmatrix = function(x, p, q = NULL, n = 1) {
+
+  order_max <- p
+  if (!is.null(q)) {
+    order_max <- max(p, q)
+  }
+
+  X <- matrix(NA, length(x), order_max + 1)
+  for (i in 0:(order_max)) {
+    X[, i + 1] <- x ^ i
+  }
+
+  XBeta <- X[, 1:(p + 1)]
+  # design matrix for Beta (the polynomial regressions)
+  if (!is.null(q)) {
+    Xw <- X[, 1:(q + 1)]
+    Xw <- repmat(Xw, n, 1)
+    # design matrix for w (the logistic regression)
+  } else {
+    Xw <- NULL
+  }
+
+  XBeta <- repmat(XBeta, n, 1)
+
+  return(list(Xw = Xw, XBeta = XBeta))
+}
+
+ones <- function(n, d, g = 1) {
+  if (g == 1) {
+    return(matrix(1, n, d))
+  }
+  else{
+    return(array(1, dim = c(n, d, g)))
+  }
+}
+
+zeros <- function(n, d, g = 1) {
+  if (g == 1) {
+    return(matrix(0, n, d))
+  }
+  else{
+    return(array(0, dim = c(n, d, g)))
+  }
+}
+
+rand <- function(n, d, g = 1) {
+  if (g == 1) {
+    return(matrix(stats::runif(n * d), n, d))
+  }
+  else{
+    return(array(stats::runif(n * d), dim = c(n, d, g)))
+  }
+}
+
+repmat <- function(M, n, d) {
+  return(kronecker(matrix(1, n, d), M))
+}
+
+drnorm <- function(n, d, mean, sd) {
+  A <- matrix(nrow = n, ncol = d)
+  for (i in 1:d) {
+    A[, i] <- stats::rnorm(n, mean, sd)
+  }
+  return(A)
+}
+
+lognormalize <- function(M) {
+  if (!is.matrix(M)) {
+    M <- matrix(M)
+  }
+  n <- nrow(M)
+  d <- ncol(M)
+  a <- apply(M, 1, max)
+  return(M - repmat(a + log(rowSums(exp(M - repmat(a, 1, d)))), 1, d))
+}
+
+normalize <- function(A, dim) {
+  # Normalize makes the entries of a (multidimensional <= 2) array sum to 1.
+  # Input
+  # A = Array to be normalized
+  # dim = dimension is specified to normalize.
+  # Output
+  # M = Array after normalize.
+  # z is the normalize constant
+  # Note:
+  # If dim is specified, we normalize the specified dimension only,
+  # Otherwise we normalize the whole array.
+  # Dim = 1 normalize each column
+  # Dim = 2 normalize each row
+
+  if (nargs() < 2) {
+    z <- sum(A)
+    # Set any zeros to one before dividing
+    # This is valid, since c = 0 ==> all i.A[i] = 0 ==> the anser should be 0/1 = 0.
+    s <- z + (z == 0)
+    M <- A / s
+  } else if (dim == 1) {
+    # normalize each column
+    z <- colSums(A)
+    s <- z + (z == 0)
+    M <- A / matrix(s, nrow = dim(A)[1], ncol = length(s), byrow = TRUE)
+  } else{
+    z <- rowSums(A)
+    s <- z + (z == 0)
+    M <- A / matrix(s, ncol = dim(A)[2], nrow = length(s), byrow = FALSE)
+  }
+  output <- list(M = M, z = z)
+  return(output)
+}
diff --git a/README.Rmd b/README.Rmd
new file mode 100644
index 0000000..12e2e11
--- /dev/null
+++ b/README.Rmd
@@ -0,0 +1,67 @@
+---
+output: github_document
+---
+
+<!-- README.md is generated from README.Rmd. Please edit that file -->
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  fig.align = "center",
+  fig.path = "man/figures/README-"
+)
+```
+## Overview
+
+<!-- badges: start -->
+<!-- badges: end -->
+
+User-friendly and flexible algorithm for modeling, sampling, inference, and clustering heteregenous data with the Skew-Normal Mixture-of-Experts (SNMoE) model.
+
+## Installation
+
+You can install the development version of STMoE from [GitHub](https://github.com/) with:
+
+```{r, eval = FALSE}
+# install.packages("devtools")
+devtools::install_github("fchamroukhi/SNMoE")
+```
+
+To build *vignettes* for examples of usage, type the command below instead:
+
+```{r, eval = FALSE}
+# install.packages("devtools")
+devtools::install_github("fchamroukhi/SNMoE", 
+                         build_opts = c("--no-resave-data", "--no-manual"), 
+                         build_vignettes = TRUE)
+```
+
+Use the following command to display vignettes:
+
+```{r, eval = FALSE}
+browseVignettes("SNMoE")
+```
+
+## Usage
+
+```{r, message = FALSE}
+library(SNMoE)
+
+data("simulatedstructureddata")
+
+K <- 2 # Number of regimes (mixture components)
+p <- 1 # Dimension of beta (order of the polynomial regressors)
+q <- 1 # Dimension of w (order of the logistic regression: to be set to 1 for segmentation)
+
+n_tries <- 1
+max_iter <- 1500
+threshold <- 1e-6
+verbose <- TRUE
+verbose_IRLS <- FALSE
+
+snmoe <- emSNMoE(simulatedstructureddata$X, matrix(simulatedstructureddata$Y), 
+                 K, p, q, n_tries, max_iter, threshold, verbose, verbose_IRLS)
+
+snmoe$plot()
+```
diff --git a/data/simulatedstructureddata.rda b/data/simulatedstructureddata.rda
new file mode 100644
index 0000000..d73dd12
Binary files /dev/null and b/data/simulatedstructureddata.rda differ
diff --git a/inst/CITATION b/inst/CITATION
new file mode 100644
index 0000000..3c2f17a
--- /dev/null
+++ b/inst/CITATION
@@ -0,0 +1,97 @@
+citHeader("To cite meteorit in publications use:")
+
+citEntry(
+  entry    = "Article",
+  title    = "Practical and theoretical aspects of mixture-of-experts modeling: An overview",
+  author = personList(person(given="Hien D.", family="Nguyen"),
+                      person(given="F.", family="Chamroukhi")),
+  journal  = "Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery",
+  publisher = "Wiley Periodicals, Inc",
+  year     = "2018",
+  pages = "e1246--n/a",
+  doi = "10.1002/widm.1246",
+  url      = "https://chamroukhi.com/papers/Nguyen-Chamroukhi-MoE-DMKD-2018",
+  textVersion = "Nguyen H and Chamroukhi F (2018). \u201cPractical and theoretical aspects
+  of mixture-of-experts modeling: An overview.\u201d Wiley Interdisciplinary Reviews:
+  Data Mining and Knowledge Discovery, pp. e1246-n/a.
+  doi: 10.1002/widm.1246 (URL:https://doi.org/10.1002/widm.1246),
+  <URL: https://chamroukhi.com/papers/Nguyen-Chamroukhi-MoE-DMKD-2018>."
+)
+
+citEntry(
+  entry    = "InProceedings",
+  title    = "Skew-Normal Mixture of Experts",
+  author = personList(person(given="F.", family="Chamroukhi")),
+  booktitle  = "The International Joint Conference on Neural Networks (IJCNN)",
+  year     = "2016",
+  url      = "https://chamroukhi.com/papers/Chamroukhi-SNMoE-IJCNN2016.pdf",
+  textVersion = "Chamroukhi F (2016). \u201cSkew-Normal Mixture of Experts.\u201d In
+  The International Joint Conference on Neural Networks (IJCNN).
+  <URL:https://chamroukhi.com/papers/Chamroukhi-SNMoE-IJCNN2016.pdf>."
+)
+
+citEntry(
+  entry    = "Misc",
+  title    = "Non-Normal Mixtures of Experts",
+  author = personList(person(given="F.", family="Chamroukhi")),
+  year     = "2015",
+  url      = "http://arxiv.org/pdf/1506.06707.pdf",
+  textVersion = "Chamroukhi F (2015). \u201cNon-Normal Mixtures of Experts.\u201d
+                <URL: http://arxiv.org/pdf/1506.06707.pdf>."
+)
+
+citEntry(
+  entry = "Article",
+  title = "Model-Based Clustering and Classification of Functional Data",
+  author = personList(person(given="F.", family="Chamroukhi"),
+                      person(given="Hien D.", family="Nguyen")),
+  journal = "Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery",
+  year = "2019",
+  url = "https://chamroukhi.com/papers/MBCC-FDA.pdf",
+  doi = "10.1002/widm.1298",
+  textVersion = "Chamroukhi F, Nguyen H (2019). \u201cModel-Based Clustering and
+  Classification of Functional Data.\u201d Wiley Interdisciplinary Reviews: Data Mining and
+  Knowledge Discovery. doi: 10.1002/widm.1298 (URL: https://doi.org/10.1002/widm.1298),
+  <URL:https://chamroukhi.com/papers/MBCC-FDA.pdf>."
+)
+
+citEntry(
+  entry    = "Article",
+  title    = "Skew t mixture of experts",
+  author   = personList(person(given="F.", family="Chamroukhi")),
+  journal  = "Neurocomputing - Elsevier",
+  year     = "2017",
+  volume   = "266",
+  pages    = "390-408",
+  url      = "https://chamroukhi.com/papers/STMoE.pdf",
+  textVersion = "Chamroukhi F (2017). \u201cSkew t mixture of experts.\u201d
+  Neurocomputing - Elsevier,
+  266, pp. 390-408. <URL:https://chamroukhi.com/papers/STMoE.pdf>."
+)
+
+citEntry(
+  entry = "Article",
+  title = "Robust mixture of experts modeling using the t-distribution",
+  author = personList(person(given="F.", family="Chamroukhi")),
+  journal = "Neural Networks - Elsevier",
+  year = "2016",
+  volume = "79",
+  pages = "20--36",
+  url = "https://chamroukhi.com/papers/TMoE.pdf",
+  textVersion = "Chamroukhi F (2016). \u201cRobust mixture of experts modeling using the
+  t-distribution.\u201d Neural Networks - Elsevier, 79, pp. 20-36.
+  <URL: https://chamroukhi.com/papers/TMoE.pdf>."
+)
+
+citEntry(
+  entry = "PhdThesis",
+  title = "Statistical learning of latent data models for complex data analysis",
+  author = person(given="F.", family="Chamroukhi"),
+  school = "Universit\u00e9 de Toulon",
+  year = "2015",
+  type = "Habilitation Thesis (HDR)",
+  url = "https://chamroukhi.com/Dossier/FChamroukhi-Habilitation.pdf",
+  textVersion = "Chamroukhi F (2015). Statistical learning of latent data models for
+  complex data analysis. Habilitation Thesis (HDR), Universit\u00e9 de Toulon.
+  <URL: https://chamroukhi.com/Dossier/FChamroukhi-Habilitation.pdf>."
+)
diff --git a/man/FData-class.Rd b/man/FData-class.Rd
new file mode 100644
index 0000000..eba5cf5
--- /dev/null
+++ b/man/FData-class.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/FData.R
+\docType{class}
+\name{FData-class}
+\alias{FData-class}
+\alias{FData}
+\title{A Reference Class which represents functional data.}
+\description{
+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 \code{X}
+represents the time.
+}
+\section{Fields}{
+
+\describe{
+\item{\code{X}}{Numeric vector of length \emph{m}.}
+
+\item{\code{Y}}{Matrix of size \eqn{(n, m)} representing \emph{n}
+functions of \code{X} observed at points \eqn{1,\dots,m}.}
+}}
+
+
diff --git a/man/ModelSNMoE-class.Rd b/man/ModelSNMoE-class.Rd
new file mode 100644
index 0000000..a9777af
--- /dev/null
+++ b/man/ModelSNMoE-class.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ModelSNMoE.R
+\docType{class}
+\name{ModelSNMoE-class}
+\alias{ModelSNMoE-class}
+\alias{ModelSNMoE}
+\title{A Reference Class which represents a fitted SNMoE model.}
+\description{
+ModelSNMoE represents a \link[=ModelSNMoE]{SNMoE} model for which parameters have
+been estimated.
+}
+\section{Fields}{
+
+\describe{
+\item{\code{param}}{A \link{ParamSNMoE} object. It contains the estimated values of the parameters.}
+
+\item{\code{stat}}{A \link{StatSNMoE} object. It contains all the statistics associated to the SNMoE model.}
+}}
+
+
+\seealso{
+\link{ParamSNMoE}, \link{StatSNMoE}
+}
diff --git a/man/ModelStMoE-class.Rd b/man/ModelStMoE-class.Rd
new file mode 100644
index 0000000..fd91cfe
--- /dev/null
+++ b/man/ModelStMoE-class.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ModelStMoE.R
+\docType{class}
+\name{ModelStMoE-class}
+\alias{ModelStMoE-class}
+\alias{ModelStMoE}
+\title{A Reference Class which represents a fitted StMoE model.}
+\description{
+ModelStMoE represents a \link[=ModelStMoE]{StMoE} model for which parameters have
+been estimated.
+}
+\section{Fields}{
+
+\describe{
+\item{\code{param}}{A \link{ParamStMoE} object. It contains the estimated values of the parameters.}
+
+\item{\code{stat}}{A \link{StatStMoE} object. It contains all the statistics associated to the StMoE model.}
+}}
+
+
+\seealso{
+\link{ParamStMoE}, \link{StatStMoE}
+}
diff --git a/man/ModelTMoE-class.Rd b/man/ModelTMoE-class.Rd
new file mode 100644
index 0000000..592c72b
--- /dev/null
+++ b/man/ModelTMoE-class.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ModelTMoE.R
+\docType{class}
+\name{ModelTMoE-class}
+\alias{ModelTMoE-class}
+\alias{ModelTMoE}
+\title{A Reference Class which represents a fitted TMoE model.}
+\description{
+ModelMRHLP represents a \link[=ModelTMoE]{TMoE} model for which parameters have
+been estimated.
+}
+\section{Fields}{
+
+\describe{
+\item{\code{param}}{A \link{ParamTMoE} object. It contains the estimated values of the parameters.}
+
+\item{\code{stat}}{A \link{StatTMoE} object. It contains all the statistics associated to the TMoE model.}
+}}
+
+
+\seealso{
+\link{ParamTMoE}, \link{StatTMoE}
+}
diff --git a/man/ParamSNMoE-class.Rd b/man/ParamSNMoE-class.Rd
new file mode 100644
index 0000000..4d14e01
--- /dev/null
+++ b/man/ParamSNMoE-class.Rd
@@ -0,0 +1,41 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ParamSNMoE.R
+\docType{class}
+\name{ParamSNMoE-class}
+\alias{ParamSNMoE-class}
+\alias{ParamSNMoE}
+\title{A Reference Class which contains parameters of a SNMoE model.}
+\description{
+ParamSNMoE contains all the parameters of a SNMoE model.
+}
+\section{Fields}{
+
+\describe{
+\item{\code{fData}}{\link{FData} object representing the sample.}
+
+\item{\code{K}}{The number of mixture components.}
+
+\item{\code{p}}{The order of the polynomial regression.}
+
+\item{\code{q}}{The dimension of the logistic regression. For the purpose of
+segmentation, it must be set to 1.}
+
+\item{\code{nu}}{degree of freedom}
+
+\item{\code{alpha}}{is the parameter vector of the logistic model with \eqn{alpha_K} being the null vector.}
+
+\item{\code{beta}}{is the vector of regression coefficients of component k,
+the updates for each of the expert component parameters consist in analytically solving a weighted
+Gaussian linear regression problem.}
+
+\item{\code{sigma}}{The variances for the \emph{K} mixture component.}
+
+\item{\code{lambda}}{skewness parameter}
+
+\item{\code{delta}}{the skewness parameter lambda (by equivalence delta)}
+}}
+
+
+\seealso{
+\link{FData}
+}
diff --git a/man/ParamStMoE-class.Rd b/man/ParamStMoE-class.Rd
new file mode 100644
index 0000000..f916929
--- /dev/null
+++ b/man/ParamStMoE-class.Rd
@@ -0,0 +1,52 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ParamStMoE.R
+\docType{class}
+\name{ParamStMoE-class}
+\alias{ParamStMoE-class}
+\alias{ParamStMoE}
+\title{A Reference Class which contains parameters of a MRHLP model.}
+\description{
+ParamMRHLP contains all the parameters of a MRHLP model.
+}
+\section{Fields}{
+
+\describe{
+\item{\code{fData}}{\link{FData} object representing the sample.}
+
+\item{\code{K}}{The number of mixture components.}
+
+\item{\code{p}}{The order of the polynomial regression.}
+
+\item{\code{q}}{The dimension of the logistic regression. For the purpose of
+segmentation, it must be set to 1.}
+
+\item{\code{nu}}{degree of freedom}
+
+\item{\code{alpha}}{is the parameter vector of the logistic model with \eqn{alpha_K} being the null vector.}
+
+\item{\code{beta}}{is the vector of regression coefficients of component k,
+the updates for each of the expert component parameters consist in analytically solving a weighted
+Gaussian linear regression problem.}
+
+\item{\code{sigma}}{The variances for the \emph{K} mixture component.}
+
+\item{\code{lambda}}{skewness parameter}
+
+\item{\code{delta}}{the skewness parameter lambda (by equivalence delta)}
+
+\item{\code{nuk}}{degrees of freedom}
+}}
+
+\section{Methods}{
+
+\describe{
+\item{\code{initParam(try_EM, segmental = FALSE)}}{Method to initialize parameters \code{alpha}, \code{beta} and
+\code{sigma}.}
+
+\item{\code{MStep(statStMoE, verbose_IRLS)}}{Method used in the EM algorithm to learn the parameters of the StMoE model
+based on statistics provided by \code{statStMoE}.}
+}}
+
+\seealso{
+\link{FData}
+}
diff --git a/man/ParamTMoE-class.Rd b/man/ParamTMoE-class.Rd
new file mode 100644
index 0000000..335e665
--- /dev/null
+++ b/man/ParamTMoE-class.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ParamTMoE.R
+\docType{class}
+\name{ParamTMoE-class}
+\alias{ParamTMoE-class}
+\alias{ParamTMoE}
+\title{A Reference Class which contains parameters of a TMoE model.}
+\description{
+ParamTMoE contains all the parameters of a TMoE model.
+}
+\section{Fields}{
+
+\describe{
+\item{\code{fData}}{\link{FData} object representing the sample.}
+
+\item{\code{K}}{The number of mixture components.}
+
+\item{\code{p}}{The order of the polynomial regression.}
+
+\item{\code{q}}{The dimension of the logistic regression. For the purpose of
+segmentation, it must be set to 1.}
+
+\item{\code{nu}}{degree of freedom}
+
+\item{\code{alpha}}{is the parameter vector of the logistic model with \eqn{alpha_K} being the null vector.}
+
+\item{\code{beta}}{is the vector of regression coefficients of component k,
+the updates for each of the expert component parameters consist in analytically solving a weighted
+Gaussian linear regression problem.}
+
+\item{\code{sigma}}{The variances for the \emph{K} mixture component.}
+
+\item{\code{delta}}{the skewness parameter lambda (by equivalence delta)}
+}}
+
+
+\seealso{
+\link{FData}
+}
diff --git a/man/StatSNMoE-class.Rd b/man/StatSNMoE-class.Rd
new file mode 100644
index 0000000..8c57d82
--- /dev/null
+++ b/man/StatSNMoE-class.Rd
@@ -0,0 +1,95 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/StatSNMoE.R
+\docType{class}
+\name{StatSNMoE-class}
+\alias{StatSNMoE-class}
+\alias{StatSNMoE}
+\title{A Reference Class which contains statistics of a SNMoE model.}
+\description{
+StatMRHLP contains all the parameters of a \link[=ParamSNMoE]{SNMoE} model.
+}
+\section{Fields}{
+
+\describe{
+\item{\code{piik}}{Matrix of size \eqn{(n, K)} representing the probabilities
+\eqn{P(zi = k; W) = P(z_{ik} = 1; W)}{P(zi = k; W) = P(z_ik = 1; W)} of the
+latent variable \eqn{zi,\ i = 1,\dots,m}{zi, i = 1,\dots,n}.}
+
+\item{\code{z_ik}}{Hard segmentation logical matrix of dimension \eqn{(n, K)}
+obtained by the Maximum a posteriori (MAP) rule:
+\eqn{z_{ik} = 1 \ \textrm{if} \ z_{ik} = \textrm{arg} \ \textrm{max}_{k} \
+P(z_i = k | Y, W, \beta);\ 0 \ \textrm{otherwise}}{z_ik = 1 if z_ik =
+arg max_k P(z_i = k | Y, W, \beta); 0 otherwise}, \eqn{k = 1,\dots,K}.}
+
+\item{\code{klas}}{Column matrix of the labels issued from \code{z_ik}. Its elements are
+\eqn{klas(i) = k}, \eqn{k = 1,\dots,K}.}
+
+\item{\code{Ey_k}}{Matrix of dimension \emph{(n,K)}.}
+
+\item{\code{Ey}}{Column matrix of dimension \emph{n}.}
+
+\item{\code{Var_yk}}{Column matrix of dimension \emph{K}.}
+
+\item{\code{Vary}}{Column matrix of dimension \emph{n}.}
+
+\item{\code{log\_lik}}{Numeric. Log-likelihood of the SNMoE model.}
+
+\item{\code{com_loglik}}{Numeric. Complete log-likelihood of the SNMoE model.}
+
+\item{\code{stored_loglik}}{List. Stored values of the log-likelihood at each EM
+iteration.}
+
+\item{\code{BIC}}{Numeric. Value of the BIC (Bayesian Information Criterion)
+criterion. The formula is \eqn{BIC = log\_lik - nu \times
+\textrm{log}(n) / 2}{BIC = log\_lik - nu x log(n) / 2} with \emph{nu} the
+degree of freedom of the SNMoE model.}
+
+\item{\code{ICL}}{Numeric. Value of the ICL (Integrated Completed Likelihood)
+criterion. The formula is \eqn{ICL = com\_loglik - nu \times
+\textrm{log}(n) / 2}{ICL = com_loglik - nu x log(n) / 2} with \emph{nu} the
+degree of freedom of the SNMoE model.}
+
+\item{\code{AIC}}{Numeric. Value of the AIC (Akaike Information Criterion)
+criterion. The formula is \eqn{AIC = log\_lik - nu}{AIC = log\_lik - nu}.}
+
+\item{\code{log_piik_fik}}{Matrix of size \eqn{(n, K)} giving the values of the
+logarithm of the joint probability
+\eqn{P(Y_{i}, \ zi = k)}{P(Yi, zi = k)}, \eqn{i = 1,\dots,n}.}
+
+\item{\code{log_sum_piik_fik}}{Column matrix of size \emph{n} giving the values of
+\eqn{\sum_{k = 1}^{K} \textrm{log} P(Y_{i}, \ zi = k)}{\sum_{k = 1}^{K} log
+P(Yi, zi = k)}, \eqn{i = 1,\dots,n}.}
+
+\item{\code{tik}}{Matrix of size \eqn{(n, K)} giving the posterior probability that
+\eqn{Y_{i}}{Yi} originates from the \eqn{k}-th regression model
+\eqn{P(zi = k | Y, W, \beta)}.}
+
+\item{\code{E1ik}}{To define.}
+
+\item{\code{E2ik}}{To define.}
+}}
+
+\section{Methods}{
+
+\describe{
+\item{\code{EStep(paramSNMoE)}}{Method used in the EM algorithm to update statistics based on parameters
+provided by \code{paramSNMoE} (prior and posterior probabilities).}
+
+\item{\code{MAP()}}{calcule une partition d'un echantillon par la regle du Maximum A Posteriori ?? partir des probabilites a posteriori
+Entrees : post_probas , Matrice de dimensions [n x K] des probabibiltes a posteriori (matrice de la partition floue)
+n : taille de l'echantillon
+K : nombres de classes
+klas(i) = arg   max (post_probas(i,k)) , for all i=1,...,n
+1<=k<=K
+= arg   max  p(zi=k|xi;theta)
+1<=k<=K
+= arg   max  p(zi=k;theta)p(xi|zi=k;theta)/sum{l=1}^{K}p(zi=l;theta) p(xi|zi=l;theta)
+1<=k<=K
+Sorties : classes : vecteur collones contenant les classe (1:K)
+Z : Matrice de dimension [nxK] de la partition dure : ses elements sont zik, avec zik=1 si xi
+appartient ?? la classe k (au sens du MAP) et zero sinon.}
+}}
+
+\seealso{
+\link{ParamSNMoE}, \link{FData}
+}
diff --git a/man/StatStMoE-class.Rd b/man/StatStMoE-class.Rd
new file mode 100644
index 0000000..2f1e1a3
--- /dev/null
+++ b/man/StatStMoE-class.Rd
@@ -0,0 +1,103 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/StatStMoE.R
+\docType{class}
+\name{StatStMoE-class}
+\alias{StatStMoE-class}
+\alias{StatStMoE}
+\title{A Reference Class which contains statistics of a StMoE model.}
+\description{
+StatMRHLP contains all the parameters of a \link[=ParamStMoE]{StMoE} model.
+}
+\section{Fields}{
+
+\describe{
+\item{\code{piik}}{Matrix of size \eqn{(n, K)} representing the probabilities
+\eqn{P(zi = k; W) = P(z_{ik} = 1; W)}{P(zi = k; W) = P(z_ik = 1; W)} of the
+latent variable \eqn{zi,\ i = 1,\dots,m}{zi, i = 1,\dots,n}.}
+
+\item{\code{z_ik}}{Hard segmentation logical matrix of dimension \eqn{(n, K)}
+obtained by the Maximum a posteriori (MAP) rule:
+\eqn{z_{ik} = 1 \ \textrm{if} \ z_{ik} = \textrm{arg} \ \textrm{max}_{k} \
+P(z_i = k | Y, W, \beta);\ 0 \ \textrm{otherwise}}{z_ik = 1 if z_ik =
+arg max_k P(z_i = k | Y, W, \beta); 0 otherwise}, \eqn{k = 1,\dots,K}.}
+
+\item{\code{klas}}{Column matrix of the labels issued from \code{z_ik}. Its elements are
+\eqn{klas(i) = k}, \eqn{k = 1,\dots,K}.}
+
+\item{\code{Ey_k}}{Matrix of dimension \emph{(n,K)}.}
+
+\item{\code{Ey}}{Column matrix of dimension \emph{n}.}
+
+\item{\code{Var_yk}}{Column matrix of dimension \emph{K}.}
+
+\item{\code{Var_y}}{Column matrix of dimension \emph{n}.}
+
+\item{\code{log\_lik}}{Numeric. Log-likelihood of the StMoE model.}
+
+\item{\code{com_loglik}}{Numeric. Complete log-likelihood of the StMoE model.}
+
+\item{\code{stored_loglik}}{List. Stored values of the log-likelihood at each EM
+iteration.}
+
+\item{\code{BIC}}{Numeric. Value of the BIC (Bayesian Information Criterion)
+criterion. The formula is \eqn{BIC = log\_lik - nu \times
+\textrm{log}(n) / 2}{BIC = log\_lik - nu x log(n) / 2} with \emph{nu} the
+degree of freedom of the StMoE model.}
+
+\item{\code{ICL}}{Numeric. Value of the ICL (Integrated Completed Likelihood)
+criterion. The formula is \eqn{ICL = com\_loglik - nu \times
+\textrm{log}(n) / 2}{ICL = com_loglik - nu x log(n) / 2} with \emph{nu} the
+degree of freedom of the StMoE model.}
+
+\item{\code{AIC}}{Numeric. Value of the AIC (Akaike Information Criterion)
+criterion. The formula is \eqn{AIC = log\_lik - nu}{AIC = log\_lik - nu}.}
+
+\item{\code{log_piik_fik}}{Matrix of size \eqn{(n, K)} giving the values of the
+logarithm of the joint probability
+\eqn{P(Y_{i}, \ zi = k)}{P(Yi, zi = k)}, \eqn{i = 1,\dots,n}.}
+
+\item{\code{log_sum_piik_fik}}{Column matrix of size \emph{n} giving the values of
+\eqn{\sum_{k = 1}^{K} \textrm{log} P(Y_{i}, \ zi = k)}{\sum_{k = 1}^{K} log
+P(Yi, zi = k)}, \eqn{i = 1,\dots,n}.}
+
+\item{\code{tik}}{Matrix of size \eqn{(n, K)} giving the posterior probability that
+\eqn{Y_{i}}{Yi} originates from the \eqn{k}-th regression model
+\eqn{P(zi = k | Y, W, \beta)}.}
+
+\item{\code{wik}}{To define.}
+
+\item{\code{dik}}{To define.}
+
+\item{\code{stme_pdf}}{skew-t mixture of experts density}
+
+\item{\code{E1ik}}{To define.}
+
+\item{\code{E2ik}}{To define.}
+
+\item{\code{E3ik}}{To define.}
+}}
+
+\section{Methods}{
+
+\describe{
+\item{\code{EStep(paramStMoE)}}{Method used in the EM algorithm to update statistics based on parameters
+provided by \code{paramStMoE} (prior and posterior probabilities).}
+
+\item{\code{MAP()}}{calcule une partition d'un echantillon par la regle du Maximum A Posteriori ?? partir des probabilites a posteriori
+Entrees : post_probas , Matrice de dimensions [n x K] des probabibiltes a posteriori (matrice de la partition floue)
+n : taille de l'echantillon
+K : nombres de classes
+klas(i) = arg   max (post_probas(i,k)) , for all i=1,...,n
+1<=k<=K
+= arg   max  p(zi=k|xi;theta)
+1<=k<=K
+= arg   max  p(zi=k;theta)p(xi|zi=k;theta)/sum{l=1}^{K}p(zi=l;theta) p(xi|zi=l;theta)
+1<=k<=K
+Sorties : classes : vecteur collones contenant les classe (1:K)
+Z : Matrice de dimension [nxK] de la partition dure : ses elements sont zik, avec zik=1 si xi
+appartient et la classe k (au sens du MAP) et zero sinon.}
+}}
+
+\seealso{
+\link{ParamStMoE}, \link{FData}
+}
diff --git a/man/StatTMoE-class.Rd b/man/StatTMoE-class.Rd
new file mode 100644
index 0000000..27d44c0
--- /dev/null
+++ b/man/StatTMoE-class.Rd
@@ -0,0 +1,90 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/StatTMoE.R
+\docType{class}
+\name{StatTMoE-class}
+\alias{StatTMoE-class}
+\alias{StatTMoE}
+\title{A Reference Class which contains statistics of a TMoE model.}
+\description{
+StatTMoE contains all the parameters of a \link[=ParamTMoE]{TMoE} model.
+}
+\section{Fields}{
+
+\describe{
+\item{\code{piik}}{Matrix of size \eqn{(n, K)} representing the probabilities
+\eqn{P(zi = k; W) = P(z_{ik} = 1; W)}{P(zi = k; W) = P(z_ik = 1; W)} of the
+latent variable \eqn{zi,\ i = 1,\dots,m}{zi, i = 1,\dots,n}.}
+
+\item{\code{z_ik}}{Hard segmentation logical matrix of dimension \eqn{(n, K)}
+obtained by the Maximum a posteriori (MAP) rule:
+\eqn{z_{ik} = 1 \ \textrm{if} \ z_{ik} = \textrm{arg} \ \textrm{max}_{k} \
+P(z_i = k | Y, W, \beta);\ 0 \ \textrm{otherwise}}{z_ik = 1 if z_ik =
+arg max_k P(z_i = k | Y, W, \beta); 0 otherwise}, \eqn{k = 1,\dots,K}.}
+
+\item{\code{klas}}{Column matrix of the labels issued from \code{z_ik}. Its elements are
+\eqn{klas(i) = k}, \eqn{k = 1,\dots,K}.}
+
+\item{\code{Wik}}{Matrix of dimension \emph{(nm,K)}.}
+
+\item{\code{Ey_k}}{Matrix of dimension \emph{(n,K)}.}
+
+\item{\code{Ey}}{Column matrix of dimension \emph{n}.}
+
+\item{\code{Var_yk}}{Column matrix of dimension \emph{K}.}
+
+\item{\code{Vary}}{Column matrix of dimension \emph{n}.}
+
+\item{\code{log\_lik}}{Numeric. Log-likelihood of the StMoE model.}
+
+\item{\code{com_loglik}}{Numeric. Complete log-likelihood of the StMoE model.}
+
+\item{\code{stored_loglik}}{List. Stored values of the log-likelihood at each EM
+iteration.}
+
+\item{\code{BIC}}{Numeric. Value of the BIC (Bayesian Information Criterion)
+criterion. The formula is \eqn{BIC = log\_lik - nu \times
+\textrm{log}(n) / 2}{BIC = log\_lik - nu x log(n) / 2} with \emph{nu} the
+degree of freedom of the StMoE model.}
+
+\item{\code{ICL}}{Numeric. Value of the ICL (Integrated Completed Likelihood)
+criterion. The formula is \eqn{ICL = com\_loglik - nu \times
+\textrm{log}(n) / 2}{ICL = com_loglik - nu x log(n) / 2} with \emph{nu} the
+degree of freedom of the StMoE model.}
+
+\item{\code{AIC}}{Numeric. Value of the AIC (Akaike Information Criterion)
+criterion. The formula is \eqn{AIC = log\_lik - nu}{AIC = log\_lik - nu}.}
+
+\item{\code{log_piik_fik}}{Matrix of size \eqn{(n, K)} giving the values of the
+logarithm of the joint probability
+\eqn{P(Y_{i}, \ zi = k)}{P(Yi, zi = k)}, \eqn{i = 1,\dots,n}.}
+
+\item{\code{log_sum_piik_fik}}{Column matrix of size \emph{n} giving the values of
+\eqn{\sum_{k = 1}^{K} \textrm{log} P(Y_{i}, \ zi = k)}{\sum_{k = 1}^{K} log
+P(Yi, zi = k)}, \eqn{i = 1,\dots,n}.}
+
+\item{\code{tik}}{Matrix of size \eqn{(n, K)} giving the posterior probability that
+\eqn{Y_{i}}{Yi} originates from the \eqn{k}-th regression model
+\eqn{P(zi = k | Y, W, \beta)}.}
+}}
+
+\section{Methods}{
+
+\describe{
+\item{\code{MAP()}}{calcule une partition d'un echantillon par la regle du Maximum A Posteriori ?? partir des probabilites a posteriori
+Entrees : post_probas , Matrice de dimensions [n x K] des probabibiltes a posteriori (matrice de la partition floue)
+n : taille de l'echantillon
+K : nombres de classes
+klas(i) = arg   max (post_probas(i,k)) , for all i=1,...,n
+1<=k<=K
+= arg   max  p(zi=k|xi;theta)
+1<=k<=K
+= arg   max  p(zi=k;theta)p(xi|zi=k;theta)/sum{l=1}^{K}p(zi=l;theta) p(xi|zi=l;theta)
+1<=k<=K
+Sorties : classes : vecteur collones contenant les classe (1:K)
+Z : Matrice de dimension [nxK] de la partition dure : ses elements sont zik, avec zik=1 si xi
+appartient ?? la classe k (au sens du MAP) et zero sinon.}
+}}
+
+\seealso{
+\link{ParamStMoE}, \link{FData}
+}
diff --git a/man/emSNMoE.Rd b/man/emSNMoE.Rd
new file mode 100644
index 0000000..c190bb3
--- /dev/null
+++ b/man/emSNMoE.Rd
@@ -0,0 +1,62 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/emSNMoE.R
+\name{emSNMoE}
+\alias{emSNMoE}
+\title{emSNMoE is used to fit a SNMoE model.}
+\usage{
+emSNMoE(X, Y, K, p = 3, q = 1, n_tries = 1, max_iter = 1500,
+  threshold = 1e-06, verbose = FALSE, verbose_IRLS = FALSE)
+}
+\arguments{
+\item{X}{Numeric vector of length \emph{m} representing the covariates.}
+
+\item{Y}{Matrix of size \eqn{(n, m)} representing \emph{n} functions of \code{X}
+observed at points \eqn{1,\dots,m}.}
+
+\item{K}{The number of mixture components.}
+
+\item{p}{The order of the polynomial regression.}
+
+\item{q}{The dimension of the logistic regression. For the purpose of
+segmentation, it must be set to 1.}
+
+\item{n_tries}{Number of times EM algorithm will be launched.
+The solution providing the highest log-likelihood will be returned.
+
+If \code{n_tries} > 1, then for the first pass, parameters are initialized
+by uniformly segmenting the data into K segments, and for the next passes,
+parameters are initialized by randomly segmenting the data into K contiguous
+segments.}
+
+\item{max_iter}{The maximum number of iterations for the EM algorithm.}
+
+\item{threshold}{A numeric value specifying the threshold for the relative
+difference of log-likelihood between two steps  of the EM as stopping
+criteria.}
+
+\item{verbose}{A logical value indicating whether values of the
+log-likelihood should be printed during EM iterations.}
+
+\item{verbose_IRLS}{A logical value indicating whether values of the
+criterion optimized by IRLS should be printed at each step of the EM
+algorithm.}
+}
+\value{
+EM returns an object of class \link{ModelSNMoE}.
+}
+\description{
+emSNMoE is used to fit a SNMoE model. The estimation method is performed by
+the Expectation-Maximization algorithm.
+}
+\details{
+emSNMoE function is based on the EM algorithm. This functions starts
+with an initialization of the parameters done by the method \code{initParam} of
+the class \link{ParamSNMoE}, then it alternates between a E-Step
+(method of the class \link{StatSNMoE}) and a M-Step (method of the class
+\link{ParamSNMoE}) until convergence (until the absolute difference of
+log-likelihood between two steps of the EM algorithm is less than the
+\code{threshold} parameter).
+}
+\seealso{
+\link{ModelSNMoE}, \link{ParamSNMoE}, \link{StatSNMoE}
+}
diff --git a/man/emStMoE.Rd b/man/emStMoE.Rd
new file mode 100644
index 0000000..dea6c07
--- /dev/null
+++ b/man/emStMoE.Rd
@@ -0,0 +1,65 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/emStMoE.R
+\name{emStMoE}
+\alias{emStMoE}
+\title{emStMoE is used to fit a StMoE model.}
+\usage{
+emStMoE(X, Y, K, p = 3, q = 1, n_tries = 1, max_iter = 1500,
+  threshold = 1e-06, verbose = FALSE, verbose_IRLS = FALSE)
+}
+\arguments{
+\item{X}{Numeric vector of length \emph{m} representing the covariates.}
+
+\item{Y}{Matrix of size \eqn{(n, m)} representing \emph{n} functions of \code{X}
+observed at points \eqn{1,\dots,m}.}
+
+\item{K}{The number of mixture components.}
+
+\item{p}{The order of the polynomial regression.}
+
+\item{q}{The dimension of the logistic regression. For the purpose of
+segmentation, it must be set to 1.}
+
+\item{n_tries}{Number of times EM algorithm will be launched.
+The solution providing the highest log-likelihood will be returned.
+
+If \code{n_tries} > 1, then for the first pass, parameters are initialized
+by uniformly segmenting the data into K segments, and for the next passes,
+parameters are initialized by randomly segmenting the data into K contiguous
+segments.}
+
+\item{max_iter}{The maximum number of iterations for the EM algorithm.}
+
+\item{threshold}{A numeric value specifying the threshold for the relative
+difference of log-likelihood between two steps  of the EM as stopping
+criteria.}
+
+\item{verbose}{A logical value indicating whether values of the
+log-likelihood should be printed during EM iterations.}
+
+\item{verbose_IRLS}{A logical value indicating whether values of the
+criterion optimized by IRLS should be printed at each step of the EM
+algorithm.}
+
+\item{variance_type}{Numeric indicating if the model is homoskedastic
+(\code{variance_type} = 1) or heteroskedastic (\code{variance_type} = 2).}
+}
+\value{
+EM returns an object of class \link{ModelStMoE}.
+}
+\description{
+emStMoE is used to fit a StMoE model. The estimation method is performed by
+the Expectation-Maximization algorithm.
+}
+\details{
+emStMoE function is based on the EM algorithm. This functions starts
+with an initialization of the parameters done by the method \code{initParam} of
+the class \link{ParamStMoE}, then it alternates between a E-Step
+(method of the class \link{StatStMoE}) and a M-Step (method of the class
+\link{ParamStMoE}) until convergence (until the absolute difference of
+log-likelihood between two steps of the EM algorithm is less than the
+\code{threshold} parameter).
+}
+\seealso{
+\link{ModelStMoE}, \link{ParamStMoE}, \link{StatStMoE}
+}
diff --git a/man/emTMoE.Rd b/man/emTMoE.Rd
new file mode 100644
index 0000000..f78fbbb
--- /dev/null
+++ b/man/emTMoE.Rd
@@ -0,0 +1,65 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/emTMoE.R
+\name{emTMoE}
+\alias{emTMoE}
+\title{emTMoE is used to fit a TMoE model.}
+\usage{
+emTMoE(X, Y, K, p = 3, q = 1, n_tries = 1, max_iter = 1500,
+  threshold = 1e-06, verbose = FALSE, verbose_IRLS = FALSE)
+}
+\arguments{
+\item{X}{Numeric vector of length \emph{m} representing the covariates.}
+
+\item{Y}{Matrix of size \eqn{(n, m)} representing \emph{n} functions of \code{X}
+observed at points \eqn{1,\dots,m}.}
+
+\item{K}{The number of mixture components.}
+
+\item{p}{The order of the polynomial regression.}
+
+\item{q}{The dimension of the logistic regression. For the purpose of
+segmentation, it must be set to 1.}
+
+\item{n_tries}{Number of times EM algorithm will be launched.
+The solution providing the highest log-likelihood will be returned.
+
+If \code{n_tries} > 1, then for the first pass, parameters are initialized
+by uniformly segmenting the data into K segments, and for the next passes,
+parameters are initialized by randomly segmenting the data into K contiguous
+segments.}
+
+\item{max_iter}{The maximum number of iterations for the EM algorithm.}
+
+\item{threshold}{A numeric value specifying the threshold for the relative
+difference of log-likelihood between two steps  of the EM as stopping
+criteria.}
+
+\item{verbose}{A logical value indicating whether values of the
+log-likelihood should be printed during EM iterations.}
+
+\item{verbose_IRLS}{A logical value indicating whether values of the
+criterion optimized by IRLS should be printed at each step of the EM
+algorithm.}
+
+\item{variance_type}{Numeric indicating if the model is homoskedastic
+(\code{variance_type} = 1) or heteroskedastic (\code{variance_type} = 2).}
+}
+\value{
+EM returns an object of class \link{ModelTMoE}.
+}
+\description{
+emTMoE is used to fit a TMoE model. The estimation method is performed by
+the Expectation-Maximization algorithm.
+}
+\details{
+emStMoE function is based on the EM algorithm. This functions starts
+with an initialization of the parameters done by the method \code{initParam} of
+the class \link{ParamTMoE}, then it alternates between a E-Step
+(method of the class \link{StatTMoE}) and a M-Step (method of the class
+\link{ParamTMoE}) until convergence (until the absolute difference of
+log-likelihood between two steps of the EM algorithm is less than the
+\code{threshold} parameter).
+}
+\seealso{
+\link{ModelTMoE}, \link{ParamTMoE}, \link{StatTMoE}
+}
diff --git a/man/figures/README-unnamed-chunk-5-1.png b/man/figures/README-unnamed-chunk-5-1.png
new file mode 100644
index 0000000..0289968
Binary files /dev/null and b/man/figures/README-unnamed-chunk-5-1.png differ
diff --git a/man/figures/README-unnamed-chunk-5-2.png b/man/figures/README-unnamed-chunk-5-2.png
new file mode 100644
index 0000000..5b319be
Binary files /dev/null and b/man/figures/README-unnamed-chunk-5-2.png differ
diff --git a/man/figures/README-unnamed-chunk-5-3.png b/man/figures/README-unnamed-chunk-5-3.png
new file mode 100644
index 0000000..ff97775
Binary files /dev/null and b/man/figures/README-unnamed-chunk-5-3.png differ
diff --git a/man/figures/README-unnamed-chunk-5-4.png b/man/figures/README-unnamed-chunk-5-4.png
new file mode 100644
index 0000000..9d65284
Binary files /dev/null and b/man/figures/README-unnamed-chunk-5-4.png differ
diff --git a/man/figures/README-unnamed-chunk-5-5.png b/man/figures/README-unnamed-chunk-5-5.png
new file mode 100644
index 0000000..ae8b8b3
Binary files /dev/null and b/man/figures/README-unnamed-chunk-5-5.png differ
diff --git a/meteorit.Rproj b/meteorit.Rproj
new file mode 100644
index 0000000..0ef8ced
--- /dev/null
+++ b/meteorit.Rproj
@@ -0,0 +1,20 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: ASCII
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+AutoAppendNewline: Yes
+StripTrailingWhitespace: Yes
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source
diff --git a/src/.gitignore b/src/.gitignore
new file mode 100644
index 0000000..22034c4
--- /dev/null
+++ b/src/.gitignore
@@ -0,0 +1,3 @@
+*.o
+*.so
+*.dll
diff --git a/src/IRLS.cpp b/src/IRLS.cpp
new file mode 100644
index 0000000..78d236d
--- /dev/null
+++ b/src/IRLS.cpp
@@ -0,0 +1,284 @@
+// [[Rcpp::depends(RcppArmadillo)]]
+
+#include <RcppArmadillo.h>
+#include "multinomialLogit.h"
+
+using namespace Rcpp;
+
+// This is a simple example of exporting a C++ function to R. You can
+// source this function into an R session using the Rcpp::sourceCpp
+// function (or via the Source button on the editor toolbar). Learn
+// more about Rcpp at:
+//
+//   http://www.rcpp.org/
+//   http://adv-r.had.co.nz/Rcpp.html
+//   http://gallery.rcpp.org/
+//
+
+// [[Rcpp::export]]
+List IRLS(arma::mat& X, arma::mat& Tau, arma::mat& Gamma, arma::mat& Winit, bool verbose = false) {
+
+  // res = IRLS_MixFRHLP(X, Tau, Gamma, Winit, verbose) : an efficient Iteratively Reweighted Least-Squares (IRLS) algorithm for estimating
+  // the parameters of a multinomial logistic regression model given the
+  // "predictors" X and a partition (hard or smooth) Tau into K>=2 segments,
+  // and a cluster weights Gamma (hard or smooth)
+  //
+  //
+  // Inputs :
+  //
+  //         X : desgin matrix for the logistic weights.  dim(X) = [nx(q+1)]
+  //                            X = [1 t1 t1^2 ... t1^q
+  //                                 1 t2 t2^2 ... t2^q
+  //                                      ..
+  //                                 1 ti ti^2 ... ti^q
+  //                                      ..
+  //                                 1 tn tn^2 ... tn^q]
+  //            q being the number of predictors
+  //         Tau : matrix of a hard or fauzzy partition of the data (here for
+  //         the RHLP model, Tau is the fuzzy partition represented by the
+  //         posterior probabilities (responsibilities) (tik) obtained at the E-Step).
+  //
+  //         Winit : initial parameter values W(0). dim(Winit) = [(q+1)x(K-1)]
+  //         verbose : 1 to print the loglikelihood values during the IRLS
+  //         iterations, 0 if not
+  //
+  // Outputs :
+  //
+  //          res : structure containing the fields:
+  //              W : the estimated parameter vector. matrix of dim [(q+1)x(K-1)]
+  //                  (the last vector being the null vector)
+  //              piigk : the logistic probabilities (dim [n x K])
+  //              loglik : the value of the maximized objective
+  //              LL : stored values of the maximized objective during the
+  //              IRLS training
+  //
+  //        Probs(i,gk) = Pro(segment k|cluster g;W)
+  //                    = \pi_{ik}(W)
+  //                           exp(wgk'vi)
+  //                    =  ---------------------------
+  //                      1+sum_{l=1}^{K-1} exp(wgl'vi)
+  //
+  //       with :
+  //            * Probs(i,gk) is the prob of component k at time t_i in
+  //            cluster g
+  //            i=1,...n,j=1...m,  k=1,...,K,
+  //            * vi = [1,ti,ti^2,...,ti^q]^T;
+  //       The parameter vecrots are in the matrix W=[w1,...,wK] (with wK is the null vector);
+  //
+
+  //// References
+  // Please cite the following papers for this code:
+  //
+  //
+  // @INPROCEEDINGS{Chamroukhi-IJCNN-2009,
+  //   AUTHOR =       {Chamroukhi, F. and Sam\'e,  A. and Govaert, G. and Aknin, P.},
+  //   TITLE =        {A regression model with a hidden logistic process for feature extraction from time series},
+  //   BOOKTITLE =    {International Joint Conference on Neural Networks (IJCNN)},
+  //   YEAR =         {2009},
+  //   month = {June},
+  //   pages = {489--496},
+  //   Address = {Atlanta, GA},
+  //  url = {https://chamroukhi.users.lmno.cnrs.fr/papers/chamroukhi_ijcnn2009.pdf}
+  // }
+  //
+  // @article{chamroukhi_et_al_NN2009,
+  // 	Address = {Oxford, UK, UK},
+  // 	Author = {Chamroukhi, F. and Sam\'{e}, A. and Govaert, G. and Aknin, P.},
+  // 	Date-Added = {2014-10-22 20:08:41 +0000},
+  // 	Date-Modified = {2014-10-22 20:08:41 +0000},
+  // 	Journal = {Neural Networks},
+  // 	Number = {5-6},
+  // 	Pages = {593--602},
+  // 	Publisher = {Elsevier Science Ltd.},
+  // 	Title = {Time series modeling by a regression approach based on a latent process},
+  // 	Volume = {22},
+  // 	Year = {2009},
+  // 	url  = {https://chamroukhi.users.lmno.cnrs.fr/papers/Chamroukhi_Neural_Networks_2009.pdf}
+  // 	}
+  // @article{Chamroukhi-FDA-2018,
+  // 	Journal = {},
+  // 	Author = {Faicel Chamroukhi and Hien D. Nguyen},
+  // 	Volume = {},
+  // 	Title = {Model-Based Clustering and Classification of Functional Data},
+  // 	Year = {2018},
+  // 	eprint ={arXiv:1803.00276v2},
+  // 	url =  {https://chamroukhi.users.lmno.cnrs.fr/papers/MBCC-FDA.pdf}
+  // 	}
+  //
+  //
+  //////////////////////////////////////////////////////////////////////////////////////////////
+
+  int n = Tau.n_rows;
+  int K = Tau.n_cols;
+
+  // Checker nrow(X) = nrow(Tau) = n
+
+  int q = X.n_cols; // q here is (q+1)
+
+  /* Handle NULL Winit and Gamma
+   *
+   * if nargin<4; Winit = zeros(q,K-1);end % if there is no a specified initialization
+   * if nargin<3; Gamma = ones(n,1);end % for standard weighted multinomial logistic regression
+   *
+   */
+
+
+
+  double lambda = 1e-9; // if a MAP regularization (a gaussian prior on W) (L_2 penalization); lambda is a positive hyperparameter
+  arma::mat I(q * (K - 1), q * (K - 1), arma::fill::eye); // Define an identity matrix
+
+  arma::mat piik_old(n, K, arma::fill::zeros);
+  arma::mat piik(n, K, arma::fill::zeros);
+  arma::colvec gwk(n, arma::fill::zeros);
+  arma::colvec vq(n, arma::fill::zeros);
+
+  arma::mat Hkl(q, q, arma::fill::zeros);
+  arma::colvec vqa(n, arma::fill::zeros);
+  arma::colvec vqb(n, arma::fill::zeros);
+  double hwk;
+  arma::colvec w(q * (K - 1), arma::fill::zeros);
+
+
+  // IRLS Initialization (iter = 0)
+  arma::mat W_old = Winit;
+
+  List out = multinomialLogit(W_old, X, Tau, Gamma);
+
+  double loglik_old = out["loglik"];
+  arma::mat tmp = out["piik"];
+  piik_old = tmp;
+  piik = tmp;
+
+
+  loglik_old = loglik_old - lambda * std::pow(arma::norm(W_old), 2.0);
+
+  int iter = 0;
+  int max_iter = 300;
+  double loglik = 0;
+
+  // List ret;
+  List LL;
+
+  bool converge1 = false;
+  bool converge2 = false;
+  bool converge = false;
+
+  arma::mat W(q, K - 1, arma::fill::zeros);
+
+  if (verbose == true) {
+    Rcout << "IRLS : Iteration : "  + std::to_string(iter) + "; Log-lik : " + std::to_string(loglik_old) + "\n";
+  }
+
+  // IRLS
+
+  // Hw_old a squared matrix of dimensions  q*(K-1) x  q*(K-1)
+  int hx = q * (K - 1);
+  arma::mat Hw_old(hx, hx, arma::fill::zeros);
+  arma::mat gw_old(q, K - 1, arma::fill::zeros);
+
+
+  while (not(converge) && (iter < max_iter)) {
+
+    // Gradient :
+    gw_old.reshape(q, (K - 1));
+    for (int k = 0; k < K - 1; k++) {
+      gwk = Gamma % (Tau.col(k) - piik_old.col(k));
+
+      for (int qq = 0; qq < q; qq++) {
+        vq = X.col(qq);
+        gw_old(qq, k) = arma::as_scalar(gwk.t() * vq);
+      }
+    }
+    gw_old.reshape(q * (K - 1), 1);
+
+    // Hessian
+    for (int k = 0; k < K - 1; k++) {
+      for (int ell = 0; ell < K - 1; ell++) {
+        bool delta_kl = (k == ell); // kronecker delta
+        gwk = Gamma % (piik_old.col(k) % (arma::mat(n, 1, arma::fill::ones) * delta_kl - piik_old.col(ell)));
+        for (int qqa = 0; qqa < q; qqa++) {
+          vqa = X.col(qqa);
+          for (int qqb = 0; qqb < q; qqb++) {
+            vqb = X.col(qqb);
+            hwk = as_scalar(vqb.t() * (gwk % vqa));
+            Hkl(qqa, qqb) = hwk;
+
+          }
+        }
+        Hw_old.submat(k * q, ell * q, (k + 1) * q - 1, (ell + 1) * q - 1) = -Hkl;
+      }
+    }
+
+    // If a gaussian prior on W (lambda ~=0)
+    Hw_old = Hw_old + lambda * I;
+    gw_old = gw_old - lambda * vectorise(W_old, 0);
+
+    // Newton Raphson : W(t+1) = W(t) - H(W(t))^(-1)g(W(t))
+    w = vectorise(W_old, 0) - arma::inv(Hw_old) * gw_old; // [(q+1)x(K-1),1]
+    W.reshape(1, q * (K - 1));
+    W = w;
+    W.reshape(q, K - 1); // [(q+1)*(K-1)]
+
+    // Mise a jour des probas et de la loglik
+    out = multinomialLogit(W, X, Tau, Gamma);
+    loglik = out["loglik"];
+    loglik = loglik - lambda * std::pow(arma::norm(W), 2.0);
+    arma::mat tmp = out["piik"];
+    piik = tmp;
+
+    // Check if Qw1(w^(t+1),w^(t))> Qw1(w^(t),w^(t))
+    // (adaptive stepsize in case of troubles with stepsize 1) Newton Raphson : W(t+1) = W(t) - stepsize*H(W)^(-1)*g(W)
+    double stepsize = 1; // Initialisation pas d'adaptation de l'algo Newton raphson
+    double alpha = 2;
+    while (loglik < loglik_old) {
+
+      stepsize = stepsize / alpha;
+
+      // Recalculate the parameter W and the "loglik"
+      w = vectorise(W_old, 0) - stepsize * arma::inv(Hw_old) * gw_old;
+      W.reshape(1, q * (K - 1));
+      W = w;
+      W.reshape(q, K - 1);
+      out = multinomialLogit(W, X, Tau, Gamma);
+      loglik = out["loglik"];
+
+      loglik = loglik - lambda * std::pow(arma::norm(W), 2.0);
+
+      arma::mat tmp = out["piik"];
+      piik = tmp;
+    }
+
+    converge1 = std::abs((loglik - loglik_old) / loglik_old) <= 1e-7;
+    converge2 = std::abs(loglik - loglik_old) <= 1e-6;
+    converge = converge1 || converge2;
+
+    piik_old = piik;
+    W_old = W;
+    iter = iter + 1;
+    LL[std::to_string(iter)] = loglik_old;
+    loglik_old = loglik;
+
+    if (verbose == true) {
+      Rcout << "IRLS : Iteration : "  + std::to_string(iter) + "; Log-lik" + std::to_string(loglik_old) + "\n";
+    }
+  } // End of IRLS
+
+  if (converge == true) {
+    if (verbose == true) {
+      Rcout << "\n";
+      Rcout << "IRLS convergence OK ; nbr of iteration "  + std::to_string(iter) + "\n";
+      Rcout << "\n";
+    }
+  } else {
+    Rcout << "\n";
+    Rcout << "IRLS : doesn't converged (augment the number of iterations > "  + std::to_string(iter) + "\n";
+  }
+
+  double reg_irls = 0;
+  if (lambda != 0) { // Calculate the value of the regularization part to calculate the value of the MAP criterion in case of regularization
+    reg_irls = - lambda * std::pow(arma::norm(W), 2.0); // bayesian l2 regularization
+  }
+
+  return List::create(Named("W") = W, Named("LL") = LL, Named("loglik") = loglik, Named("piik") = piik, Named("reg_irls") = reg_irls);
+  // return ret;
+}
diff --git a/src/Makevars b/src/Makevars
new file mode 100644
index 0000000..05ee448
--- /dev/null
+++ b/src/Makevars
@@ -0,0 +1,2 @@
+CXX_STD = CXX11
+PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
diff --git a/src/Makevars.win b/src/Makevars.win
new file mode 100644
index 0000000..d715e49
--- /dev/null
+++ b/src/Makevars.win
@@ -0,0 +1,3 @@
+CXX_STD = CXX11
+PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS)
+PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
new file mode 100644
index 0000000..2974da8
--- /dev/null
+++ b/src/RcppExports.cpp
@@ -0,0 +1,48 @@
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
+// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#include <RcppArmadillo.h>
+#include <Rcpp.h>
+
+using namespace Rcpp;
+
+// IRLS
+List IRLS(arma::mat& X, arma::mat& Tau, arma::mat& Gamma, arma::mat& Winit, bool verbose);
+RcppExport SEXP _meteorit_IRLS(SEXP XSEXP, SEXP TauSEXP, SEXP GammaSEXP, SEXP WinitSEXP, SEXP verboseSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type Tau(TauSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type Gamma(GammaSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type Winit(WinitSEXP);
+    Rcpp::traits::input_parameter< bool >::type verbose(verboseSEXP);
+    rcpp_result_gen = Rcpp::wrap(IRLS(X, Tau, Gamma, Winit, verbose));
+    return rcpp_result_gen;
+END_RCPP
+}
+// multinomialLogit
+List multinomialLogit(arma::mat& W, arma::mat& X, arma::mat& Y, arma::mat& Gamma);
+RcppExport SEXP _meteorit_multinomialLogit(SEXP WSEXP, SEXP XSEXP, SEXP YSEXP, SEXP GammaSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::mat& >::type W(WSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type Y(YSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type Gamma(GammaSEXP);
+    rcpp_result_gen = Rcpp::wrap(multinomialLogit(W, X, Y, Gamma));
+    return rcpp_result_gen;
+END_RCPP
+}
+
+static const R_CallMethodDef CallEntries[] = {
+    {"_meteorit_IRLS", (DL_FUNC) &_meteorit_IRLS, 5},
+    {"_meteorit_multinomialLogit", (DL_FUNC) &_meteorit_multinomialLogit, 4},
+    {NULL, NULL, 0}
+};
+
+RcppExport void R_init_meteorit(DllInfo *dll) {
+    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+    R_useDynamicSymbols(dll, FALSE);
+}
diff --git a/src/multinomialLogit.cpp b/src/multinomialLogit.cpp
new file mode 100644
index 0000000..b47a1a2
--- /dev/null
+++ b/src/multinomialLogit.cpp
@@ -0,0 +1,117 @@
+// [[Rcpp::depends(RcppArmadillo)]]
+
+#include <RcppArmadillo.h>
+
+using namespace Rcpp;
+
+// This is a simple example of exporting a C++ function to R. You can
+// source this function into an R session using the Rcpp::sourceCpp
+// function (or via the Source button on the editor toolbar). Learn
+// more about Rcpp at:
+//
+//   http://www.rcpp.org/
+//   http://adv-r.had.co.nz/Rcpp.html
+//   http://gallery.rcpp.org/
+//
+
+// [[Rcpp::export]]
+List multinomialLogit(arma::mat& W, arma::mat& X, arma::mat& Y, arma::mat& Gamma) {
+
+  //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+  // function [probs, loglik] = logit_model_MixRHLP(W, X, Y, Gamma)
+  //
+  // calculates the pobabilities according to multinomial logistic model:
+  //
+  // probs(i,k) = p(zi=k;W)= \pi_{ik}(W)
+  //                                  exp(wk'vi)
+  //                        =  ----------------------------
+  //                          1 + sum_{l=1}^{K-1} exp(wl'vi)
+  // for i=1,...,n and k=1...K
+  //
+  // Inputs :
+  //
+  //         1. W : parametre du modele logistique ,Matrice de dimensions
+  //         [(q+1)x(K-1)]des vecteurs parametre wk. W = [w1 .. wk..w(K-1)]
+  //         avec les wk sont des vecteurs colonnes de dim [(q+1)x1], le dernier
+  //         est suppose nul (sum_{k=1}^K \pi_{ik} = 1 -> \pi{iK} =
+  //         1-sum_{l=1}^{K-1} \pi{il}. vi : vecteur colonne de dimension [(q+1)x1]
+  //         qui est la variable explicative (ici le temps): vi = [1;ti;ti^2;...;ti^q];
+  //         2. M : Matrice de dimensions [nx(q+1)] des variables explicatives.
+  //            M = transpose([v1... vi ....vn])
+  //              = [1 t1 t1^2 ... t1^q
+  //                 1 t2 t2^2 ... t2^q
+  //                       ..
+  //                 1 ti ti^2 ... ti^q
+  //                       ..
+  //                 1 tn tn^2 ... tn^q]
+  //           q : ordre de regression logistique
+  //           n : nombre d'observations
+  //        3. Y Matrice de la partition floue (les probs a posteriori tik)
+  //           tik = p(zi=k|xi;theta^m); Y de dimensions [nxK] avec K le nombre de classes
+  // Sorties :
+  //
+  //        1. probs : Matrice de dim [nxK] des probabilites p(zi=k;W) de la vaiable zi
+  //          (i=1,...,n)
+  //        2. loglik : logvraisemblance du parametre W du modele logistique
+  //           loglik = Q1(W) = E(l(W;Z)|X;theta^m) = E(p(Z;W)|X;theta^m)
+  //                  = logsum_{i=1}^{n} sum_{k=1}^{K} tik log p(zi=k;W)
+  //
+  // Cette fonction peut egalement ?tre utilis?e pour calculer seulement les
+  // probs de la fa?oc suivante : probs = modele_logit(W,X)
+  //
+  // Faicel Chamroukhi 31 Octobre 2008 (mise ? jour)
+  /////////////////////////////////////////////////////////////////////////////////////////
+
+  unsigned n = X.n_rows;
+  unsigned q = X.n_cols;
+
+  unsigned K = Y.n_cols;
+
+  // Handle different q
+  if (q != W.n_rows) {
+    stop("W must have q + 1 rows and X must have q + 1 columns.");
+  }
+
+  arma::mat Wc = W;
+  // Handle size of K issues
+  if (Wc.n_cols == (K - 1)) { // W doesnt contain the null vector associated with the last class
+    Wc = join_rows(Wc, arma::mat(q, 1, arma::fill::zeros)); // Add the null vector wK for the last component probability
+  }
+  if (Wc.n_cols != K) {
+    stop("W must have K - 1 or K columns.");
+  }
+
+  // Handle different n
+  if ((n != Y.n_cols) && (n != Gamma.n_rows)) {
+    stop("X, Y and Gamma must have the same number of rows which is n.");
+  }
+
+  arma::mat XW(n, K, arma::fill::zeros);
+  arma::colvec maxm(n, arma::fill::zeros);
+  arma::mat expXW(n, K, arma::fill::zeros);
+  arma::mat piik(n, K, arma::fill::zeros);
+  arma::mat GammaMat(n, K, arma::fill::ones);
+
+  GammaMat = Gamma * arma::rowvec(K, arma::fill::ones);
+
+  double loglik;
+
+  XW = X * Wc;
+  maxm = arma::max(XW, 1);
+
+  XW = XW - maxm * arma::rowvec(K, arma::fill::ones); // To avoid overfolow
+
+  double minvalue = -745.1;
+  XW = arma::max(XW, minvalue * arma::mat(XW.n_rows, XW.n_cols, arma::fill::ones));
+  double maxvalue = 709.78;
+  XW = arma::min(XW, maxvalue * arma::mat(XW.n_rows, XW.n_cols, arma::fill::ones));
+  expXW = arma::exp(XW);
+
+  piik = expXW / (arma::sum(expXW, 1) * arma::rowvec(K, arma::fill::ones));
+
+  // log-likelihood
+  loglik = sum(sum((GammaMat % (Y % XW)) - ((GammaMat % Y) % arma::log(arma::sum(expXW, 1) * arma::rowvec(K, arma::fill::ones)))));
+
+  return List::create(Named("loglik") = loglik, Named("piik") = piik);
+
+}
diff --git a/src/multinomialLogit.h b/src/multinomialLogit.h
new file mode 100644
index 0000000..b9e88e3
--- /dev/null
+++ b/src/multinomialLogit.h
@@ -0,0 +1,9 @@
+// Defines a header file containing function signatures for functions in src/
+
+// Protect signatures using an inclusion guard.
+#ifndef multinomialLogit_H
+#define multinomialLogit_H
+
+Rcpp::List multinomialLogit(arma::mat& W, arma::mat& X, arma::mat& Y, arma::mat& Gamma);
+
+#endif
diff --git a/vignettes/.gitignore b/vignettes/.gitignore
new file mode 100644
index 0000000..097b241
--- /dev/null
+++ b/vignettes/.gitignore
@@ -0,0 +1,2 @@
+*.html
+*.R
diff --git a/vignettes/A-quick-tour-of-SNMoE.Rmd b/vignettes/A-quick-tour-of-SNMoE.Rmd
new file mode 100644
index 0000000..51d6936
--- /dev/null
+++ b/vignettes/A-quick-tour-of-SNMoE.Rmd
@@ -0,0 +1,94 @@
+---
+title: "A-quick-tour-of-SNMoE"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{A-quick-tour-of-SNMoE}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r setup, include=FALSE}
+library(knitr)
+knitr::opts_chunk$set(
+	fig.align = "center",
+	fig.height = 5.5,
+	fig.width = 6,
+	warning = FALSE,
+	collapse = TRUE,
+	dev.args = list(pointsize = 10),
+	out.width = "90%",
+	par = TRUE
+)
+knit_hooks$set(par = function(before, options, envir)
+  { if (before && options$fig.show != "none") 
+       par(family = "sans", mar = c(4.1,4.1,1.1,1.1), mgp = c(3,1,0), tcl = -0.5)
+})
+```
+
+```{r, message = FALSE, echo = FALSE}
+library(meteorit)
+```
+
+# Introduction
+
+**MEteorit** 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. 
+It was written in R Markdown, using the [knitr](https://cran.r-project.org/package=knitr) package for production. 
+See `help(package="meteorit")` for further details and references provided by `citation("meteorit")`.
+
+# Load data
+
+```{r}
+data("simulatedstructureddata")
+```
+
+# Set up SNMoE model parameters
+
+```{r}
+K <- 2 # number of regimes (mixture components)
+p <- 1 # dimension of beta (order of the polynomial regressors)
+q <- 1 # dimension of w (order of the logistic regression: to be set to 1 for segmentation)
+```
+
+# Set up EM parameters
+
+```{r}
+n_tries <- 1
+max_iter = 1500
+threshold <- 1e-6
+verbose <- TRUE
+verbose_IRLS <- FALSE
+```
+
+# Estimation
+
+```{r}
+snmoe <- emSNMoE(simulatedstructureddata$X, matrix(simulatedstructureddata$Y), 
+                 K, p, q, n_tries, max_iter, threshold, verbose, verbose_IRLS)
+```
+
+
+# Plots
+
+## Mean curve
+
+```{r}
+snmoe$plot(what = "meancurve")
+```
+
+## Confidence regions
+
+```{r}
+snmoe$plot(what = "confregions")
+```
+
+## Clusters
+
+```{r}
+snmoe$plot(what = "clusters")
+```
+
+## Log-likelihood
+
+```{r}
+snmoe$plot(what = "loglikelihood")
+```
diff --git a/vignettes/A-quick-tour-of-StMoE.Rmd b/vignettes/A-quick-tour-of-StMoE.Rmd
new file mode 100644
index 0000000..a9a9504
--- /dev/null
+++ b/vignettes/A-quick-tour-of-StMoE.Rmd
@@ -0,0 +1,93 @@
+---
+title: "A-quick-tour-of-StMoE"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{A-quick-tour-of-StMoE}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r setup, include=FALSE}
+library(knitr)
+knitr::opts_chunk$set(
+	fig.align = "center",
+	fig.height = 5.5,
+	fig.width = 6,
+	warning = FALSE,
+	collapse = TRUE,
+	dev.args = list(pointsize = 10),
+	out.width = "90%",
+	par = TRUE
+)
+knit_hooks$set(par = function(before, options, envir)
+  { if (before && options$fig.show != "none") 
+       par(family = "sans", mar = c(4.1,4.1,1.1,1.1), mgp = c(3,1,0), tcl = -0.5)
+})
+```
+
+```{r, message = FALSE, echo = FALSE}
+library(meteorit)
+```
+
+# Introduction
+
+**MEteorit** 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. 
+It was written in R Markdown, using the [knitr](https://cran.r-project.org/package=knitr) package for production. 
+See `help(package="meteorit")` for further details and references provided by `citation("meteorit")`.
+
+# Load data
+
+```{r}
+data("simulatedstructureddata")
+```
+
+# Set up StMoE model parameters
+
+```{r}
+K <- 2 # Number of regimes (mixture components)
+p <- 1 # Dimension of beta (order of the polynomial regressors)
+q <- 1 # Dimension of w (order of the logistic regression: to be set to 1 for segmentation)
+
+```
+
+# Set up EM parameters
+
+```{r}
+n_tries <- 1
+max_iter <- 1500
+threshold <- 1e-5
+verbose <- TRUE
+verbose_IRLS <- FALSE
+```
+
+# Estimation
+
+```{r}
+stmoe <- emStMoE(simulatedstructureddata$X, matrix(simulatedstructureddata$Y), K, p, q, n_tries, max_iter, threshold, verbose, verbose_IRLS)
+```
+
+# Plots
+
+## Mean curve
+
+```{r}
+stmoe$plot(what = "meancurve")
+```
+
+## Confidence regions
+
+```{r}
+stmoe$plot(what = "confregions")
+```
+
+## Clusters
+
+```{r}
+stmoe$plot(what = "clusters")
+```
+
+## Log-likelihood
+
+```{r}
+stmoe$plot(what = "loglikelihood")
+```
diff --git a/vignettes/A-quick-tour-of-tMoE.Rmd b/vignettes/A-quick-tour-of-tMoE.Rmd
new file mode 100644
index 0000000..9797a1b
--- /dev/null
+++ b/vignettes/A-quick-tour-of-tMoE.Rmd
@@ -0,0 +1,92 @@
+---
+title: "A-quick-tour-of-tMoE"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{A-quick-tour-of-tMoE}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r setup, include=FALSE}
+library(knitr)
+knitr::opts_chunk$set(
+	fig.align = "center",
+	fig.height = 5.5,
+	fig.width = 6,
+	warning = FALSE,
+	collapse = TRUE,
+	dev.args = list(pointsize = 10),
+	out.width = "90%",
+	par = TRUE
+)
+knit_hooks$set(par = function(before, options, envir)
+  { if (before && options$fig.show != "none") 
+       par(family = "sans", mar = c(4.1,4.1,1.1,1.1), mgp = c(3,1,0), tcl = -0.5)
+})
+```
+
+```{r, message = FALSE, echo = FALSE}
+library(meteorit)
+```
+
+# Introduction
+
+**MEteorit** 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. 
+It was written in R Markdown, using the [knitr](https://cran.r-project.org/package=knitr) package for production. 
+See `help(package="meteorit")` for further details and references provided by `citation("meteorit")`.E")`.
+
+# Load data
+
+```{r}
+data("simulatedstructureddata")
+```
+
+# Set up tMoE model parameters
+
+```{r}
+K <- 2 # Number of regimes (mixture components)
+p <- 1 # Dimension of beta (order of the polynomial regressors)
+q <- 1 # Dimension of w (order of the logistic regression: to be set to 1 for segmentation)
+```
+
+# Set up EM parameters
+
+```{r}
+n_tries <- 1
+max_iter <- 1500
+threshold <- 1e-5
+verbose <- TRUE
+verbose_IRLS <- FALSE
+```
+
+# Estimation
+
+```{r}
+tmoe <- emTMoE(simulatedstructureddata$X, matrix(simulatedstructureddata$Y), K, p, q, n_tries, max_iter, threshold, verbose, verbose_IRLS)
+```
+
+# Plots
+
+## Mean curve
+
+```{r}
+tmoe$plot(what = "meancurve")
+```
+
+## Confidence regions
+
+```{r}
+tmoe$plot(what = "confregions")
+```
+
+## Clusters
+
+```{r}
+tmoe$plot(what = "clusters")
+```
+
+## Log-likelihood
+
+```{r}
+tmoe$plot(what = "loglikelihood")
+```