Nothing
forward <- function(X, reduced, tol, max.iter){
if(!is.matrix(X)){
if(is.vector(X)){
n <- length(X)
X <- matrix(X, n, 1)
}
}
K <- dim(reduced$Mu)[1]
p <- dim(reduced$Mu)[2]
n <- dim(X)[1]
la.reduced <- reduced$la
ind <- as.vector(la.reduced)==0
M <- K - 1 + 2 * K * p + K * p * (p + 1) / 2 - sum(ind)
reduced.bic <- Manly.bic(reduced$ll, n, M)
forward.list <- list()
forward.bic <- rep(NA, sum(ind))
tick <- 0
for(i in (1:length(ind))[ind]){
tick <- tick + 1
la <- la.reduced
for(step in seq(0.01, 0.1, 0.01)){
la[i] <- step
res1 <- try(forward <- Manly.EM(X, id = reduced$id, la = la, tol = tol, max.iter = max.iter))
if(!inherits(res1, "try-error")){
if(!is.null(forward$ll)){
break
}
}
else {
la[i] <- -step
res1 <- try(forward <- Manly.EM(X, id = reduced$id, la = la, tol = tol, max.iter = max.iter))
if(!inherits(res1, "try-error")){
if(!is.null(forward$ll)){
break
}
}
}
}
if(inherits(res1, "try-error") || (is.null(forward$ll))){
forward.list[[tick]] <- NULL
forward.bic[tick] <- Inf
}
else{
forward.list[[tick]] <- forward
forward.bic[tick] <- Manly.bic(forward$ll, n, M + 1)
}
}
return(list(reduced.bic = reduced.bic, forward.bic = forward.bic, forward.list = forward.list))
}
backward <- function(X, full, tol, max.iter){
if(!is.matrix(X)){
if(is.vector(X)){
n <- length(X)
X <- matrix(X, n, 1)
}
}
K <- dim(full$Mu)[1]
p <- dim(full$Mu)[2]
n <- dim(X)[1]
la.full <- full$la
ind <- as.vector(la.full)!=0
M <- K - 1 + K * p + K * p * (p + 1) / 2 + sum(ind)
full.bic <- Manly.bic(full$ll, n, M)
reduced.list <- list()
reduced.bic <- rep(NA, sum(ind))
tick <- 0
for(i in (1:length(ind))[ind]){
tick <- tick + 1
la <- la.full
la[i] <- 0.0
res1 <- try(reduced <- Manly.EM(X, id = full$id, la = la, tol = tol, max.iter = max.iter))
if(inherits(res1, "try-error") || (is.null(reduced$ll))){
reduced.list[[tick]] <- NULL
reduced.bic[tick] <- Inf
}
else{
reduced.list[[tick]] <- reduced
reduced.bic[tick] <- Manly.bic(reduced$ll, n, M - 1)
}
}
return(list(full.bic = full.bic, reduced.bic = reduced.bic, reduced.list = reduced.list))
}
Manly.select <- function(X, model, method, tol = 1e-5, max.iter = 1000, silent = FALSE){
step <- 0
if (tol <= 0) stop("Wrong value of tol...\n")
if (max.iter < 1) stop("Wrong number of iterations iter...\n")
if(method == "backward"){
if(!silent){
repeat{
step <- step + 1
C <- backward(X, model, tol, max.iter)
comparison <- C$reduced.bic < C$full.bic
if (length(C$reduced.bic != 0)){
cat("step", step, ":\n\tcurrent BIC =", C$full.bic, "\n\talternative BICs =", C$reduced.bic, "\n")
} else {
cat("step", step, ":\n\tcurrent BIC =", C$full.bic, "\n\talternative BICs =", NA, "\n")
}
if(sum(comparison) != 0){
index <- which.min(C$reduced.bic)
model <- C$reduced.list[[index]]
BIC <- C$reduced.bic[index]
} else {
BIC <- C$full.bic
break
}
}
}
else{
repeat{
step <- step + 1
C <- backward(X, model, tol, max.iter)
comparison <- C$reduced.bic < C$full.bic
if(sum(comparison)!=0){
index <- which.min(C$reduced.bic)
model <- C$reduced.list[[index]]
BIC <- C$reduced.bic[index]
} else {
BIC <- C$full.bic
break
}
}
}
}
else if(method == "forward"){
if(!silent){
repeat{
step <- step + 1
C <- forward(X, model, tol, max.iter)
comparison <- C$forward.bic < C$reduced.bic
if (length(C$forward.bic != 0)){
cat("step", step, ":\n\tcurrent BIC =", C$reduced.bic, "\n\talternative BICs =", C$forward.bic, "\n")
} else {
cat("step", step, ":\n\tcurrent BIC =", C$reduced.bic, "\n\talternative BICs =", NA, "\n")
}
if(sum(comparison) != 0){
index <- which.min(C$forward.bic)
model <- C$forward.list[[index]]
BIC <- C$forward.bic[index]
} else {
BIC <- C$reduced.bic
break
}
}
}
else{
repeat{
step <- step + 1
C <- forward(X, model, tol, max.iter)
comparison <- C$forward.bic < C$reduced.bic
if(sum(comparison) != 0){
index <- which.min(C$forward.bic)
model <- C$forward.list[[index]]
BIC <- C$forward.bic[index]
} else {
BIC <- C$reduced.bic
break
}
}
}
}
return(list(la = model$la, tau = model$tau, Mu = model$Mu, S = model$S, gamma = model$gamma, id = model$id, ll = model$ll, bic = BIC, iter = model$iter, flag = model$flag))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.