C++ CData::get_logUnif_y_tilde方法代码示例

void CParam::init_logUnif_y_tilde(CData &Data,CFeasibilityMap &FM, Uniform &randUnif, int n_simul) {
  logUnif_y_tilde = ColumnVector(Data.n_faulty);
  for (int i = 1; i <= Data.n_faulty; i++) {
    ColumnVector s_i = S_Mat.column(i);
    int i_original = Data.Faulty2Original[i-1];
    // double log_f_y_tilde_q = 0;  Commented by HANG on 5/16/2015
    if (FM.useMap && Data.is_case(i_original,2)) {
      logUnif_y_tilde(i) = FM.Simulate_logUnif_case2(CFeasibilityMap::s_to_tau_fn(s_i),i_original,n_simul,randUnif);
    } else {
      logUnif_y_tilde(i) = Data.get_logUnif_y_tilde(s_i, CFeasibilityMap::s_to_tau_fn(s_i),i_original); 
    } 
    //logUnif_y_tilde(i) = Data.get_logUnif_y_tilde(s_i, CFeasibilityMap::s_to_tau_fn(s_i),i_original);
  }
}

void CParam::S1(int iter, Uniform &randUnif, CData &Data, CFeasibilityMap &FM, int n_simul) {
  is_accept(1) = 0;
  for (int i_faulty=1; i_faulty<=Data.n_faulty; i_faulty++) {
    if (FM.useMap) {FM.pmm->iter = iter;}
    int i_original = Data.Faulty2Original[i_faulty-1];
    double g_option_q, g_mode_q;
    int what_type_move;
    ColumnVector s_i = S_Mat.column(i_faulty);
    int tau_q = FM.EvaluateMove(i_original, Data, s_i, iter, what_type_move, g_option_q, g_mode_q, true);
    // calculate g_mode_i & g_option_i
    ColumnVector s_q = CFeasibilityMap::tau_to_s_fn(tau_q, n_var);
    double g_option_i, g_mode_i;
    FM.EvaluateMove(i_original, Data, s_q, iter, what_type_move, g_option_i, g_mode_i, false);
    // calculate g_s_q & g_s_i
    double g_s_q = g_mode_q * g_option_q;
    double g_s_i = g_mode_i * g_option_i;

    // whether drawn by item_by_rnorm and balance edits
    ColumnVector item_by_bal;
    ColumnVector item_by_rnorm;
    if (Data.is_case(i_original,1)) {
      item_by_rnorm = Data.get_item_by_norm_indicator(s_q,item_by_bal);
    } else {
      item_by_rnorm = Data.copy_non_balance_edit(s_q);
    }
    
    //Generate from normal distribution
    ColumnVector mu_z_i = Mu.column(z_in(i_original)) ;
    ColumnVector tilde_y_i = Data.log_D_Observed.row(i_original).t();
    ColumnVector y_q(n_var);
    LowerTriangularMatrix LSigma_i;
    double log_cond_norm_q = calculate_log_cond_norm(Data, i_original, item_by_rnorm, tilde_y_i, y_q, true, LSigma_i, s_q); // MODIFIED 2015/02/16
    
    // Put values from balance edits
    ColumnVector x_q = exp_ColumnVector(y_q);
    if (Data.is_case(i_original,1)) {
      Data.set_balance_edit_values_for_x_q(s_q, x_q, item_by_bal); // CHANGED by Hang, 2014/12/29
    } else {
      Data.update_full_x_for_balance_edit(x_q);
    }

    // Acceptance/Rejection
    if (Data.PassEdits(x_q)) {  // Check constraints
      ColumnVector item_i_by_rnorm;
      if (Data.is_case(i_original,1)) {
        item_i_by_rnorm = Data.get_item_by_norm_indicator(s_i,item_by_bal);
      } else {
        item_i_by_rnorm = Data.copy_non_balance_edit(s_i);
      }
    
    double log_cond_norm_i = calculate_log_cond_norm(Data, i_original, item_i_by_rnorm, tilde_y_i, y_q, false, LSigma_i, s_i); // MODIFIED 2015/02/16

      y_q = log_ColumnVector(x_q) ;
      ColumnVector y_i = (Y_in.row(i_original)).t();
  
      // Calculate acceptance ratio
      ColumnVector y_compact_q = Data.get_compact_vector(y_q);
      ColumnVector y_compact_i = Data.get_compact_vector(y_i);
      double log_full_norm_q = log_MVN_fn(y_compact_q,mu_z_i,LSIGMA_i[z_in(i_original)-1],logdet_and_more(z_in(i_original)));
      double log_full_norm_i = log_MVN_fn(y_compact_i,mu_z_i,LSIGMA_i[z_in(i_original)-1],logdet_and_more(z_in(i_original)));
      
    
      double log_f_y_tilde_q = 0;
      if (FM.useMap && Data.is_case(i_original,2)) {
        log_f_y_tilde_q = FM.Simulate_logUnif_case2(CFeasibilityMap::s_to_tau_fn(s_q),i_original,n_simul,randUnif);
      } else {
        log_f_y_tilde_q = Data.get_logUnif_y_tilde(s_q, CFeasibilityMap::s_to_tau_fn(s_q),i_original); 
      }
      double logNum = log_full_norm_q + log_f_y_tilde_q - log_cond_norm_q - log(g_s_q);
      double logDen = log_full_norm_i + logUnif_y_tilde(i_faulty) - log_cond_norm_i - log(g_s_i);
    
      accept_rate(1) = exp( logNum - logDen ) ;
      if ( randUnif.Next() < accept_rate(1) ){
        Y_in.row(i_original) = y_q.t() ;
        S_Mat.column(i_faulty) = s_q;
        logUnif_y_tilde(i_faulty) = log_f_y_tilde_q;
        is_accept(1)++;
      }
    }
    
  }
  is_accept(1) = is_accept(1) / Data.n_faulty;
}