当前位置: 首页>>代码示例>>C++>>正文


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

本文整理汇总了C++中CData::E_forcein_mat方法的典型用法代码示例。如果您正苦于以下问题:C++ CData::E_forcein_mat方法的具体用法?C++ CData::E_forcein_mat怎么用?C++ CData::E_forcein_mat使用的例子?那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在CData的用法示例。


在下文中一共展示了CData::E_forcein_mat方法的1个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。

示例1: log

void CParam::S6_G_beta_ij(CData &Data) {
  
  arma::mat G_prop = G_mat ; arma::mat Beta_prop = Beta ; 
  double logQ_numer, logQ_denom ; 
  double logP_numer, logP_denom ;
  double normC_prop ;
  
  bool is_forcein_edge_selected = false ;
  
	// double P_G_q, P_G ; // For updated E(i,j)
	
  // Step 1
  int w_cur = 0 ; int w_max = 0 ; int w_prop ; 
  for (int i=0; i<(n_pheno-1); i++){
    for (int j=(i+1); j<n_pheno; j++){
      w_max ++ ; w_cur = w_cur + G_mat(i,j) ; 
    }
  }
  if ( w_cur==0 ){ w_prop = 1 ; logQ_numer = log(0.5) ; logQ_denom = log(1.0) ; } // q(w|w^q) // q(w^q|w)
  if ( w_cur==w_max ){ w_prop = w_max - 1 ; logQ_numer = log(0.5) ;  logQ_denom = log(1.0) ; } // q(w|w^q) // q(w^q|w)
  if ( (w_cur > 0) && (w_cur < w_max) ){
    RandVec = Rcpp::runif(1, 0, 1) ; 
    if ( RandVec(0) < 0.5 ){ w_prop = w_cur + 1 ; } else { w_prop = w_cur - 1 ; }
    logQ_denom = log(0.5) ; // q(w^q|w)
    if ( (w_prop==0) || (w_prop==w_max) ){ logQ_numer = log(1.0) ; } else { logQ_numer = log(0.5) ; } // q(w|w^q)
  }
  
  // Step 2
  if ( w_prop > w_cur  ){ // step 2-a
    
		int id_added = rDiscrete(w_max-w_cur) + 1 ; // 1 ~ total. of empty edges
    int count_empty_edge = 0 ;
		
    for (int i=0; i<(n_pheno-1); i++){
      for (int j=(i+1); j<n_pheno; j++){
        if ( G_mat(i,j)==0 ){
          count_empty_edge++;
          if ( id_added==count_empty_edge ){
            G_prop(i,j) = 1 ; G_prop(j,i) = G_prop(i,j) ; 
            RandVec = Rcpp::rgamma(1, Data.a_betaG, 1.0/Data.b_betaG) ; // q(beta_ij^q)
            Beta_prop(i,j) = RandVec(0) ; 
            Beta_prop(j,i) = Beta_prop(i,j) ; 
            id_added = -9 ; 
            
						// if ( Data.PriorSetting==2 ){
						//	P_G = 1.0 - Data.priorprob_G(i,j) ; // Bernoulli(E_ij=0; p_ij)
						//	P_G_q = Data.priorprob_G(i,j) ; // Bernoulli(E_ij^q=1; p_ij)
						// }
												
            logP_numer = R::dgamma(Beta_prop(i,j),Data.a_beta,(1.0/Data.b_beta),1) ; // f(beta_ij^q|G_ij^q)
            logP_denom = 0 ; // f(beta_ij|G_ij) cancelled with q(beta_ij)
            normC_prop = normC_fn(Beta_prop, Data) ;
            for (int i_SNP=0; i_SNP<n_SNP; i_SNP++){
              double e_it = E_mat(i,i_SNP) ;
              double e_jt = E_mat(j,i_SNP) ; 
              logP_numer = logP_numer - log(normC_prop) + Beta_prop(i,j) * e_it * e_jt ; // f(e_t|alpha,beta^q,G^q)
              logP_denom = logP_denom - log(normC) + Beta(i,j) * e_it * e_jt ; // f(e_t|alpha,beta,G)
            }
            logQ_numer = logQ_numer + 0 ; // q(beta_ij) cancelled with f(beta_ij|G_ij)  
            logQ_denom = logQ_denom + R::dgamma(Beta_prop(i,j),Data.a_betaG,(1.0/Data.b_betaG),1) ; // q(beta_ij^q)
          } // if ( id_added==count_empty_edge ) 
        } // for ( G_mat(i,j)==0 )
      } // for (j)
    } // for (i)
		
    logQ_numer = logQ_numer - log(w_cur) ; // q(G|G^q,w)
    logQ_denom = logQ_denom - log(w_max-w_cur)  ; // q(G^q|G,w^q)
  
	} else { // step 2-b
	
    int id_deleted = rDiscrete(w_cur)+1 ; // 1 ~ total no. of connected edges
    int count_connected_edge = 0 ; 
    for (int i=0; i<(n_pheno-1); i++){
      for (int j=(i+1); j<n_pheno; j++){
        if ( G_mat(i,j)==1 ){
          count_connected_edge++;
          if ( id_deleted==count_connected_edge ){
            G_prop(i,j) = 0 ; G_prop(j,i) = G_prop(i,j) ; 
            Beta_prop(i,j) = 0 ; // q(beta_ij^q)
            Beta_prop(j,i) = Beta_prop(i,j) ; 
            id_deleted = -9 ; 
            
            if (Data.isforcein==true){
              if (Data.E_forcein_mat(i,j)==1) is_forcein_edge_selected = true ; 
            }
            
						// if ( Data.PriorSetting==2 ){
						//	P_G = Data.priorprob_G(i,j) ; // Bernoulli(E_ij=1; p_ij)
						//	P_G_q = 1.0 - Data.priorprob_G(i,j) ; // Bernoulli(E_ij^q=0; p_ij)
						// }
						
            logP_numer = 0 ; // f(beta_ij^q|G_ij^q) cancelled with q(beta_ij^q)
            logP_denom = R::dgamma(Beta(i,j),Data.a_beta,(1.0/Data.b_beta),1) ; // f(beta_ij|G_ij)
            normC_prop = normC_fn(Beta_prop, Data) ;
            for (int i_SNP=0; i_SNP<n_SNP; i_SNP++){
              double e_it = E_mat(i,i_SNP) ;
              double e_jt = E_mat(j,i_SNP) ; 
              logP_numer = logP_numer - log(normC_prop) + Beta_prop(i,j) * e_it * e_jt ; // f(e_t|alpha,beta^q,G^q)
              logP_denom = logP_denom - log(normC) + Beta(i,j) * e_it * e_jt ; // f(e_t|alpha,beta,G)
            }
//.........这里部分代码省略.........
开发者ID:dongjunchung,项目名称:GGPA,代码行数:101,代码来源:3_Param.cpp


注:本文中的CData::E_forcein_mat方法示例由纯净天空整理自Github/MSDocs等开源代码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。