机器学习——建立因果连系（传递熵）

因果联系的建立（Causality Inference)

学习目标：

1·了解什么是信息熵；
2·如何建立因果联系；
3·代码coding.

信息熵

由信息论之父——克劳德·艾尔伍德·香农提出，并首次用数学公式阐明了概率与信息冗余度之间的关系。为了方便理解，我们可以参照热力学里面的熵来理解信息熵。在热力学中，我们用熵来度量体系混乱程度，熵增代表物质向无规则方向发展即体系越来越混乱。其次，我们可以通过降温等手段来减熵即恢复有序性。同样，在信息论中，我们也可以这样理解：信息熵是用来度量体系不确定度的量，被观测的系统越不确定信息熵越大，系统越稳定信息熵越小。其次，我们可以通过被告知某个信息来降低系统的不确定性即‘减熵’。
但是值得注意的一点是，这个信息不一定可以降低不确定度，有时可能不会给你提供任何有效信息甚至可能会增大这个系统的不确定度。打个比方：玩投色子的游戏，点数1~6，投完色子让你去猜点数。在你猜之前呢，我告诉你的信息“这个点数是个位数”。那么对你来说这条信息就是无用的，他不会提高你猜对的概率。
下面是Shannon Entropy的公式：

我们可以看到信息熵是度量单个变量的不确定度的，我们在之后的运用就是建立在此基础上的。

因果关系

就是字面上的意思，获取原因的信息可以降低我们观测结果的不确定性。当我们知道了单个变量的信息熵之后，我们如何去度量变量之间是如何彼此影响的呢？我们应该如何去建立多个变量之间的因果关系呢？这个时候我们就需要引入另一个概念——传递熵(Tansfer Entropy)

传递熵

信息熵是用来研究变量与变量之间的信息的传递，可以计算这个信息传递能减少多少被观测系统的不确定度。当 X 对 Y 的传递熵 > Y 对 X 的传递熵时，我们就把 X 称为因，把 Y 称为果，并以此来建立两个变量之间的因果关系。

为了方便理解，我在这里给出了一张图，用比较形象的方式来阐述：
中间小三角形的部分就代表了传递熵 T(Xt>Yt,t) ，而其它三个圆形就代表了在 t-1 时刻 X 与 Y 信息熵和 t 时刻 Y 的信息熵。首先可以明确的是，Y 的历史信息可以降低 Y 的不确定性，而对于 X 的信息能否降低 Y 的不确定度，我们可以通过计算 Transfer Entropy 的大小来判断。那么怎么计算传递熵呢？
我们给出推导（推导较为复杂这里不一一赘述）：
我们可以把三个圆形的面积当作信息熵，通过对圆形面积做加减法最后算出中间那个“三角形”，如下：

以下是对用到的四个信息熵分别做计算：


最后带入、化简得出公式如下：

我们可以看到（1）式只有变量之间的联合概率分布P，所以在写代码的时候我们运用（1）式去编码是最方便的。
到此我们已经掌握了传递熵的计算，理论知识已经准备好了，接下来就是代码方面的补充：

气候科学实战

我们今天还是去研究：降水与气温之间的关系。：）

#读入数据
file1 = 'pr_Amon_GFDL-ESM2M_historical_r1i1p1_186101-200512.nc'
file2 = 'tas_Amon_GFDL-ESM2M_historical_r1i1p1_186101-200512.nc'

ds1 = Dataset(file1,'r')
ds2 = Dataset(file2,'r')

pr = ds1.variables['pr'][:,:,:]
tas = ds2.variables['tas'][:,:,:]
lon = ds1.variables['lon'][:]
lat = ds1.variables['lat'][:]

全球范围太大了，我们先以 New York 为例：

# site New York
# 40.7, -74
lon_NY_index = int((360-74) / (360/144))
lat_NY_index = int( 40.7 / (180/90) + 45)

tas_NY = tas[:,lat_NY_index,lon_NY_index]
pr_NY = pr[:,lat_NY_index,lon_NY_index]

计算联合概率分布：

# calculate transfer entropy
#首先确定变量——温度的历史、降水的历史和降水
Data = np.hstack((tas_NY[:-1].reshape(-1,1), pr_NY[1:].reshape(-1,1), pr_NY[:-1].reshape(-1,1)))
#为了计算概率，我们需要把这三条时间序列数据分别以相应的min和max为界分成 x（随你选，不要太大也不要太小）个box
#去看一下分别有多少个数值掉进相应的box这样就可以计算概率了。
#这里我们分了11个box 
nBins = 11
#nSignals代表有几条时间序列
[nData, nSignals] = np.shape(Data)
#找出每条时间序列的最大值和最小值
# classify data
binEdges = np.full([nSignals, nBins], np.nan)
minEdge = np.full([nSignals],np.nan)
maxEdge = np.full([nSignals],np.nan)
for s in range(nSignals):
    minEdge[s] = np.min(Data[:,s])
    maxEdge[s] = np.max(Data[:,s])
    binEdges_all = np.linspace(minEdge[s], maxEdge[s], num=nBins+1)   #等差分割
    binEdges[s,0:nBins] = binEdges_all[1:]
# 切分box
classifiedData = np.full([nData,nSignals], np.nan)
for s in range(nSignals):
    # 保护原始数据
    smat = Data[:,s].copy()
    # 设定分类矩阵
    cmat = np.full([nData], np.nan)
    # loop over local bins
    for e in range(nBins):
        if np.where( smat<= binEdges[s,e])[0].size > 0:
            cmat[ np.where( smat<= binEdges[s,e])[0] ] = e
            smat[ np.where( smat<= binEdges[s,e])[0] ] = maxEdge[s] + 9999   #为了防止前面的数被重复覆盖
    classifiedData[:,s] = cmat
    
#统计有多少个数值掉入了相应的box
C = np.full([nBins, nBins, nBins], 0)
for i in range(nData):
    C[ int(classifiedData[i,0]), int(classifiedData[i,1]),  int(classifiedData[i,2]) ] = C[ int(classifiedData[i,0]), int(classifiedData[i,1]), int(classifiedData[i,2]) ] + 1

#计算联合概率分布
pXtYwYf = (C+1e-20)/np.sum(C)   #sum函数可以对相应的变量求和，也即积分
# Marginal PDFs
pYw=np.sum(np.sum(pXtYwYf,axis=0),axis=1)
# Joint PDFs
pXtYw=np.sum(pXtYwYf,axis=2)
pYwYf=np.sum(pXtYwYf,axis=0)

最后可以计算出温度对降水的传递熵：

# transfer Shannon information entropy 
a = pXtYwYf * pYw
b = pXtYw * pYwYf
H = np.sum(pXtYwYf * np.log2(a/b))
print(H)                 #温度 -> 降水  H=1.6054489

至此我们发现已知温度信息可以降低降水1.6的信息熵，由此也就可以建立一个最简单的因果关系啦！
实际上建立因果关系要远比这个例子要复杂得多。因为我们只计算了温度对降水的传递熵，我们还可以去计算一下降水对温度的传递熵，比较一下哪个因素才是在这个系统里占据主导地位。当然我们还可以把很多的变量放到一个系统里去看看他们之间的因果关系。