diff --git a/R/ParamMRHLP.R b/R/ParamMRHLP.R
index 370691e..9ee7cb7 100755
--- a/R/ParamMRHLP.R
+++ b/R/ParamMRHLP.R
@@ -180,7 +180,7 @@ ParamMRHLP <- setRefClass(
         beta[, , k] <<- solve(M, tol = FALSE) %*% t(Xk) %*% yk # Maximization w.r.t betak
         z <- (mData$Y - phi$XBeta %*% beta[, , k]) * (sqrt(weights) %*% ones(1, mData$d))
 
-        # Maximisation w.r.t sigmak (the variances)
+        # Maximisation w.r.t sigma2k (the variances)
         sk <- t(z) %*% z
 
         if (variance_type == "homoskedastic") {
diff --git a/R/emHMMR.R b/R/emHMMR.R
index b1ea1ed..6e74593 100755
--- a/R/emHMMR.R
+++ b/R/emHMMR.R
@@ -95,7 +95,7 @@ emHMMR <- function(X, Y, K, p = 3, variance_type = c("heteroskedastic", "homoske
       stat$loglik <- stat$loglik + log(lambda)
 
       if (verbose) {
-        cat(paste0("EM: Iteration : ", iter, " || log-likelihood : "  , stat$loglik, "\n"))
+        cat(paste0("EM: Iteration : ", iter, " | log-likelihood : "  , stat$loglik, "\n"))
       }
 
       if ((prev_loglik - stat$loglik) > 1e-4) {
@@ -108,7 +108,7 @@ emHMMR <- function(X, Y, K, p = 3, variance_type = c("heteroskedastic", "homoske
       converged <- (abs(stat$loglik - prev_loglik) / abs(prev_loglik) < threshold)
       if (is.na(converged)) {
         converged <- FALSE
-      } # Basically for the first iteration when prev_loglik is Inf
+      } # Basically for the first iteration when prev_loglik is -inf
 
       prev_loglik <- stat$loglik
       stat$stored_loglik <- c(stat$stored_loglik, stat$loglik)
@@ -145,7 +145,7 @@ emHMMR <- function(X, Y, K, p = 3, variance_type = c("heteroskedastic", "homoske
   # Smoothing state sequences : argmax(smoothing probs), and corresponding binary allocations partition
   statSolution$MAP()
 
-  # Finish the computation of statistics
+  # Finish the computation of the model statistics
   statSolution$computeStats(paramSolution)
 
   return(ModelHMMR(param = paramSolution, stat = statSolution))
diff --git a/R/emMHMMR.R b/R/emMHMMR.R
index 6e9bb3b..a565f46 100755
--- a/R/emMHMMR.R
+++ b/R/emMHMMR.R
@@ -111,7 +111,7 @@ emMHMMR <- function(X, Y, K, p = 3, variance_type = c("heteroskedastic", "homosk
       converged <- (abs(stat$loglik - prev_loglik) / abs(prev_loglik) < threshold)
       if (is.na(converged)) {
         converged <- FALSE
-      } # Basically for the first iteration when prev_loglik is Inf
+      } # Basically for the first iteration when prev_loglik is -inf
 
       prev_loglik <- stat$loglik
       stat$stored_loglik <- c(stat$stored_loglik, stat$loglik)
@@ -147,7 +147,7 @@ emMHMMR <- function(X, Y, K, p = 3, variance_type = c("heteroskedastic", "homosk
   # Smoothing state sequences : argmax(smoothing probs), and corresponding binary allocations partition
   statSolution$MAP()
 
-  # Finish the computation of statistics
+  # Finish the computation of model statistics
   statSolution$computeStats(paramSolution)
 
   return(ModelMHMMR(param = paramSolution, stat = statSolution))
diff --git a/R/emMRHLP.R b/R/emMRHLP.R
index 02117aa..2901971 100644
--- a/R/emMRHLP.R
+++ b/R/emMRHLP.R
@@ -103,7 +103,7 @@ emMRHLP <- function(X, Y, K, p = 3, q = 1, variance_type = c("heteroskedastic",
       converge <- abs((stat$loglik - prev_loglik) / prev_loglik) <= threshold
       if (is.na(converge)) {
         converge <- FALSE
-      } # Basically for the first iteration when prev_loglik is Inf
+      } # Basically for the first iteration when prev_loglik is -inf
 
       prev_loglik <- stat$loglik
       stat$stored_loglik <- c(stat$stored_loglik, stat$loglik)
diff --git a/R/emRHLP.R b/R/emRHLP.R
index 50e5077..b85b243 100755
--- a/R/emRHLP.R
+++ b/R/emRHLP.R
@@ -98,7 +98,7 @@ emRHLP <- function(X, Y, K, p = 3, q = 1, variance_type = c("heteroskedastic", "
       converge <- abs((stat$loglik - prev_loglik) / prev_loglik) <= threshold
       if (is.na(converge)) {
         converge <- FALSE
-      } # Basically for the first iteration when prev_loglik is Inf
+      } # Basically for the first iteration when prev_loglik is -inf
 
       prev_loglik <- stat$loglik
       stat$stored_loglik <- c(stat$stored_loglik, stat$loglik)
diff --git a/README.Rmd b/README.Rmd
index fedcd93..ac31ef7 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -299,18 +299,17 @@ mhmmr$plot(what = c("smoothed", "regressors", "loglikelihood"))
 
 # Model selection
 
-In this package, it is possible to select models based on information criteria 
-such as **BIC**, **AIC** and **ICL**.
+samurais also implements model selection procedures to select an optimal model based on information criteria including **BIC**, **AIC** and **ICL**.
 
 The selection can be done for the two following parameters:
 
- * **K**: The number of regimes;
+ * **K**: The number of regimes (segments);
  * **p**: The order of the polynomial regression.
-
+Instructions below can be used to illustrate the model on provided simulated and real data sets.
 <details>
   <summary>RHLP</summary>
 
-Let's select a RHLP model for the following time series **Y**:
+Let's select a RHLP model for the following time series:
 
 ```{r, message = FALSE}
 data("univtoydataset")
@@ -331,7 +330,7 @@ selectedrhlp$plot(what = "estimatedsignal")
 <details>
   <summary>HMMR</summary>
 
-Let's select a HMMR model for the following time series **Y**:
+Let's select a HMMR model for the following time series:
 
 ```{r, message = FALSE}
 data("univtoydataset")
@@ -352,7 +351,7 @@ selectedhmmr$plot(what = "smoothed")
 <details>
 <summary>MRHLP</summary>
 
-Let's select a MRHLP model for the following multivariate time series **Y**:
+Let's select a MRHLP model for the following multivariate time series:
 
 <br />
 
@@ -374,7 +373,7 @@ selectedmrhlp$plot(what = "estimatedsignal")
 <details>
   <summary>MHMMR</summary>
 
-Let's select a MHMMR model for the following multivariate time series **Y**:
+Let's select a MHMMR model for the following multivariate time series:
 
 ```{r}
 data("multivtoydataset")