R语言互信息和条件互信息

墨蓝 2023-07-04 04:40 20阅读 0赞

### 互信息和条件互信息 ###

最近看了信息熵的相关知识，考虑用信息熵来计算互信息和条件互信息。

MI(X,Y)=H(X)-H(X|Y)

H(X):熵，一种不确定性的度量 H(X,Y):联合熵，两个元素同时发生的不确定度

MI(X,Y)即是信息增益：

H(X,Y)=H(X)+H(Y|X)->H(X|Y)=H(X,Y)-H(Y)

则MI(X,Y)=H(X)+H(Y)-H(X,Y)

MI(X1，X2|Y)=H(X1|Y)-H(X1|X2，Y)=H(X1，Y)-H(Y)-(H(X1，X2，Y)-H(X2，Y)) =H(X1，Y)+H(X2，Y)-H(Y)-H(X1，X2，Y)

正好在R语言当中有直接计算信息熵的包，代码如下：

library(infotheo)
    data<-read.csv("test.csv",header = F)
    
    #互信息
    mi<-function(mylist1,mylist2){
      return(entropy(mylist1)+entropy(mylist2)-entropy(cbind(mylist1,mylist2)))
    }
    pvalue<-mi(data[,1],data[,2])#测试
    
    
    #条件互信息
    cmi<-function(mylist1,mylist2,mylist3){
      return(entropy(cbind(mylist1,mylist3))+entropy(cbind(mylist2,mylist3))-entropy(mylist3)-entropy(cbind(mylist1,mylist2,mylist3)))
    }
    pvalue1<-cmi(data[,5],data[,13],data[,c(7,8)])

## 案例结果 ##

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70]

setwd("C:\\Users\\Administrator\\Desktop\\MI_cntry 3")
    par(mfrow = c(2, 2), mar = c(4, 4, 1, 1))  ###图形分栏布局
    #---------------------------------------------------------GDP
    data1 <- read.csv("GDP.csv")
    str(data1)
    zdata1<- t(data1)
    mdata1 <- zdata1[-c(1,2),]
    m=seq(from = 2001, to = 2015, by=1)  ## 先生成等差数列的刻度，间距为2.5
    x1=seq(from =1, to = 15, by=3)
    x2=seq(from = 2001, to = 2015, by=3)
    plot(mdata1[,1],type='l',main="CHN",xlab="",ylab="GDP",xaxt="n")
    axis(side=1,at=x1,labels=x2)
    #mtext("*10^13",adj=0,font=0.0001,font.axis=6)
    
    #----------------------------------------------------------NP
    data2 <- read.csv("NP.csv")
    str(data2)
    zdata2<- t(data2)
    mdata2 <- zdata2[-c(1,2),]
    m=seq(from = 2001, to = 2015, by=1)  ## 先生成等差数列的刻度，间距为2.5
    x1=seq(from =1, to = 15, by=3)
    x2=seq(from = 2001, to = 2015, by=3)
    plot(mdata2[,1],type='l',main="CHN",xlab="",ylab="NP",xaxt="n")
    axis(side=1,at=x1,labels=x2)
    #mtext("*10^13",adj=0,font=0.0001,font.axis=6)
    
    #----------------------------------------------------------U
    
    data3 <- read.csv("U.csv")
    str(data3)
    zdata3<- t(data3)
    mdata3 <- zdata3[-c(1,2),]
    m=seq(from = 2001, to = 2015, by=1)  ## 先生成等差数列的刻度，间距为2.5
    x1=seq(from =1, to = 15, by=3)
    x2=seq(from = 2001, to = 2015, by=3)
    plot(mdata3[,1],type='l',main="CHN",xlab="",ylab="U",xaxt="n")
    axis(side=1,at=x1,labels=x2)
    #mtext("*10^13",adj=0,font=0.0001,font.axis=6)
    
    
    
    #----------------------------------------------------------
    ############互信息值计算
    #install.packages("infotheo")
    library(infotheo)
    #条件互信息
    cmi<-function(mylist1,mylist2,mylist3){
    return(entropy(cbind(mylist1,mylist3))+entropy(cbind(mylist2,mylist3))-entropy(mylist3)-entropy(cbind(mylist1,mylist2,mylist3)))
     }
    
    #pvalue1<-cmi(mdata1[,1],mdata2[,1],mdata3[,1])
    
    #plot(pvalue1)
    
    #互信息参考网站：https://www.cnblogs.com/jiangyaling/p/8040024.html
    
    mi<-function(mylist1,mylist2){
    return(entropy(mylist1)+entropy(mylist2)-entropy(cbind(mylist1,mylist2)))
    }
    pvalue01<-mi(mdata1[,2],mdata2[,2])#测试
    pvalue02<-mi(mdata2[,2],mdata3[,2])#测试
    
    plot(pvalue01,ylin=c(2,3),yaxt="n",ylab="",pch=3)
    par(new=TRUE)
    plot(pvalue02,ylim=c(2,3),pch=6)

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70]: https://img-blog.csdnimg.cn/20200207220137762.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu,size_16,color_FFFFFF,t_70