update NMoE model with vigniette and readme

fchamroukhi · Jul 16, 2019 · 587e4da · 587e4da
1 parent 19765f1
commit 587e4da
Show file tree

Hide file tree

Showing 29 changed files with 416 additions and 48 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -1,5 +1,22 @@
 # Generated by roxygen2: do not edit by hand
 
+export(emNMoE)
+export(emSNMoE)
+export(emStMoE)
+export(emTMoE)
+exportClasses(FData)
+exportClasses(ModelNMoE)
+exportClasses(ModelSNMoE)
+exportClasses(ModelStMoE)
+exportClasses(ModelTMoE)
+exportClasses(ParamNMoE)
+exportClasses(ParamSNMoE)
+exportClasses(ParamStMoE)
+exportClasses(ParamTMoE)
+exportClasses(StatNMoE)
+exportClasses(StatSNMoE)
+exportClasses(StatStMoE)
+exportClasses(StatTMoE)
 import(methods)
 importFrom(Rcpp,sourceCpp)
 importFrom(pracma,fzero)

diff --git a/R/ParamNMoE.R b/R/ParamNMoE.R
@@ -15,10 +15,9 @@
 #' @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
 ParamNMoE <- setRefClass(
-  "ParamTMoE",
+  "ParamNMoE",
   fields = list(
     fData = "FData",
     phiBeta = "list",
@@ -106,23 +105,23 @@ ParamNMoE <- setRefClass(
 
     MStep = function(statNMoE, verbose_IRLS) {
       # M-Step
-      res_irls <- IRLS(phiAlpha$XBeta, statTMoE$tik, ones(nrow(statTMoE$tik), 1), alpha, verbose_IRLS)
-      statTMoE$piik <- res_irls$piik
+      res_irls <- IRLS(phiAlpha$XBeta, statNMoE$tik, ones(nrow(statNMoE$tik), 1), alpha, verbose_IRLS)
+      statNMoE$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] %*% ones(1, p + 1)))
-        yk <- fData$Y * sqrt(statTMoE$tik[,k])
+        Xbeta <- phiBeta$XBeta * sqrt(statNMoE$tik[,k] %*% ones(1, p + 1))
+        yk <- fData$Y * sqrt(statNMoE$tik[,k])
 
         #update the regression coefficients
         beta[, k] <<- solve((t(Xbeta) %*% Xbeta)) %*% (t(Xbeta) %*% yk)
 
         # update the variances sigma2k
-        sigma[k] <<- sum(statTMoE$tik[, k] * ((fData$Y - phiBeta$XBeta %*% beta[, k])^2)) / sum(statTMoE$tik[,k])
+        sigma[k] <<- sum(statNMoE$tik[, k] * ((fData$Y - phiBeta$XBeta %*% beta[, k])^2)) / sum(statNMoE$tik[,k])
 
       }
 

diff --git a/R/StatNMoE.R b/R/StatNMoE.R
@@ -1,6 +1,6 @@
-#' A Reference Class which contains statistics of a TMoE model.
+#' A Reference Class which contains statistics of a NMoE model.
 #'
-#' StatTMoE contains all the parameters of a [TMoE][ParamTMoE] model.
+#' StatNoE contains all the parameters of a [NMoE][ParamNMoE] 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
@@ -40,10 +40,10 @@
 #' @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]
+#' @seealso [ParamNMoE], [FData]
 #' @export
