R的极客理想系列文章
,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。
R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。
要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。
关于作者:
-
张丹(Conan), 程序员/Quant: Java,R,Nodejs
-
blog: http://blog.fens.me
-
email: bsspirit@gmail.com
转载请注明出处:
http://blog.fens.me/r-entropy
香农的《通信的数学理论》是20世纪非常伟大的著作,被认为是现代信息论研究的开端。信息论定义了信息熵,用于把信息进行度量,以比特(bit)作为量纲单位,为如今发达的信息产业和互联网产业奠定了基础。本文接上一篇文章
R语言实现46种距离算法
,继续philentropy包的介绍,包括信息度量函数的使用。
-
信息熵介绍
-
关键概念
-
信息度量函数
-
应用举例
1.信息熵介绍
信息论(Information Theory)是概率论与数理统计的一个分枝,用于研究信息处理、信息熵、通信系统、数据传输、率失真理论、密码学、信噪比、数据压缩等问题的应用数学学科。信息论将信息的传递作为一种统计现象来考虑,给出了估算通信信道容量的方法。信息传输和信息压缩是信息论研究中的两大领域。
香农被称为是“信息论之父”,香农于1948年10月发表的A Mathematical Theory of Communication,
通信的数学理论(中文版)
,通常被认为是现代信息论研究的开端。
信息熵,是对信息随机性的量度,又指信息能被压缩的极限,用bit作为衡量信息的最小单位。一切信息所包含的信息量,都是1bit的正整数倍。计算机系统中常采用二进制编码,一个0或1就是1bit。
举例来说明一下信息熵的计算原理,假设小明最喜欢5种水果,苹果、香蕉、西瓜、草莓、樱桃中的一种,如果小明没有偏爱,选择每种水果的概率都是20%,那么这一信息的信息熵为
H(A) = -1*(0.2*log2(0.2)*5)
= 2.321928 bits
如果小明偏爱香蕉,选择这5种水果的概率分别是10%,20%,45%,15%,10%,那么这一信息信息熵为
H(B)=-1*(0.1*log2(0.1)+0.2*log2(0.2)+0.45*log2(0.45)+0.15*log2(0.15)+0.1*log2(0.1))
= 2.057717 bits
从结果得到H(A)大于H(B),信息熵越大表示越不确定。对于B的情况,对某一种水果的偏好,比A增加了确定性的因素,所以H(B)小于H(A)是符合对于信息熵的定义的。
2. 关键概念
我们从一幅图来认识信息熵,图中显示了随机变量X和Y的2个集合,在信息熵的概念里的所有可能逻辑关系。两个圆所包含的面积为
联合熵H(X,Y)
, 左边的整个圆表示X的熵H(X),左边半圆是
条件熵H(X|Y)
。 右边的整个圆表示Y的熵H(Y),右边半圆
条件熵H(Y|X)
,中间交集的部分是
互信息I(X; Y)
。
信息熵(Entropy):
是对信息随机性的量度,用于计算信息能被压缩的极限。对随机变量X,不确定性越大,X的信息熵H(X)也就越大。
公式定义:
H(x)的取值范围,0<=H(x)<=log(n), 其中n是随机变量x取值的种类数。需要注意的是,熵只依赖于随机变量的分布,与随机变量取值无关。
条件熵(Conditional Entropy):
表示两个随机变量X和Y,在已知Y的情况下对随机变量X的不确定性,称之为条件熵H(X|Y),
公式定义:
联合熵(Joint Entropy):
表示为两个随机事件X和Y的熵的并集,联合熵解决将一维随机变量分布推广到多维随机变量分布。
公式定义:
互信息(Mutual Information, 信息增益):
两个随机变量X和Y,Y对X的互信息,为后验概率与先验概率比值的对数,即原始的熵H(X)和已知Y的情况下的条件熵H(X|Y)的比值的对数,信息增益越大表示条件Y对于确定性的贡献越大。互信息,也可以用来衡量相似性。
公式定义:
当MI(X,Y)=0时,表示两个事件X和Y完全不相关。决策树ID3算法就是使用信息增益来划分特征,信息增益大时,说明对数据划分帮助很大,优先选择该特征进行决策树的划分。
信息增益比率:
是信息增益与该特征的信息熵之比,用于解决信息增益对多维度特征的选择,决策树C4.5算法使用信息增益比率进行特征划分。
KL散度(Kullback–Leibler Divergence, 相对熵):
随机变量x取值的两个概率分布p和q,用来衡量这2个分布的差异,通常用p表示真实分布,用q表示预测分布。
公式定义:
n为事件的所有可能性,如果两个分布完全相同,那么它们的相关熵为0。如果相对熵KL越大,说明它们之间的差异越大,反之相对熵KL越小,说明它们之间的差异越小。
交叉熵(Cross Entropy):
是对KL散度的一种变型,把KL散度log(p(x)/q(x))进行拆分,前面部分就是p的熵H(p),后面就是交叉熵H(p,q)。
公式定义:
交叉熵可以用来计算学习模型分布与训练分布之间的差异,一般在机器学习中直接用交叉熵做损失函数,用于评估模型。
信息论是通信理论的基础,也是xx的基础,关于信息论的理论,等后面有时时间再做分享,本文重要研究信息熵的函数计算问题。
3. 信息度量函数
philentropy包的函数,主要分为3种类别的函数,第一类是距离测量的函数,第二类是相关性分析,第三类是信息度量函数,本文重点介绍这些信息度量的函数。有关于距离测量函数和相关性分析函数,请参考文章
R语言实现46种距离算法
我们来看一下,philentropy包里信息度量的函数:
-
H(): 香农熵, Shannon’s Entropy H(X)
-
JE() : 联合熵, Joint-Entropy H(X,Y)
-
CE() : 条件熵, Conditional-Entropy H(X|Y)
-
MI() : 互信息, Shannon’s Mutual Information I(X,Y)
-
KL() : KL散度, Kullback–Leibler Divergence
-
JSD() : JS散度,Jensen-Shannon Divergence
-
gJSD() : 通用JS散度,Generalized Jensen-Shannon Divergence
本文的系统环境为:
-
Win10 64bit
-
R: 3.4.2 x86_64-w64-mingw32
3.1 H()香农熵
H()函数,可用于快速计算任何给定概率向量的香农熵。
H()函数定义:
H (x, unit = "log2")
参数列表:
- x, 概率向量
- unit,对数化的单位,默认为log2
函数使用:
# 创建数据x
> x<-1:10;x
[1] 1 2 3 4 5 6 7 8 9 10
> px<-x/sum(x);x1
[1] 0.01818182 0.03636364 0.05454545 0.07272727
[5] 0.09090909 0.10909091 0.12727273 0.14545455
[9] 0.16363636 0.18181818
> H(px)
[1] 3.103643
同样地,我们也可以用程序实现公式自己算一下。
# 创建数据x
> x<-1:10
#计算x的概率密度px
> px<-x/sum(x)
> -1sum(pxlog2(px))
[1] 3.103643
我们动手的计算结果,用于H()函数的计算结果是一致的。
3.2 CE()条件熵
CE()函数,基于给定的联合概率向量P(X,Y)和概率向量P(Y),根据公式 H(X|Y)= H(X,Y)-H(Y)计算香农的条件熵。
函数定义:
CE(xy, y, unit = "log2")
参数列表:
- xy, 联合概率向量
- y, 概率向量,必须是随机变量y的概率分布
- unit,对数化的单位,默认为log2
函数使用:
> x3<- 1:10/sum(1:10)
> y3<- 30:40/sum(30:40)
> CE(x3, y3)
[1] -0.3498852
3.3 JE()联合熵
JE()函数,基于给定的联合概率向量P(X,Y)计算香农的联合熵H(X,Y)。
JE()函数定义:
JE (x, unit = "log2")
参数列表:
- x, 联合概率向量
- unit,对数化的单位,默认为log2
函数使用:
# 创建数据x
> x2 <- 1:100/sum(1:100)
> JE(x2)
[1] 6.372236
3.4 MI()互信息
MI()函数,根据给定联合概率向量P(X,Y)、概率向量P(X)和概率向量P(X),按公式I(X,Y)= H(X)+ H(Y)-H(X,Y)计算。
函数定义:
MI(x, y, xy, unit = "log2")
参数列表:
- x, 概率向量
- x, 概率向量
- xy, 联合概率向量
- unit,对数化的单位,默认为log2
函数使用:
# 创建数据集
> x3 <- 1:10/sum(1:10)
> y3<- 20:29/sum(20:29)
> xy3 <- 1:10/sum(1:10)
> MI(x3, y3, xy3)
[1] 3.311973
3.5 KL()散度
KL()函数,计算两个概率分布P和Q的Kullback-Leibler散度。
函数定义:
KL(x, test.na = TRUE, unit = "log2", est.prob = NULL)
参数列表:
- x, 概率向量或数据框
- test.na, 是否检查NA值
- unit,对数化的单位,默认为log2
- est.prob, 用计数向量估计概率的方法,默认值NULL。
函数使用:
# 创建数据集
> df4 <- rbind(x3,y3);df4
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
x3 0.01818182 0.03636364 0.05454545 0.07272727 0.09090909 0.1090909 0.1272727 0.1454545 0.1636364
y3 0.08163265 0.08571429 0.08979592 0.09387755 0.09795918 0.1020408 0.1061224 0.1102041 0.1142857
[,10]
x3 0.1818182
y3 0.1183673
> KL(df4, unit = “log2”) # Default
kullback-leibler
0.1392629
> KL(df4, unit = “log10”)
kullback-leibler
0.0419223
> KL(df4, unit = “log”)
kullback-leibler
0.09652967
3.5 JSD()散度
JSD()函数,基于具有相等权重的Jensen-Shannon散度,计算距离矩阵或距离值。
公式定义:
函数定义:
JSD(x, test.na = TRUE, unit = "log2", est.prob = NULL)
参数列表:
- x, 概率向量或数据框
- test.na, 是否检查NA值
- unit, 对数化的单位,默认为log2
- est.prob, 用计数向量估计概率的方法,默认值NULL。
# 创建数据
> x5 <- 1:10
> y5 <- 20:29
> df5 <- rbind(x5,y5)
> JSD(df5,unit=‘log2’)
jensen-shannon
50.11323
> JSD(df5,unit=‘log’)
jensen-shannon
34.73585
> JSD(df5,unit=‘log10’)
jensen-shannon
15.08559
> JSD(df5, est.prob = “empirical”)
jensen-shannon
0.03792749
3.6 gJSD()散度
gJSD()函数,计算概率矩阵的广义Jensen-Shannon散度。
公式定义:
函数定义:
gJSD(x, unit = "log2", weights = NULL)
参数列表:
- x, 概率矩阵
- unit, 对数化的单位,默认为log2
- weights, 指定x中每个值的权重,默认值NULL。
# 创建数据
> Prob <- rbind(1:10/sum(1:10), 20:29/sum(20:29), 30:39/sum(30:39))
> gJSD(Prob)
[1] 0.023325
4. 应用举例
在我们了解了熵的公式原理和使用方法后,我们就可以做一个案例来试一下。我们定义一个场景的目标:通过用户的看书行为,预测用户是否爱玩游戏。通过我们一步一步地推倒,我们计算出熵,条件熵,联合熵,互信息等指标。
第一步,创建数据集为2列,X列用户看书的类型,包括旅游(Tourism)、美食(Food)、IT技术(IT),Y列用户是否喜欢打游戏,喜欢(Y),不喜欢(N)。
Tourism,Y
Food,N
Tourism,N
Tourism,N
Food,N
Tourism,Y
第二步,建立联合概率矩阵,分别计算H(X),Y(X)。
X | Y | N | p(X) |
Tourism | 2/8=0.25 | 2/8=0.25 | 0.25+0.25=0.5 |
Food | 0/8=0 | 2/8=0.25 | 0+0.25=0.25 |
IT | 2/8=0.25 | 0/8=0 | 0.25+0=0.25 |
p(Y) | 0.25+0+0.25=0.5 | 0.25+0.25+0=0.5 | |
# 分别计算每种情况的概率
p(X=Tourism) = 2/8 + 2/8 = 0.5
p(X=Food) = 2/8 + 0/8 = 0.25
p(X=IT) = 0/8 + 2/8 = 0.25
p(Y=Y) = 4/8 = 0.5
p(Y=N) = 4/8 = 0.5
H(X) = -∑p(xi)*log2(p(xi))
= -p(X=Tourism)log2(p(X=Tourism) ) -p(X=Food)log2(p(X=Food) ) -p(X=IT)log2(p(X=IT) )
= -0.5log(0.5) -0.25log(0.25) - 0.25log(0.25)
= 1.5
H(Y) = -∑p(yi)*log2(p(yi))
= -p(Y=Y)log2(p(Y=Y)) -p(Y=N)log2(p(Y=N))
= -0.5log(0.5) -0.5log(0.5)
= 1
第三步,计算每一项的条件熵,H(Y|X=Tourism),H(Y|X=Food),H(Y|X=IT)。
H(Y|X=Tourism) = -p(Y|X=Tourism)*log(p(Y|X=Tourism)) - p(N|X=Tourism)*log(p(N|X=Tourism))
= -0.5*log(0.5) -0.5*log(0.5)
H(Y|X=Food) = -p(Y|X=Food)log(p(Y|X=Food)) -p(N|X=Food)log(p(N|X=Food))
= -0log(0) -1log(1)
= 0
H(Y|X=IT) = -p(Y|X=IT)log(p(Y|X=IT)) -p(N|X=IT)log(p(N|X=IT))
= -1log(1) -0log(0)
= 0
第四步,计算条件熵H(Y|X)
H(Y|X) = ∑p(xi)*H(Y|xi)
= p(X=Tourism)*H(Y|X=Tourism) + p(X=Food)*H(Y|X=Food) + p(X=IT)*H(Y|X=IT)
= 0.5*1 + 0.25*0 + 0.25*0
= 0.5
第五步,计算联合熵H(X,Y)
H(X,Y) = −∑p(x,y)log(p(x,y))
= H(X) + H(Y|X)
= 1.5 + 0.5
第六步,计算互信息I(X;Y)
I(X;Y) = H(Y) - H(Y|X) = 1 - 0.5 = 0.5
= H(X) + H(Y) - H(X,Y) = 1.5 + 1 - 2 = 0.5
我们把上面的推到过程,用程序来实现一下。
# 创建数据集
> X<-c('Tourism','Food','IT','Tourism','Tourism','IT','Food','Tourism')
> Y<-c('Y','N','Y','N','N','Y','N','Y')
> df<-cbind(X,Y);df
X Y
[1,] "Tourism" "Y"
[2,] "Food" "N"
[3,] "IT" "Y"
[4,] "Tourism" "N"
[5,] "Tourism" "N"
[6,] "IT" "Y"
[7,] "Food" "N"
[8,] "Tourism" "Y
变型为频率矩阵
> tf<-table(df[,1],df[,2]);tf
Food 2 0
IT 0 2
Tourism 2 2
计算概率矩阵
> pX<-margin.table(tf,1)/margin.table(tf);pX
Tourism Food IT
0.50 0.25 0.25
> pY<-margin.table(tf,2)/margin.table(tf);pY
Y N
0.5 0.5
> pXY<-prop.table(tf);pXY
Y N
Tourism 0.25 0.25
Food 0.00 0.25
IT 0.25 0.00
> H(pX)
[1] 1.5
> H(pY)
[1] 1
> CE(pX,pY)
[1] 0.5
> JE(pXY)
[1] 2
> MI(pX,pY,pXY)
[1] 0.5
计算原理是复杂的,用R语言的程序实现却是很简单的,几行代码就搞定了,
本文只是对的信息论的初探,重点还是在信息度量方法的R语言实现。信息熵作为信息度量的基本方法,对各种主流的机器学习的算法都有支撑,是我们必须要掌握的知识。了解本质才能发挥数据科学的潜力,学习的路上不断积累和前进。
转载请注明出处:
原文链接:http://blog.fens.me/r-entropy/
matlab提取文件要素代码传递熵部分信息分解
针对单个试验的时间序列输入矩阵或包含与多个试验对应的多个矩阵的输入像元,计算传递熵的部分信息分解(PID)。
冗余部分信息项由Timme等人(2016年)描述的最小信息函数给出。
MATLAB
R2019a:此处找到的所有函数均为.m文件。
调用了各种MATLAB内置函数。
我们用时间序列/峰值训练来识别神经元。
要计算所有可能的神经元三元组的传递熵PID,请使用3个必填参数调用TE_PID.m
:输出文件名,矩阵或单元格以及正整数时延。
对于包含用于多个试验的多个矩阵的输入像元,输入像元必须为一维。
每个矩阵或单元格列应包含单个神经元的整个时间序列,即,列应代表神经元,而行则代表递增时间的观察值。
(可选)提供要为其计算PID的神经元三重态索引的列表。
否则,将为所有可能的三元组计算PID。
(可选)提供正整数时间分辨率,以便对输入数据进行分时。
输出被写入一个单独的文件。
7列按升序表示:
target_index
source1_index
source2_index
synergy
作者:张丹,R语言中文社区专栏特邀作者,《R的极客理想》系列图书作者,民生银行大数据中心数据分析师,前况客创始人兼CTO。个人博客 http://fens.me, Ale...
转移熵@TOC
转移熵 Transfer entropy(也可译为传递熵)是衡量两个随机过程之间有向(时间不对称)信息传递量的非参数统计量。13过程X到过程Y的转移熵是指在给定过去值Y得到过去值X时,Y值不确定性的减少量。更具体地,如果Xt和Yt(t∈N)表示两个随机过程,且信息量用香农熵 Shannon entropy测量,则转移熵可以写为:
其中 H (x)是 x 的香农熵。上述转移熵的定义已被其他类型的熵测度(如 Rényi熵 Rényi entropy)所扩展。3
转移熵是条件[5][6]互信息
香农编码是是采用信源符号的累计概率分布函数来分配字码的。香农编码是根据香农第一定理直接得出的,指出了平均码长与信息之间的关系,同时也指出了可以通过编码使平均码长达到极限值。香农第一定理是将原始信源符号转化为新的码符号,使码符号尽量服从等概分布,从而每个码符号所携带的信息量达到最大,进而可以用尽量少的码符号传输信源信息。
香农编码属于不等长编码,通常将经常出现的
目录α多样性简介R语言的安装R依赖包及需要命令数据导入数据导入注意事项数据塑形Alpha多样性指数的计算-计算和储存数据可视化可视化-数据导入可视化-数据塑形(合并)可视化-箱线图可视化-添加字母显著性标记可视化-图片的修饰可视化-图片保存总结
α多样性简介
α多样性描述的是样本内的多样性,主要有两个维度
种类数量(丰富度)
种类数量的均匀度
种类数量越多,种类间的分布越均衡,α多样性指数越高
Observed species指数:指样本中实际测定得到的OTU数量,亨利OTU丰富度指数
chaol指数:预测
16S扩增子测序相对已经比较普遍啦,16S数据分析也慢慢的变得常态化,更多的时候我们会根据测序公司给出的原始数据进行分析,或者采用murther等软件,在线的分析工具microbiomeanalyst的功能也比较强大 (https://www.microbiomeanalyst.ca/MicrobiomeAnalyst/home.xhtml)。能够满足基本的一些分析需要,本文以vegan包为基础,为大家提供了分析微生物数据α多样性的一种方法。以下是具体的内容:
#安装R包
install.p
复杂的世界
我们生活在一个极其复杂的世界,不管是小到分子原子亦或是大到整个宇宙,其复杂程度都是超乎想象。或许你不曾深入去思考过身边事物的复杂性,那是因为你已经对你日常所见习以为常。所有在你出生之前发明的事物都是这个世界的自然组成部分,所以很多事物给人的感觉都是这个世界本来的样子。
像人类这样复杂的事物时如何出现的?像电脑这般复杂的事物是如何出现的?像大河山川那样复杂的事物又是怎样出现的呢?某些事...
香农编码概念:香农编码是是采用信源符号的累计概率分布函数来分配字码的。香农编码是根据香农第一定理直接得出的,指出了平均码长与信息之间的关系,同时也指出了可以通过编码使平均码长达到极限值。香农第一定理是将原始信源符号转化为新的码符号,使码符号尽量服从等概分布,从而每个码符号所携带的信息量达到最大,进而可以用尽量少的码符号传输信源信息。香农编码属于不等长编码,通常将经常出现的消息变成短码,不经常出现的...
香农编码(Shannon-Fano coding)是一种编码方式,用于将信源符号(例如字符或单词)转换为二进制位序列。香农编码是基于每个符号的出现频率来构建编码表的。符号出现频率越高,对应的编码就越短。
费诺编码(Huffman coding)是另一种用于将信源符号转换为二进制位序列的编码方式。与香农编码类似,费诺编码也是基于每个符号的出现频率来构建编码表的。不同的是,费诺编码使用了一种名为最小堆...