互信息条件互信息的计算

未分类

公式

互信息计算

协方差矩阵计算函数：

getDoubleHang_value<-function(mylist1,mylist2){
  Dx1=var(mylist1)
  Dx2=var(mylist2)
  double_h_value=Dx1*Dx2-cov(mylist1,mylist2)*cov(mylist1,mylist2)  // cov 计算C（X）
}

计算整个互信息：

for(xi in 1:nrow(z)){
  mylist_a<-as.double(c(z[xi,]))
  Dxa=var(mylist_a)                                 //计算 行xi 的方差
  for(yi in 1:nrow(z)){
    mylist_b<-as.double(c(z[yi,]))
    Dxb=var(mylist_b)                               //计算 行yi 的方差
    cov_ab=getDoubleHang_value(mylist_a,mylist_b)   //计算 MI(X,Y)
    mi=0.5*log(Dxa*Dxb/cov_ab)                      //套公式计算上图的公式(左边)
    mylist_info<-c(mylist_info,mi)                  //mylist_info 是最后的结果集
  }
}

条件互信息计算

info <- Info(Translat(graE[i]),Translat(graE[i+1]),Translat(brek),gene)

条件互信息的计算函数
输入：gxi源点；gyi入点；cutset割集；gen源数据矩阵
输出：CMI

Info <- function(gxi,gyi,cutset,gen){
  cutset<-data.frame(cutset,stringsAsFactors = F)
  x<-c()
  y<-c()
  cutset_table<-data.frame()
  for(xi in 1:nrow(gen)){
    if(gxi==gen[xi,1]){
      x<-as.numeric(c(gen[xi,2:length(gen)]))
    }
    if(gyi==gen[xi,]){
      y<-as.numeric(gen[xi,2:length(gen)])
    }
    for(yi in 1:nrow(cutset)){
      if(gen[xi,1]==cutset[yi,1]){  
        if(nrow(cutset_table)==0)
          cutset_table<-data.frame(t(gen[xi,]),stringsAsFactors = F)
        else
          cutset_table<-data.frame(cutset_table,t(gen[xi,]),stringsAsFactors = F)
      }
    }
  }
  cutset_table<-data.frame(t(cutset_table),stringsAsFactors = F) 
  cutset_table<-cutset_table[,-1]
  to_list<-c()
  cut<-data.frame()
  for(xi in 1:nrow(cutset_table)){
    to_list<-as.numeric(c(cutset_table[xi,]))
    if(nrow(cut)==0)
      cut<-data.frame(to_list)
    else
      cut<-data.frame(cut,to_list)
  }
  // 上边计算完割集cut
  cut_x<-data.frame(cut,x)
  cut_y<-data.frame(cut,y)
  cut_x_y<-data.frame(cut,x,y)
  t=det(cov(cut_x))*det(cov(cut_y))/(det(cov(cut))*det(cov(cut_x_y))) // 套公式（图右）
  cmi=0.5*log(t)
}

调用阶段

for(i in  seq(1,length(graE),2)){
      g_d <- g - edge(paste(graE[i],"|",graE[i+1],sep = ""))  
      if( edge_connectivity(g_d, source = graE[i], target = graE[i+1], checks = TRUE) > 0){
            shortpa<-shortest_paths(g_d, from = graE[i], to = graE[i+1], mode = c("all"))$vpath
            one_path<-names(V(g))[as.integer(shortpa[[1]])]
            brek <- one_path[-c(1,length(one_path))]
            if(length(brek) > 0){
                  info <- Info(Translat(graE[i]),Translat(graE[i+1]),Translat(brek),gene) 
                  // 条件互信息CMI 计算
                  if(info < weight_3){ g <- g_d }
            }else{
                  print(i)
            }
      }
}

互信息计算的替换

# 原先的核心逻辑
for(xi in 1:nrow(z)){
  mylist_a<-as.double(c(z[xi,]))
  Dxa=var(mylist_a)
  for(yi in 1:nrow(z)){
    print(yi) 
    #每次计算瞬间完成，主要归功于as.double(c(z[xi,]))，整体计算
    mylist_b<-as.double(c(z[yi,]))
    Dxb=var(mylist_b)
    cov_ab=getDoubleHang_value(mylist_a,mylist_b)
    mi=0.5*log(Dxa*Dxb/cov_ab)
    mylist_info<-c(mylist_info,mi)
  }
}
# 现阶段是计算的元素，所以非常的慢
for(xi in 1:ncol(data)){
  mylist_a<-data[,xi]
  for(yi in 1:ncol(data)){
    mylist_b<-data[,yi]
    mi <- myfunction3(mylist_a,mylist_b)
    mylist_info<-c(mylist_info,mi)
  }
}

区别：
原先快，现在慢
原先的结果矩阵对角为Inf，现在的对角有值，即有（本身，本身）的互信息
原先的结果34/28正确，现在只替换互信息计算后45/18正确

条件互信息计算的替换

发生在三阶段的第二三阶段，原先都是先计算割集，然后将割集作为计算条件
现在直接更换公式，不计算割集

110/34正确