-StatTMoE <- setRefClass(
-  "StatTMoE",
+StatNMoE <- setRefClass(
+  "StatNMoE",
   fields = list(
     piik = "matrix",
     z_ik = "matrix",
@@ -65,24 +65,24 @@ StatTMoE <- setRefClass(
     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)
+    initialize = function(paramNMoE = ParamNMoE()) {
+      piik <<- matrix(NA, paramNMoE$fData$n, paramNMoE$K)
+      z_ik <<- matrix(NA, paramNMoE$fData$n, paramNMoE$K)
+      klas <<- matrix(NA, paramNMoE$fData$n, 1)
+      Wik <<- matrix(0, paramNMoE$fData$n * paramNMoE$fData$m, paramNMoE$K)
+      Ey_k <<- matrix(NA, paramNMoE$fData$n, paramNMoE$K)
+      Ey <<- matrix(NA, paramNMoE$fData$n, 1)
+      Var_yk <<- matrix(NA, 1, paramNMoE$K)
+      Vary <<- matrix(NA, paramNMoE$fData$n, 1)
       log_lik <<- -Inf
       com_loglik <<- -Inf
       stored_loglik <<- numeric()
       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)
+      log_piik_fik <<- matrix(0, paramNMoE$fData$n, paramNMoE$K)
+      log_sum_piik_fik <<- matrix(NA, paramNMoE$fData$n, 1)
+      tik <<- matrix(0, paramNMoE$fData$n, paramNMoE$K)
     },
 
     MAP = function() {
@@ -124,55 +124,55 @@ StatTMoE <- setRefClass(
     #######
     # compute the final solution stats
     #######
-    computeStats = function(paramTMoE) {
+    computeStats = function(paramNMoE) {
 
       # E[yi|zi=k]
-      Ey_k <<- paramTMoE$phiBeta$XBeta[1:paramTMoE$fData$n, ] %*% paramTMoE$beta
+      Ey_k <<- paramNMoE$phiBeta$XBeta[1:paramNMoE$fData$n, ] %*% paramNMoE$beta
 
       # E[yi]
       Ey <<- matrix(apply(piik * Ey_k, 1, sum))
 
       # Var[yi|zi=k]
-      Var_yk <<- paramTMoE$sigma
+      Var_yk <<- paramNMoE$sigma
 
       # Var[yi]
-      Vary <<- apply(piik * (Ey_k ^ 2 + ones(paramTMoE$fData$n, 1) %*% Var_yk), 1, sum) - Ey ^ 2
+      Vary <<- apply(piik * (Ey_k ^ 2 + ones(paramNMoE$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
+      BIC <<- log_lik - (paramNMoE$nu * log(paramNMoE$fData$n * paramNMoE$fData$m) / 2)
+      AIC <<- log_lik - paramNMoE$nu
       ## CL(theta) : complete-data loglikelihood
-      zik_log_piik_fk <- (repmat(z_ik, paramTMoE$fData$m, 1)) * log_piik_fik
+      zik_log_piik_fk <- (repmat(z_ik, paramNMoE$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)
+      ICL <<- com_loglik - (paramNMoE$nu * log(paramNMoE$fData$n * paramNMoE$fData$m) / 2)
       # solution.XBeta = XBeta(1:m,:);
       # solution.XAlpha = XAlpha(1:m,:);
     },
     #######
     # EStep
     #######
     EStep = function(paramNMoE) {
-      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)
+      piik <<- multinomialLogit(paramNMoE$alpha, paramNMoE$phiAlpha$XBeta, ones(paramNMoE$fData$n, paramNMoE$K), ones(paramNMoE$fData$n, 1))$piik
+      piik_fik <- zeros(paramNMoE$fData$m * paramNMoE$fData$n, paramNMoE$K)
 
-      for (k in (1:paramTMoE$K)) {
-        muk <- paramTMoE$phiBeta$XBeta %*% paramTMoE$beta[, k]
+      for (k in (1:paramNMoE$K)) {
+        muk <- paramNMoE$phiBeta$XBeta %*% paramNMoE$beta[, k]
 
-        sigma2k <- paramTMoE$sigma[k]
+        sigma2k <- paramNMoE$sigma[k]
 
-        log_piik_fik[, k] <<- log(piik[, k]) -0.5 * log(2 * pi) - 0.5 * log(sigma2k) - 0.5 * ((paramTMoE$fData$Y - muk)^2)/sigma2k
+        log_piik_fik[, k] <<- log(piik[, k]) -0.5 * log(2 * pi) - 0.5 * log(sigma2k) - 0.5 * ((paramNMoE$fData$Y - muk)^2)/sigma2k
       }
 
       log_sum_piik_fik <<- matrix(log(rowSums(exp(log_piik_fik))))
 
-      log_tik <-  log_piik_fik - log_sum_piik_fik %*% ones(1,paramTMoE$K)
+      log_tik <-  log_piik_fik - log_sum_piik_fik %*% ones(1,paramNMoE$K)
       ttik <- exp(log_tik)
 
-      tik <<-  ttik / (rowSums(ttik) %*% ones(1,paramTMoE$K))
+      tik <<-  ttik / (rowSums(ttik) %*% ones(1,paramNMoE$K))
     }
   )
 )
diff --git a/README.Rmd b/README.Rmd
@@ -134,5 +134,31 @@ tmoe$plot()
 
 </details>
 
+
+
+<details>
+<summary>NMoE</summary>
+```{r, message = FALSE}
+library(meteorits)
+
+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-5
+verbose <- TRUE
+verbose_IRLS <- FALSE
+
+nmoe <- emNMoE(simulatedstructureddata$X, matrix(simulatedstructureddata$Y), K, p, q, n_tries, max_iter, threshold, verbose, verbose_IRLS)
+
+nmoe$plot()
+```
+
+</details>
+
 # References
 
diff --git a/README.md b/README.md
@@ -124,6 +124,32 @@ tmoe$plot()
 
 </details>
 
+<details> <summary>NMoE</summary>
+
+``` r
+library(meteorits)
+
+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-5
+verbose <- TRUE
+verbose_IRLS <- FALSE
+
+nmoe <- emNMoE(simulatedstructureddata$X, matrix(simulatedstructureddata$Y), K, p, q, n_tries, max_iter, threshold, verbose, verbose_IRLS)
+
+nmoe$plot()
+```
+
+<img src="man/figures/README-unnamed-chunk-8-1.png" style="display: block; margin: auto;" /><img src="man/figures/README-unnamed-chunk-8-2.png" style="display: block; margin: auto;" /><img src="man/figures/README-unnamed-chunk-8-3.png" style="display: block; margin: auto;" /><img src="man/figures/README-unnamed-chunk-8-4.png" style="display: block; margin: auto;" />
+
+</details>
+
 References
 ==========
 

diff --git a/man/ModelNMoE-class.Rd b/man/ModelNMoE-class.Rd
diff --git a/man/ModelTMoE-class.Rd b/man/ModelTMoE-class.Rd
diff --git a/man/ParamNMoE-class.Rd b/man/ParamNMoE-class.Rd