R/Winitial.R

In CorReg: Linear Regression Based on Linear Structure Between Variables

# ' initialization based on a weight matrix (correlation or other)
# '
# ' @param W Weight matrix with values between 0 and 1 ( e.g.: abs(cor(X)) )
# ' @param X the dataset (if BIC=TRUE)
# ' @param Bic_null_vect the BIC of the null hypothesis (used for independent variables)
# ' @param relax TRUE: constraint relaxation, FALSE: reject mode
# ' @param p1max maximum complexity for a regression
# ' @param random if FALSE W is not used.
# ' @param mode candidates strategy
# ' @param p2max maximum number of regressions
# ' @param sorted boolean to sort the candidates based on W (best first)
# ' @param BIC boolean to use BIC for probabilities.
# ' @param nbclustmax parameter for calcul_BIC_mixmod if needed
# ' @export
Winitial <- function(W = W, X = NULL, p1max = NULL, BIC = FALSE, Bic_null_vect = NULL, relax = TRUE, random = TRUE,
                     nbclustmax = 5, sorted = TRUE, mode = c("all", "sorted", "multinom"), p2max = NULL) {
  #W est la matrice des poids (entre 0 et 1)
  #X est la matrice des donnees (ne sert que si BIC=T)
  #Bic_null_vect est le vecteur bic vide (issu de mixmod par exemple)
  #relax indique si on fait la relaxation ou le rejet
  #p1max est le nombre max de 1 sur une colonne de Z
  #random indique si on fait le tirage au sort ou pas (si on utilise BIc on a moins besoin du tirage). sans tirage = sans tenir compte de W
  p = ncol(W)
  if (is.null(p2max)) {
    p2max = p
  }
  diag(W) = 0 # on tue la diagonale au cas ou ce serait une matrice de corr?lation
  W = abs(W)
  
  mode = match.arg(mode)
  
  Z = matrix(0, ncol = p, nrow = p)
  list_j = sample(p) # melange les entiers de 1 a p
  
  if (BIC == TRUE &  !is.null(X)) {
    if (is.null(Bic_null_vect)) {
      Bic_null_vect = density_estimation(X = X)$BIC_vect
    }
    BIC_opt_vect = Bic_null_vect
    BIC_opt = sum(Bic_null_vect) # initialisation de l'optimum
  } else{
    BIC = FALSE # soit il l'etait deja soit il manque de quoi le calculer
  }
  
  lambda = 0
  if (is.null(p1max)) {
    p1max = p + 1
  } # pour etre tranquille, revient a une borne infinie
  
  if (mode == "all") {
    # on teste tous les points
    for (j in 1:p) {
      jloc = list_j[j]
      list_i = sample(p)#melange les entiers de 1 a p
      for (i in 1:p) {
        if (list_i[i] != jloc) {
          # si on n'est pas sur la diagonale
          iloc = list_i[i]
          
          if (random) {
            lambda = runif(1)
          } # on tire un nombre aleatoire entre 0 et 1
          if (!random |
              lambda < W[iloc, jloc]) {
            # si le poids est interessant ou bien pas d'alea
            if (sum(Z[, jloc]) < p1max) {
              # si on respecte la contrainte de complexite
              croisement = FALSE
              if (sum(Z[jloc, ]) != 0 | sum(Z[, iloc]) != 0) {
                croisement = TRUE
              }
              if (croisement == FALSE |
                  relax == TRUE) {
                # si pas croisement ou relax
                Zloc = Z
                Zloc[iloc, jloc] = 1
                if (croisement == TRUE) {
                  # donc croisement et relax
                  # faire la relax
                  Zloc[jloc, ] = 0
                  Zloc[, iloc] = 0
                }
                if (length(which(colSums(Zloc) != 0)) <= p2max) {
                  if (BIC) {
                    # si on veut s'appuyer sur le bic
                    
                    BICloc_vect = BicZ(X = X, Z = Zloc, Bic_null_vect = Bic_null_vect, Bic_old = BIC_opt_vect, Zold = Z)
                    BICloc = sum(BICloc_vect)
                    
                    if (BICloc < BIC_opt) {
                      # modification effective
                      Z = Zloc
                      BIC_opt_vect = BICloc_vect
                      BIC_opt = BICloc
                    }
                  } else{
                    # on se fout du bic
                    Z = Zloc
                  }
                }
              }
            }
          }
        }
      }
    }
  } else if (mode == "sorted") {
    # on teste la colonne puis on y place un 1 selon les probas de la colonne (on priorise les meilleurs donc moins de bruit)
    for (j in 1:p) {
      jloc = list_j[j]
      # on teste la colonne en fonction de sa plus grande proba
      if (random) {
        go = rbinom(1, 1, max(W[, jloc]))
      } else{
        go = 1
      }
      if (go == 1) {
        # si la colonne est jugee interessante
        i = 1
        list_i = 1:p
        list_i = list_i[-jloc] # on enleve la diagonale
        list_i = list_i[order(W[-jloc, jloc], decreasing = TRUE)] # on ordonne selon les poids
        while (go & length(list_i) > 0) {
          iloc = list_i[1]#on choisit le meilleur
          if (sum(Z[, jloc]) < p1max) {
            # si on respecte la contrainte de complexite
            croisement = FALSE
            if (sum(Z[jloc, ]) != 0 | sum(Z[, iloc]) != 0) {
              croisement = TRUE
            }
            if (croisement == FALSE |
                relax == TRUE) {
              # si pas croisement ou relax
              Zloc = Z
              Zloc[iloc, jloc] = 1
              if (croisement == TRUE) {
                #donc croisement et relax
                #faire la relax
                Zloc[jloc, ] = 0
                Zloc[, iloc] = 0
              }
              if (length(which(colSums(Zloc) != 0)) <= p2max) {
                if (BIC) {
                  # si on veut s'appuyer sur le bic
                  
                  BICloc_vect = BicZ(X = X, Z = Zloc, Bic_null_vect = Bic_null_vect, Bic_old = BIC_opt_vect, Zold = Z)
                  BICloc = sum(BICloc_vect)
                  if (BICloc < BIC_opt) {
                    # modification effective
                    Z = Zloc
                    BIC_opt_vect = BICloc_vect
                    BIC_opt = BICloc
                  }
                } else{
                  # on se fout du bic
                  Z = Zloc
                }
              }
            }
            list_i = list_i[-1]
            if (length(list_i) > 0) {
              # s'il reste des candidats
              if (random) {
                # si alea, on regarde si ca vaut la peine de continuer
                go = rbinom(1, 1, max(W[list_i, jloc], 0))#on regarde si on continue en fonction du meilleur restant
              } else{
                # si pas d'alea on continue tant qu'il y a des candidats
                go = 1
              }
            } else{
              go = 0
            }
          } else{
            go = 0
            list_i = c() # la contrainte de complexite est definitivement foutue donc on arrete
          }
        }
      }
    }
  } else{
    # multinom
    for (j in 1:p) {
      jloc = list_j[j]
      # on teste la colonne en fonction de sa plus grande proba
      if (random) {
        go = rbinom(1, 1, max(W[, jloc]))
      } else{
        go = 1
      }
      if (go == 1) {
        # si la colonne est jugee interessante
        i = 1
        list_i = 1:p
        list_i = list_i[-jloc]#on enleve la diagonale
        while (go & length(list_i) > 0) {
          wloc = W[list_i, jloc]
          wloc = wloc / wloc
          qui = which(rmultinom(1, 1, wloc) == 1)
          iloc = list_i[qui] # on choisit selon une multinomiale ponderee par les poids normalises
          if (sum(Z[, jloc]) < p1max) {
            # si on respecte la contrainte de complexite
            croisement = FALSE
            if (sum(Z[jloc, ]) != 0 | sum(Z[, iloc]) != 0) {
              croisement = TRUE
            }
            if (croisement == FALSE |
                relax == TRUE) {
              # si pas croisement ou relax
              Zloc = Z
              Zloc[iloc, jloc] = 1
              if (croisement == TRUE) {
                # donc croisement et relax
                # faire la relax
                Zloc[jloc, ] = 0
                Zloc[, iloc] = 0
              }
              if (length(which(colSums(Zloc) != 0)) <= p2max) {
                if (BIC) {
                  # si on veut s'appuyer sur le bic
                  
                  BICloc_vect = BicZ(X = X, Z = Zloc, Bic_null_vect = Bic_null_vect,  Bic_old = BIC_opt_vect, Zold = Z)
                  BICloc = sum(BICloc_vect)
                  
                  if (BICloc < BIC_opt) {
                    # modification effective
                    Z = Zloc
                    BIC_opt_vect = BICloc_vect
                    BIC_opt = BICloc
                  }
                } else{
                  # on se fout du bic
                  Z = Zloc
                }
              }
            }
            list_i = list_i[-qui]
            if (length(list_i) > 0) {
              # s'il reste des candidats
              if (random) {
                # si alea, on regarde si ca vaut la peine de continuer
                go = rbinom(1, 1, max(W[list_i, jloc], 0)) # on regarde si on continue en fonction du meilleur restant
              } else{
                # si pas d'alea on continue tant qu'il y a des candidats
                go = 1
              }
            } else{
              go = 0
            }
          } else{
            go = 0
            list_i = c() # la contrainte de complexite est definitivement foutue donc on arrete
          }
        }
        
      }
    }
    
  }
  return(Z)
}