本文整理汇总了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)
}
//.........这里部分代码省略.........