R/searchZmiss.R

In CorReg: Linear Regression Based on Linear Structure Between Variables

# ' MCMC
searchZmiss <- function(X = X, ini = NULL, maxit = 10^5, M = NULL, plot = FALSE, BIC_vrai = NULL, Z_vrai = NULL, mixmod = mixmod, nbclustmax = 10,
                        BIC_vide_vect = NULL, nett = TRUE, mode = c("relax", "rejet"), MH = TRUE, p2max = Inf, rmax = 5, nbit = 1) {
  # X est la matrice X sans la constante
  # Z est la matrice Z du formalisme priv?e de sa premi?re ligne (constante sur 1) et de sa premi?re colonne (identit?)
  # MH dit si on fait le saut vers le meilleur strict ou si on se contente de l'enregistrer
  mode = match.arg(mode)
  X = 1 * as.matrix(X)
  p = ncol(X)
  n = nrow(X)
  if (is.null(M)) {
    M = 0 * X
    M[is.na(M)] = 1
    X[is.na(X)] = 0
  }
  if (p2max < 0) {
    p2max = Inf
  }
  if (rmax < 0) {
    rmax = Inf
  }
  if (is.null(ini)) {
    Z = matrix(0, ncol = p, nrow = p)
  } else {
    Z = matrix(ini)
  }
  Z_opt = Z
  if (is.null(mixmod)) {
    mixmod = density_estimation(X = X, nbclustmax = nbclustmax, detailed = TRUE)
    mixmod = mixmod_adapter(mixmod)
    BIC_vide_vect = mixmod$BIC_vect
  }
  BIC_vide = sum(BIC_vide_vect)
  BIC_vect = BICZmiss(X = X, Z = Z_opt, Bic_vide_vect = BIC_vide_vect, intercept = TRUE, mixmod = mixmod, nbit = nbit)
  BIC_opt = sum(BIC_vect)
  message(paste("BIC ini:", BIC_opt))
  
  # on initialise
  BIC = BIC_opt
  etape = 0
  pas = 0
  etape_max = round(maxit / (2 * p - 2))
  if (plot) {
    compare = compare_struct(trueZ = Z_vrai, Zalgo = Z)
    courbe_BIC = rep(0, times = round(etape_max))
    courbe_BIC[1] = BIC
    courbe_BIC_opt = rep(0, times = round(etape_max))
    courbe_BIC_opt[1] = BIC_opt
    courbe_compl = rep(0, times = round(etape_max))
    courbe_compl[1] = sum(Z)
    courbe_nbbon1 = rep(0, times = round(etape_max))
    courbe_nbbon1[1] = compare$nbbon1
    courbe_nbtrop = rep(0, times = round(etape_max))
    courbe_nbtrop[1] = compare$nbtrop
    courbe_nbmank = rep(0, times = round(etape_max))
    courbe_nbmank[1] = compare$nbmank
    courbe_bon_gauche = rep(0, times = round(etape_max))
    courbe_bon_gauche[1] = compare$bon_gauche
    courbe_faux_gauche = rep(0, times = round(etape_max))
    courbe_faux_gauche[1] = compare$faux_gauche
  }
  ptm <- proc.time()
  candidats = (1:p^2)[-seq(1, p^2, p + 1)]
  
  nb_cand_loc = 2 * p - 2
  BIC_loc_mat = matrix(0, ncol = 3, nrow = (nb_cand_loc + 1))
  meilleur_strict = FALSE
  consecutif = TRUE
  nbstat = 0
  courbes = matrix()
  while (etape < maxit) {
    ploc = sample(1:p, 1)
    list_cand = rbind(cbind((1:p)[-ploc], ploc), cbind(ploc, (1:p)[-ploc]))
    BIC_loc_mat[, c(2, 3)] = 0
    BIC_loc_mat[, 1] = Inf
    BIC_loc_mat[1, ] = c(BIC, 0, 0)
    if (mode == "rejet") {
      for (candidat in 1:nb_cand_loc) {
        Zloc = Z
        etape = etape + 1
        i = list_cand[candidat, 1]
        j = list_cand[candidat, 2]
        Zloc[i, j] = 1 - Zloc[i, j] # on fait la modif
        if ((Zloc[i, j] & sum(Zloc %*% Zloc) != 0) | (length(which(colSums(Zloc) != 0)) > p2max) | (colSums(Zloc)[j] > rmax)) {next} # croisement ou trop de ssreg ou ssreg trop compl, on passe
        #         BIC_loc=sum(calcul_BIC2.2(Zloc=Zloc,X_appr=X,BIC_ini=BIC_vide_vect,BIC_Z=BIC_vect,Z=Z))
        BIC_loc = sum(BICZmiss(X = X, Z = Zloc, Bic_vide_vect = BIC_vide_vect, intercept = TRUE, mixmod = mixmod, nbit = nbit, BicOld = BIC_vect, Zold = Z))
        BIC_loc_mat[candidat + 1, ] = c(BIC_loc, i, j)
        #         if(is.na(BIC_loc)){message(paste("Perfect sub-regression for variable",which(is.na(calcul_BIC2.2(Zloc=Zloc,X_appr=X,BIC_ini=BIC_vide_vect,BIC_Z=BIC_vect,Z=Z)))))}
      }
    } else if (mode == "relax") {
      for (candidat in 1:nb_cand_loc) {
        Zloc = Z
        etape = etape + 1
        i = list_cand[candidat, 1]
        j = list_cand[candidat, 2]
        Zloc[i, j] = 1 - Zloc[i, j]
        if (as.numeric(Z[i, j]) == 1) { # on supprime le croisement
          Zloc[, i] = 0
          Zloc[j, ] = 0
        }
        if ((length(which(colSums(Zloc) != 0)) > p2max) | (colSums(Zloc)[j] > rmax)) {next} # si apr?s relax on est trop complexe, on passe
        #         BIC_loc=sum(calcul_BIC2.2(Zloc=Zloc,X_appr=X,BIC_ini=BIC_vide_vect,BIC_Z=BIC_vect,Z=Z))
        BIC_loc = sum(BICZmiss(X = X, Z = Zloc, Bic_vide_vect = BIC_vide_vect, intercept = TRUE, mixmod = mixmod, nbit = nbit, BicOld = BIC_vect, Zold = Z))
        BIC_loc_mat[candidat + 1, ] = c(BIC_loc, i, j)
        if (is.na(BIC_loc)) {message(paste("Perfect sub-regression for variable", which(is.na(BICZmiss(X = X, Z = Zloc, Bic_vide_vect = BIC_vide_vect, intercept = TRUE, mixmod = mixmod, nbit = nbit, BicOld = BIC_vect, Zold = Z)))))}
      }
    }
    # on termine l'?tape
    pas = pas + 1
    # on tire au sort entre stationnarit? et action
    loc = BIC_loc_mat[, 1]
    if (min(loc[loc != 0]) < BIC_opt) { # on stocke, on pr?vient, et on y va
      # qui=which(BIC_loc_mat[1:length(loc),1]==min(BIC_loc_mat[1:length(loc),1]))[1]
      meilleur_strict = TRUE
      Zloc = Z
      qui = which.min(loc)
      BIC_opt = BIC_loc_mat[qui, 1]
      i = BIC_loc_mat[qui, 2]
      j = BIC_loc_mat[qui, 3]
      if (mode == "rejet") { # donc on a pas eu de croisement
        Zloc[i, j] = 1 - Zloc[i, j] # on fait la modif
      } else if (mode == "relax") {
        Zloc[i, j] = 1 - Zloc[i, j]
        if (as.numeric(Z[i, j]) == 1) { # on supprime le croisement
          Zloc[, i] = 0
          Zloc[j, ] = 0
        }
      }
      Z_opt = Zloc
      if (MH & meilleur_strict) { # on saute
        meilleur_strict = F
        Z = Z_opt
        BIC = BIC_opt
        # cat(paste("Best BIC",BIC))
      }
    } else { # si pas de meilleur, on tire au sort
      meilleur_strict = F # normalement inutile
      loc2 = loc
      loc = loc[loc != Inf]
      # on doit maintenant ?tablir la correspondance
      loc = loc - min(loc) # rend tout positif (deltabic) . Le min (ancien max) devient nul
      loc = exp(-0.5 * loc) / sum(exp(-0.5 * loc))
      loc = cumsum(loc) / sum(loc) # on cr?e les intervalles associ?s pour le tirage au sort
      station = which(loc > runif(1))[1]
      loc2[loc2 != Inf] = 1:length(loc)
      station = which(loc2 == station) # on va cherche l'indice correspondant
      if (length(station) == 0) {
        message("No feasible candidate")
        station = 1
      }
      if (station != 1) { # alors pas de stationnarit?
        BIC = BIC_loc_mat[station, 1]
        i = BIC_loc_mat[station, 2]
        j = BIC_loc_mat[station, 3]
        if (mode == "rejet") { # donc on a pas eu de croisement
          Z[i, j] = 1 - Z[i, j] # on fait la modif
        } else if (mode == "relax") {
          Z[i, j] = 1 - Z[i, j]
          if (as.numeric(Z[i, j]) == 1) { # on supprime le croisement
            Z[, i] = 0
            Z[j, ] = 0
          }
        }
      } else {nbstat = nbstat + 1}
    }
    BIC_vect = BICZmiss(X = X, Z = Zloc, Bic_vide_vect = BIC_vide_vect, intercept = TRUE, mixmod = mixmod, nbit = nbit) # completement sous-optimal
    
    #     BIC_vect=calcul_BIC2.2(Zloc=Z,X_appr=X,BIC_ini=BIC_vide_vect)#voir si pas pr?f?rable de stocker l'etape precedente pour l'utiliser ici
    if (plot) {
      compare = compare_struct(trueZ = Z_vrai, Zalgo = Z)
      courbe_BIC[pas] = BIC
      courbe_BIC_opt[pas] = BIC_opt
      courbe_compl[pas] = sum(Z)
      courbe_nbbon1[pas] = compare$nbbon1
      courbe_nbtrop[pas] = compare$nbtrop
      courbe_nbmank[pas] = compare$nbmank
      courbe_bon_gauche[pas] = compare$bon_gauche
      courbe_faux_gauche[pas] = compare$faux_gauche
    }
  }
  message(proc.time() - ptm)
  if (nett) {
    
  }
  if (plot) {
    courbes = rbind(courbe_BIC, courbe_BIC_opt, courbe_compl, courbe_nbbon1, courbe_nbtrop, courbe_nbmank, courbe_bon_gauche, courbe_faux_gauche)
  }
  return(list(Z = Z_opt, BIC = BIC_opt, BIC_vide_vect = BIC_vide_vect, Zloc = Z, courbes = courbes, ptm = (proc.time() - ptm)[1], nbstat = nbstat))
}