Skip to content

Commit

Permalink
Change name of field 'nuk' in ParamStMoE class
Browse files Browse the repository at this point in the history 
- Change name of fields 'nuk' ('nuk' is changed into 'nu') in ParamStMoE class;
- Update vignettes (Add information about the models).
  • Loading branch information
@Florian-Lecocq
Florian-Lecocq committed Sep 16, 2019
1 parent bb60a08 commit 07c449c
Show file tree
Hide file tree
Showing 42 changed files with 749 additions and 881 deletions.
10 changes: 5 additions & 5 deletions R/ParamStMoE.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@
#' @field delta delta is equal \eqn{\delta =
#' \frac{\lambda}{\sqrt{1+\lambda^2}}}{\delta = \lambda /
#' (1+\lambda^2)^(1/2)}.
#' @field nuk The degree of freedom for the Student distribution for each
#' @field nu The degree of freedom for the Student distribution for each
#' experts (matrix of size \eqn{(1, K)}).
#' @field df The degree of freedom of the StMoE model representing the
#' complexity of the model.
Expand All @@ -53,7 +53,7 @@ ParamStMoE <- setRefClass(
sigma2 = "matrix",
lambda = "matrix",
delta = "matrix",
nuk = "matrix"
nu = "matrix"
),
methods = list(
initialize = function(X = numeric(), Y = numeric(1), K = 1, p = 3, q = 1) {
Expand All @@ -74,7 +74,7 @@ ParamStMoE <- setRefClass(
sigma2 <<- matrix(NA, 1, K)
lambda <<- matrix(NA, ncol = K)
delta <<- matrix(NA, ncol = K)
nuk <<- matrix(NA, ncol = K)
nu <<- matrix(NA, ncol = K)
},

initParam = function(segmental = FALSE) {
Expand Down Expand Up @@ -135,7 +135,7 @@ ParamStMoE <- setRefClass(
delta <<- lambda / sqrt(1 + lambda ^ 2)

# Intitialization of the degrees of freedom
nuk <<- 50 * rand(1, K)
nu <<- 50 * rand(1, K)
},

MStep = function(statStMoE, calcAlpha = FALSE, calcBeta = FALSE, calcSigma2 = FALSE, calcLambda = FALSE, calcNu = FALSE, verbose_IRLS = FALSE) {
Expand Down Expand Up @@ -195,7 +195,7 @@ ParamStMoE <- setRefClass(
if (calcNu) {

for (k in 1:K) {
try(nuk[k] <<- suppressWarnings(uniroot(f <- function(nnu) {
try(nu[k] <<- suppressWarnings(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, 100))$root), silent = TRUE)
Expand Down
22 changes: 11 additions & 11 deletions R/StatStMoE.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,7 +134,7 @@ StatStMoE <- setRefClass(
parameters provided by the object \\code{paramStMoE} of class
\\link{ParamStMoE}."

Xi_nuk = sqrt(paramStMoE$nuk / pi) * (gamma(paramStMoE$nuk / 2 - 1 / 2)) / (gamma(paramStMoE$nuk / 2))
Xi_nuk = sqrt(paramStMoE$nu / pi) * (gamma(paramStMoE$nu / 2 - 1 / 2)) / (gamma(paramStMoE$nu / 2))

# E[yi|xi,zi=k]
Ey_k <<- paramStMoE$phiBeta$XBeta[1:paramStMoE$n,] %*% paramStMoE$beta + ones(paramStMoE$n, 1) %*% (paramStMoE$delta * sqrt(paramStMoE$sigma2) * Xi_nuk)
Expand All @@ -143,7 +143,7 @@ StatStMoE <- setRefClass(
Ey <<- matrix(apply(piik * Ey_k, 1, sum))

# Var[yi|xi,zi=k]
Var_yk <<- (paramStMoE$nuk / (paramStMoE$nuk - 2) - (paramStMoE$delta ^ 2) * (Xi_nuk ^ 2)) * paramStMoE$sigma2
Var_yk <<- (paramStMoE$nu / (paramStMoE$nu - 2) - (paramStMoE$delta ^ 2) * (Xi_nuk ^ 2)) * paramStMoE$sigma2

# Var[yi|xi]
Vary <<- apply(piik * (Ey_k ^ 2 + ones(paramStMoE$n, 1) %*% Var_yk), 1, sum) - Ey ^ 2
Expand All @@ -170,8 +170,8 @@ StatStMoE <- setRefClass(

for (k in (1:paramStMoE$K)) {
dik[, k] <<- (paramStMoE$Y - paramStMoE$phiBeta$XBeta %*% paramStMoE$beta[, k]) / sqrt(paramStMoE$sigma2[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 / sqrt(paramStMoE$sigma2[k])) * dt(dik[, k], paramStMoE$nuk[k]) * pt(mik[, k], paramStMoE$nuk[k] + 1)
mik[, k] <- paramStMoE$lambda[k] %*% dik[, k] * sqrt(paramStMoE$nu[k] + 1) / (paramStMoE$nu[k] + dik[, k] ^ 2)
piik_fik[, k] <- piik[, k] * (2 / sqrt(paramStMoE$sigma2[k])) * dt(dik[, k], paramStMoE$nu[k]) * pt(mik[, k], paramStMoE$nu[k] + 1)
}

stme_pdf <<- matrix(rowSums(piik_fik)) # Skew-t mixture of experts density
Expand All @@ -196,7 +196,7 @@ StatStMoE <- setRefClass(
for (k in (1:paramStMoE$K)) {

fx = function(x) {
return((psigamma((paramStMoE$nuk[k] + 2) / 2) - psigamma((paramStMoE$nuk[k] + 1) / 2) + log(1 + (x ^ 2) / (paramStMoE$nuk[k])) + ((paramStMoE$nuk[k] + 1) * x ^ 2 - paramStMoE$nuk[k] - 1) / ((paramStMoE$nuk[k] + 1) * (paramStMoE$nuk[k] + 1 + x ^ 2))) * dt(x, paramStMoE$nuk[k] + 1))
return((psigamma((paramStMoE$nu[k] + 2) / 2) - psigamma((paramStMoE$nu[k] + 1) / 2) + log(1 + (x ^ 2) / (paramStMoE$nu[k])) + ((paramStMoE$nu[k] + 1) * x ^ 2 - paramStMoE$nu[k] - 1) / ((paramStMoE$nu[k] + 1) * (paramStMoE$nu[k] + 1 + x ^ 2))) * dt(x, paramStMoE$nu[k] + 1))
}


Expand All @@ -206,36 +206,36 @@ StatStMoE <- setRefClass(


dik[, k] <<- (paramStMoE$Y - muk) / sigmak
mik[, k] <- paramStMoE$lambda[k] %*% dik[, k] * sqrt((paramStMoE$nuk[k] + 1) / (paramStMoE$nuk[k] + dik[, k] ^ 2))
mik[, k] <- paramStMoE$lambda[k] %*% dik[, k] * sqrt((paramStMoE$nu[k] + 1) / (paramStMoE$nu[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)
wik[, k] <<- ((paramStMoE$nu[k] + 1) / (paramStMoE$nu[k] + dik[, k] ^ 2)) * pt(mik[, k] * sqrt((paramStMoE$nu[k] + 3) / (paramStMoE$nu[k] + 1)), paramStMoE$nu[k] + 3) / pt(mik[, k], paramStMoE$nu[k] + 1)


if (calcE1) {

univStMoEpdf(paramStMoE)
# E[Wi Ui |yi,zik=1]
E1ik[, k] <<- deltak * abs(paramStMoE$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)))
(sqrt(1 - deltak ^ 2) / (pi * stme_pdf)) * ((dik[, k] ^ 2 / (paramStMoE$nu[k] * (1 - deltak ^ 2)) + 1) ^ (-(paramStMoE$nu[k] / 2 + 1)))
}

if (calcE2) {

E2ik[, k] <<- deltak ^ 2 * ((paramStMoE$Y - muk) ^ 2) * wik[, k] +
(1 - deltak ^ 2) * sigmak ^ 2 + ((deltak * (paramStMoE$Y - muk) * sqrt(1 - deltak ^ 2)) / (pi * stme_pdf)) * (((dik[, k] ^ 2) / (paramStMoE$nuk[k] * (1 - deltak ^ 2)) + 1) ^ (-(paramStMoE$nuk[k] / 2 + 1)))
(1 - deltak ^ 2) * sigmak ^ 2 + ((deltak * (paramStMoE$Y - muk) * sqrt(1 - deltak ^ 2)) / (pi * stme_pdf)) * (((dik[, k] ^ 2) / (paramStMoE$nu[k] * (1 - deltak ^ 2)) + 1) ^ (-(paramStMoE$nu[k] / 2 + 1)))
}

if (calcE3) {

Integgtx[, k] <- sapply(mik[, k], function(x) try(integrate(f = fx, lower = -Inf, upper = x)$value, silent = TRUE))
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[, k]
E3ik[, k] <<- wik[, k] - log((paramStMoE$nu[k] + dik[, k] ^ 2) / 2) - (paramStMoE$nu[k] + 1) / (paramStMoE$nu[k] + dik[, k] ^ 2) + psigamma((paramStMoE$nu[k] + 1) / 2) + ((paramStMoE$lambda[k] * dik[, k] * (dik[, k] ^ 2 - 1)) / sqrt((paramStMoE$nu[k] + 1) * ((paramStMoE$nu[k] + dik[, k] ^ 2) ^ 3))) * dt(mik[, k], paramStMoE$nu[k] + 1) / pt(mik[, k], paramStMoE$nu[k] + 1) + (1 / pt(mik[, k], paramStMoE$nu[k] + 1)) * Integgtx[, k]
}

if (calcTau) {

# 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)
piik_fik[, k] <- piik[, k] * (2 / sigmak) * dt(dik[, k], paramStMoE$nu[k]) * pt(mik[, k], paramStMoE$nu[k] + 1)

}

Expand Down
12 changes: 6 additions & 6 deletions README.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,10 +31,10 @@ contains sparse mixture-of-experts models for high-dimensional data.

Our (dis-)covered meteorits are for instance the following ones:

* NMoE;
* tMoE;
* SNMoE;
* StMoE.
* NMoE (Normal Mixtures-of-Experts);
* tMoE (t Mixtures-of-Experts);
* SNMoE (Skew-Normal Mixtures-of-Experts);
* StMoE (Skew t Mixtures-of-Experts).

The models and algorithms are developped and written in Matlab by Faicel
Chamroukhi, and translated and designed into R packages by Florian Lecocq,
Expand Down Expand Up @@ -74,7 +74,7 @@ library(meteorits)
<summary>NMoE</summary>

```{r, echo = TRUE}
# Application to a simuated data set
# Application to a simulated data set
n <- 500 # Size of the sample
alphak <- matrix(c(0, 8), ncol = 1) # Parameters of the gating network
Expand Down Expand Up @@ -131,7 +131,7 @@ nmoe$plot()
<summary>TMoE</summary>

```{r, echo = TRUE}
# Application to a simuated data set
# Application to a simulated data set
n <- 500 # Size of the sample
alphak <- matrix(c(0, 8), ncol = 1) # Parameters of the gating network
Expand Down
Loading

0 comments on commit 07c449c

Please sign in to comment.