• 粉丝日志龙虎大战坐庄首页

用R语言解读统计检验-T检验

R的极客理想系列文章,涵盖了R的思想,使用,龙虎大战坐庄工具 ,创新等的一系列要点,以龙虎大战坐庄我 个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,龙虎大战坐庄互联网 ….都在使用R语言。

要成为有理想的极客,龙虎大战坐庄龙虎大战坐庄我 们 不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让龙虎大战坐庄龙虎大战坐庄我 们 一起动起来吧,开始R的极客理想。

龙虎大战坐庄关于 作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://cnc77.com
  • email: bsspirit@gmail.com

转载请注明出处:
http://cnc77.com/r-test-t

前言

做统计分析R语言是最好的,R语言对于统计检验有非常好的支持。龙虎大战坐庄我 会分7篇文章,分别介绍用R语言进行统计量和统计检验的计算过程,包括T检验F检验卡方检验P值KS检验AICBIC等常用的统计检验龙虎大战坐庄方法 和统计量。

本文先从T检验开始说起!

目录

  1. T检验介绍
  2. 数据集
  3. 单总体T检验
  4. 双总体T检验

1. T检验介绍

T检验,也称 student t检验(Student’s t test),是戈斯特为了观测酿酒质量发明的。戈斯特于1908年在Biometrika上公布T检验,但因其老板认为其为商业机密而被迫使用笔名(龙虎大战坐庄学生 )。

T检验,是用于检验两个小样本的平均值差异程度的检验龙虎大战坐庄方法 ,样本量在3-30个左右,要求样本为正态分布,但总体标准差可未知。T检验,是用T分布理论来推断差异发生的概率,从而判定两个样本平均值的差异是否显著。T检验可分为单总体T检验,双总体非配对T检验,和双总体配对T检验。下面将分别进行介绍。

2. 数据集

T检验,对于数据有比较严格的要求,所以龙虎大战坐庄龙虎大战坐庄我 们 需要先找到一个合适的数据集,作为测试数据集。龙虎大战坐庄我 发现了R语言自带的一个数据集sleep,是很好的测试数据集,本文接下来的内容,将以这个数据集进行测试,来介绍T检验。

开发环境所使用的系统环境

  • 龙虎大战坐庄Win 10 64bit
  • R: 3.4.2 x86_64-w64-mingw32/x64 b4bit

数据集,记录了10名患者使用两种催眠药物影响的对比数据,共3列,20条记录。

  • extra列,为睡眠时间的变化
  • group列,为药物编号
  • ID列,为患者编号

查看数据集。


> data(sleep)
> sleep
   extra group ID
1    0.7     1  1
2   -1.6     1  2
3   -0.2     1  3
4   -1.2     1  4
5   -0.1     1  5
6    3.4     1  6
7    3.7     1  7
8    0.8     1  8
9    0.0     1  9
10   2.0     1 10
11   1.9     2  1
12   0.8     2  2
13   1.1     2  3
14   0.1     2  4
15  -0.1     2  5
16   4.4     2  6
17   5.5     2  7
18   1.6     2  8
19   4.6     2  9
20   3.4     2 10

由于T检验的前提条件,总体要符合正态分布,只是总体方差未知。所以,龙虎大战坐庄龙虎大战坐庄我 们 需要先对选定数据集进行进行正态分布检验。使用Shapiro-Wilk和qq图,作为正态分布检验的龙虎大战坐庄方法 。

Shapiro-Wilk正态分布检验: 用来检验是否数据符合正态分布,类似于线性回归的龙虎大战坐庄方法 一样,是检验其于回归曲线的残差。该龙虎大战坐庄方法 龙虎大战坐庄推荐 在样本量很小的时候使用,样本在3到5000之间。该检验原假设为H0:数据集符合正态分布,统计量为W。统计量W最大值是1,越接近1,表示样本与正态分布匹配。p值,当p-value小于显著性水平α(0.05),则拒绝H0。


> shapiro.test(sleep$extra) # 总体正态分布

	Shapiro-Wilk normality test

data:  sleep$extra
W = 0.94607, p-value = 0.3114

检验结果,W接近1,p-value>0.05,不能拒绝原假设,所以数据集sleep是符合正态分布的。

接下来,龙虎大战坐庄龙虎大战坐庄我 们 再画出qq图,直观的看一下数据符合正态分布的情况。

图中,对角线基本能穿过数据点,也说明数据符合正态分布。

3. 单总体的T检验

目的:单总体T检验,是检验一个样本平均值与一个已知的总体平均值的差异是否显著。当总体分布是正态分布,如总体标准差未知且样本容量较小,那么样本平均数与总体平均数的离差统计量呈t分布。

适用条件:

  1. 已知总体均值Um
  2. 可得到样本均值Xm,及该样本标准误差Xs
  3. 样本来自正态或近似正态总体

计算公式:T统计量

公式解释:

  • Xm:样本均值
  • Um:总体均值
  • Xs:样本标准差
  • n:样本个数
  • 自由度:df=n-1

用R语言进行T检验,取sleep数据集做为总体,符合正态分布,总体平均值为1.54,以group=1为样本组,共10条记录,样本均值为0.75,样本标准差为1.78901。判断group=1分组的药物,是否有显著性改进睡眠?


#  总体均值
> Um<-mean(sleep$extra) 

# 取group=1的样本组
> g1<-sleep$extra[which(sleep$group==1)]

# 计算样本均值,方差,数量
> Xm<-mean(g1)
> Xs<-sd(g1)
> n<-length(g1)

# 进行T检验
> t.test(g1,mu=Um)

	One Sample t-test

data:  g1
t = -1.3964, df = 9, p-value = 0.1961
alternative hypothesis: true mean is not equal to 1.54
95 percent confidence interval:
 -0.5297804  2.0297804
sample estimates:
mean of x 
     0.75 

指标解释:

  • H0:原假设总体和样本均值无显著差异
  • t统计量:-1.3964
  • df,自由度,10-1=9
  • p-value值:0.1961
  • 95 percent confidence interval:95%的置信区间
  • mean:样本均值0.75

结果解读,以0.05为显著性水平,t=-1.3964小于临界值2.228(查表),t值不显著,不能拒绝原假设。以0.05为显著性水平,p-value=0.19大于0.05,不能拒绝原假设,group=1的药物对睡眠没有明显改进。

龙虎大战坐庄龙虎大战坐庄我 们 可以用t值进行显著性差异判断,也可以用p值进行显著性差异判断,他们的作用是一样的。t值判断时,需要用计算所得的t值,与显著性水平查表对比。p值相当于是把t值,进行一种标准化的变型,只和已经定义好的显著性水平比就行了,比如0.05, 0.01, 0.001等几个固定值。

手动计算T值和P值,龙虎大战坐庄关于 P值的详细解释,请查看文章R语言实现统计检验-P值


# 手动计算T值
> t_stat<-(Xm-Um) / (Xs/(sqrt(n)));t_stat 
[1] -1.396415

# 手动计算P值
> p_value<-2*pt(-abs(t_stat),df=n-1);p_value    
[1] 0.1960699

4. 双总体T检验

双总体T检验,是检验两个样本平均值,与其各自所代表的总体的差异是否显著。双总体T检验又分为两种情况,一种是配对的样本T检验,用于检验两种同质对象,在不同条件下所产生的数据差异性;另一种是独立样本非配对T检验,用于检验两组独立的样本的平均数差异性。

4.1 配对T检验

目标:检验两组同质样本,在不同的处理下的样本平均值,是否有显著的差异性。

配对设计:将2组样本的某些重要特征按相近的原则配成对子,消除混杂因素的影响,观察样本之间的处理因素和研究因素的差异,其它因素基本相同,把配对两组样本个体随机处理。

配对过程如下3种情况:

  • 两种同质样本,分别接受两种不同的处理,如性别、年龄、体重、病情程度相同进行配对。
  • 同一样本或同一样本的两个部分,分别接受两种不同的处理。
  • 同一样本自身对比,把同一组样本处理前后的结果进行比较。

计算公式:

公式解释:

  • D:两个样本差值
  • εD:求和
  • ε(D^2):平方和
  • n:样本个数
  • 自由度:df=n-1

使用sleep数据集,按group分成2个组,形成配对的数据集。H0原假设,g1和g2两个样本组均值没有显著性差异,对治疗睡眠的效果是一致的。


# 配对T检验
> t.test(extra ~ group, data = sleep, paired = TRUE)

	Paired t-test

data:  extra by group
t = -4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.4598858 -0.7001142
sample estimates:
mean of the differences 
                  -1.58 

结果解读,以0.05为显著性水平,p-value=0.002833<0.05,说明g1和g2有显著性的差异,可以拒绝原假设。

接下来,龙虎大战坐庄龙虎大战坐庄我 们 动手自己算一下T值


# 生成数据集
> g1<-sleep$extra[which(sleep$group==1)]
> g2<-sleep$extra[which(sleep$group==2)]

# 配对T检验,与上面的计算结果是一致的
> t.test(g1,g2,paired=TRUE)

	Paired t-test

data:  g1 and g2
t = -4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.4598858 -0.7001142
sample estimates:
mean of the differences 
                  -1.58 

# 动手计算
> n1<-length(g1)
> n2<-length(g2)
> D<-g1-g2
> Ds<-sum(D)

# 计算T值
> Ds/sqrt((n1*sum(D^2)-sum(D)^2)/(n1-1))
[1] -4.062128

4.2 非配对T检验

目标:用于检验两组独立的样本的平均数差异性。

计算公式:

  • n1 和n2 为两样本容量
  • Xm1和Xm2为两样本均值
  • Xs1和Xs2为两样本的标准差

继续使用sleep数据集,按group分成2个组,形成独立的样本,病人按随机进行配对,去掉关联关系。H0原假设,随机选取g1和g2两个样本组均值没有显著性差异,对治疗睡眠的效果是一致的。


> t.test(extra ~ group, data = sleep)

	Welch Two Sample t-test

data:  extra by group
t = -1.8608, df = 17.776, p-value = 0.07939
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -3.3654832  0.2054832
sample estimates:
mean in group 1 mean in group 2 
           0.75            2.33 

结果解读,以0.05为显著性水平,p-value=0.07939>0.05,说明g1和g2有没显著性的差异,不能拒绝原假设,两种龙虎大战坐庄方法 对于治疗效果是一致的。

接下来,龙虎大战坐庄龙虎大战坐庄我 们 动手自己算一下T值。


# 计算数量,均值,标准差
> n1<-length(g1)
> n2<-length(g2)
> Xm1<-mean(g1)
> Xm2<-mean(g2)
> Xs1<-sd(g1)
> Xs2<-sd(g2)
 
# 计算T值
> (Xm1-Xm2)/sqrt((Xs1^2/n1)+(Xs2^2/n2))
[1] -1.860813

画出箱线图,观察两组数据的均值。


> boxplot(extra~group,data=sleep)


图中,每个箱的最粗的黑色线,表示两个样本的均值。

以双总体进行T检验,当是样本是配对时,g1和g2有显著性差异;而当样本是独立时,g1和g2没有显著性差异。由于龙虎大战坐庄龙虎大战坐庄我 们 变化了初始的假设,是会很大的程度影响统计结果的,所以在使用统计学模型时,要做非常严格的条件判断,从而保证结果的可靠性,可解释性。

转载请注明出处:
http://cnc77.com/r-test-t

打赏作者

R语言中的矩阵计算

R的极客理想系列文章,涵盖了R的思想,使用,龙虎大战坐庄工具 ,创新等的一系列要点,以龙虎大战坐庄我 个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,龙虎大战坐庄互联网 ….都在使用R语言。

要成为有理想的极客,龙虎大战坐庄龙虎大战坐庄我 们 不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让龙虎大战坐庄龙虎大战坐庄我 们 一起动起来吧,开始R的极客理想。

龙虎大战坐庄关于 作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://cnc77.com
  • email: bsspirit@gmail.com

转载请注明出处:
http://cnc77.com/r-matrix

前言

R是作为统计语言,生来就对数学有良好的支持。矩阵计算作为底层的数学龙虎大战坐庄工具 ,有非常广泛的使用场景。用R语言很好地封装了,矩阵的各种计算龙虎大战坐庄方法 ,一个函数一行代码,就能完成复杂的矩阵分解等操作。让建模人员可以更专注于模型推理和业务逻辑实现,把复杂的矩阵计算交给R语言来完成。

本文总结了R语言用于矩阵的各种计算操作。

目录

  1. 基本操作
  2. 矩阵计算
  3. 矩阵性质
  4. 矩阵分解
  5. 特殊矩阵

1. 基本操作

生成一个矩阵,计算行、列。


# 生成矩阵 
> m<-matrix(1:20,4,5);m
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20

# 矩阵行
> nrow(m)
[1] 4
# 矩阵列
> ncol(m)
[1] 5

取对角线元素,生成对角矩阵,


# 对角线元素
> diag(m)
[1]  1  6 11 16

# 以对角线元素,生成对角矩阵
> diag(diag(m))
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    6    0    0
[3,]    0    0   11    0
[4,]    0    0    0   16

上三角,下三角


# 上三角
> m<-matrix(1:20,4,5)
> upper.tri(m)
      [,1]  [,2]  [,3]  [,4] [,5]
[1,] FALSE  TRUE  TRUE  TRUE TRUE
[2,] FALSE FALSE  TRUE  TRUE TRUE
[3,] FALSE FALSE FALSE  TRUE TRUE
[4,] FALSE FALSE FALSE FALSE TRUE

# 上三角值
> m[which(upper.tri(m))] 
 [1]  5  9 10 13 14 15 17 18 19 20

# 上三角矩阵
> m[!upper.tri(m)]<-0;m
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    5    9   13   17
[2,]    0    0   10   14   18
[3,]    0    0    0   15   19
[4,]    0    0    0    0   20

# 下三角矩阵
> m[!lower.tri(m)]<-0;m
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    2    0    0    0    0
[3,]    3    7    0    0    0
[4,]    4    8   12    0    0

矩阵转置


> m<-matrix(1:20,4,5);m
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20

# 转置,行列互转
> t(m)
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    5    6    7    8
[3,]    9   10   11   12
[4,]   13   14   15   16
[5,]   17   18   19   20

对角矩阵填充

# 创建方阵
> m<-matrix(1:16,4,4);m
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

# 用上三角填充下三角
> m[lower.tri(m)]<-m[upper.tri(m)]
> m
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    5    6   10   14
[3,]    9   13   11   15
[4,]   10   14   15   16

填充后,发现矩阵并不是对称的,原因是上三角取值按列取值,所以先取10后取13,导致上三角和下三角取值顺序不完全一致。


> m<-matrix(1:16,4,4);m
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

# 上三角取值
> m[upper.tri(m)]
[1]  5  9 10 13 14 15

# 下三角取值
> m[lower.tri(m)]
[1]  2  3  4  7  8 12

调整后,龙虎大战坐庄龙虎大战坐庄我 们 要先转置,再取值再填充,形成对称结构。


> m<-matrix(1:20,4,5)

# 转置后,取下三角,填充下三角
> m[lower.tri(m)]<-t(m)[lower.tri(m)]
> m
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    5    6   10   14
[3,]    9   10   11   15
[4,]   13   14   15   16

矩阵和data.frame转换,用行列形成索引结构


> m<-matrix(1:12,4,3);m
     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12

# 矩阵转行列data.frame
> row<-rep(1:nrow(m),ncol(m))              # 行
> col<-rep(1:ncol(m),each=nrow(m))         # 列
> df<-data.frame(row,col,v=as.numeric(m))  # 行列索引数据框
> df
   row col  v
1    1   1  1
2    2   1  2
3    3   1  3
4    4   1  4
5    1   2  5
6    2   2  6
7    3   2  7
8    4   2  8
9    1   3  9
10   2   3 10
11   3   3 11
12   4   3 12

# 行列索引数据框转矩阵
> m<-matrix(0,length(unique(df$row)),length(unique(df$col)) )
> apply(df,1,function(dat){
+     m[dat[1],dat[2]]<<-dat[3]  
+     invisible()
+ })
# 打印矩阵
> m
     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12

2. 矩阵计算

加法,减法


# 加载矩阵计算龙虎大战坐庄工具
包
> library(matrixcalc)

# 新建2个矩阵,行列长度相同
> m0<-matrix(1:20,4,5);m0
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20
> m1<-matrix(sample(100,20),4,5);m1
     [,1] [,2] [,3] [,4] [,5]
[1,]   40   79   97   57   78
[2,]   93   32   48    8   95
[3,]   63    6   56   12    9
[4,]   28   31   72   27   26

# 矩阵加法
> m0+m1
     [,1] [,2] [,3] [,4] [,5]
[1,]   41   84  106   70   95
[2,]   95   38   58   22  113
[3,]   66   13   67   27   28
[4,]   32   39   84   43   46

# 矩阵减法
> m0-m1
     [,1] [,2] [,3] [,4] [,5]
[1,]  -39  -74  -88  -44  -61
[2,]  -91  -26  -38    6  -77
[3,]  -60    1  -45    3   10
[4,]  -24  -23  -60  -11   -6

矩阵值相乘


> m0*m1
     [,1] [,2] [,3] [,4] [,5]
[1,]   40  395  873  741 1326
[2,]  186  192  480  112 1710
[3,]  189   42  616  180  171
[4,]  112  248  864  432  520

矩阵乘法,满足第二个矩阵的列数和第一个矩阵的行数相等,所以把上面生成的m0矩阵(4行5列)转置为(5行4列),再用m1矩阵(4行5列),进行矩阵乘法,得到一个5行5列的结果矩阵。


> t(m0)%*%m1
     [,1] [,2] [,3] [,4] [,5]
[1,]  527  285  649  217  399
[2,] 1423  877 1741  633 1231
[3,] 2319 1469 2833 1049 2063
[4,] 3215 2061 3925 1465 2895
[5,] 4111 2653 5017 1881 3727

# 通过函数实现矩阵相乘
> crossprod(m0,m1)
     [,1] [,2] [,3] [,4] [,5]
[1,]  527  285  649  217  399
[2,] 1423  877 1741  633 1231
[3,] 2319 1469 2833 1049 2063
[4,] 3215 2061 3925 1465 2895
[5,] 4111 2653 5017 1881 3727

矩阵外积


> m<-matrix(1:6,2);m
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> t(m)
     [,1] [,2]
[1,]    1    2
[2,]    3    4
[3,]    5    6

> m %o% t(m)  # 外积,同outer(m,t(m))
, , 1, 1

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2, 1

     [,1] [,2] [,3]
[1,]    3    9   15
[2,]    6   12   18

, , 3, 1

     [,1] [,2] [,3]
[1,]    5   15   25
[2,]   10   20   30

, , 1, 2

     [,1] [,2] [,3]
[1,]    2    6   10
[2,]    4    8   12

, , 2, 2

     [,1] [,2] [,3]
[1,]    4   12   20
[2,]    8   16   24

, , 3, 2

     [,1] [,2] [,3]
[1,]    6   18   30
[2,]   12   24   36

矩阵直和


> m0<-matrix(1:4,2,2);m0
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> m1<-matrix(1:9,3,3);m1
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

> m0 %s% m1       #矩阵直和,同 direct.sum(m0,m1) 
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    3    0    0    0
[2,]    2    4    0    0    0
[3,]    0    0    1    4    7
[4,]    0    0    2    5    8
[5,]    0    0    3    6    9

矩阵直积


> m0<-matrix(1:4,2,2);m0
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> m1<-matrix(1:9,3,3);m1
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

> m0 %x% m1     #矩阵直积,同 kronecker(m0,m1),direct.prod(m0,m1)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    4    7    3   12   21
[2,]    2    5    8    6   15   24
[3,]    3    6    9    9   18   27
[4,]    2    8   14    4   16   28
[5,]    4   10   16    8   20   32
[6,]    6   12   18   12   24   36

3. 矩阵性质

3.1 奇异矩阵
首先,龙虎大战坐庄龙虎大战坐庄我 们 线创建一个非奇异矩阵,判断非奇异矩阵龙虎大战坐庄方法 ,行列式不等于0,矩阵可逆,满秩。


# 创建一个非奇异矩阵
> m1<-matrix(sample(1:16),4)

# 行列式不等于0
> det(m1)     # !=0
[1] 14193

# 有逆矩阵
> solve(m1)   # 可逆
             [,1]        [,2]        [,3]         [,4]
[1,] -0.026210104  0.09666737  0.02473050 -0.119988727
[2,] -0.002325090  0.08384415 -0.07038681  0.008173043
[3,] -0.007186641 -0.13478475  0.05516804  0.176777285
[4,]  0.074614246  0.03663778 -0.01395054 -0.080462200

# 满秩
> library(matrixcalc)
> matrix.rank(m1)    # 长度=n,满秩
[1] 4

再创建一个奇异矩阵,判断奇异矩阵龙虎大战坐庄方法 包括,行列式等于0,矩阵不可逆,不是满秩。


# 奇异矩阵
> m0<-matrix(1:16,4)

# 计算行列式
> det(m0)     
[1] 0

# 不可逆
> solve(m0)   
Error in solve.default(m0) : 
  Lapack routine dgesv: system is exactly singular: U[3,3] = 0

# 不是满秩
> matrix.rank(m0)     # 长度<4
[1] 2

3.2 逆矩阵


# 创建方阵,非奇异矩阵
> m0<-matrix(sample(100,16),4,4);m0
     [,1] [,2] [,3] [,4]
[1,]   24   31   80   37
[2,]   84   13   42   71
[3,]   95   62   93   86
[4,]   69   16   94   35

# 计算矩阵的逆矩阵
> solve(m0)
            [,1]          [,2]        [,3]         [,4]
[1,] -0.03083680 -0.0076561475  0.01258023  0.017218514
[2,] -0.01710957 -0.0270246488  0.03152548 -0.004553923
[3,]  0.01384721 -0.0003070371 -0.00886117  0.007757524
[4,]  0.03142440  0.0282722871 -0.01541411 -0.024126340

# 逆矩阵的性质,逆矩阵与原矩阵进行矩阵乘法,得到对角矩阵。
> round(solve(m0) %*% m0)  # 对角矩阵
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    1    0    0
[3,]    0    0    1    0
[4,]    0    0    0    1

广义逆矩阵,将逆矩阵的概率推广到奇异矩阵和长方形矩阵上,就产生了广义逆矩阵。


# 创建奇异矩阵
> m<-matrix(1:16,4,4);m 
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

# 逆矩阵计算失败
> solve(m)
Error in solve.default(m) : 
  Lapack routine dgesv: system is exactly singular: U[3,3] = 0

# 广义逆矩阵
> library(MASS) # 加载MASS包
> ginv(m)
       [,1]    [,2]  [,3]    [,4]
[1,] -0.285 -0.1075  0.07  0.2475
[2,] -0.145 -0.0525  0.04  0.1325
[3,] -0.005  0.0025  0.01  0.0175
[4,]  0.135  0.0575 -0.02 -0.0975

用ginv函数计算非奇异矩阵,和solve()函数比较。


# 非奇异矩阵
> m0<-matrix(sample(100,16),4,4);m0
     [,1] [,2] [,3] [,4]
[1,]    5   61   75   74
[2,]   10   67   11   21
[3,]   29   32   92   17
[4,]   72   23   87   36

# 计算广义矩阵的逆矩阵,与solve()结果相同
> ginv(m0)
              [,1]         [,2]         [,3]         [,4]
[1,] -0.0080590817  0.006534814 -0.011333076  0.018105645
[2,] -0.0041284695  0.017615769  0.005333304 -0.004308072
[3,]  0.0009226026 -0.006671759  0.016312821 -0.005707878
[4,]  0.0165261737 -0.008200731 -0.020163889  0.008112906

# 计算矩阵的逆矩阵
> solve(m0)
              [,1]         [,2]         [,3]         [,4]
[1,] -0.0080590817  0.006534814 -0.011333076  0.018105645
[2,] -0.0041284695  0.017615769  0.005333304 -0.004308072
[3,]  0.0009226026 -0.006671759  0.016312821 -0.005707878
[4,]  0.0165261737 -0.008200731 -0.020163889  0.008112906

逆矩阵的特性


> m<-matrix(1:16,4,4);m 
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

# 原矩阵%*%广义逆矩阵%*%原矩阵=原矩阵
> m %*% ginv(m) %*% m 
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

# 广义逆矩阵%*%原矩阵%*%广义逆矩阵=广义逆矩阵
> ginv(m) %*% m %*% ginv(m) 
       [,1]    [,2]  [,3]    [,4]
[1,] -0.285 -0.1075  0.07  0.2475
[2,] -0.145 -0.0525  0.04  0.1325
[3,] -0.005  0.0025  0.01  0.0175
[4,]  0.135  0.0575 -0.02 -0.0975

# 原矩阵%*%广义逆矩阵=转置后的,原矩阵%*%广义逆矩阵
> m %*% ginv(m)    
     [,1] [,2] [,3] [,4]
[1,]  0.7  0.4  0.1 -0.2
[2,]  0.4  0.3  0.2  0.1
[3,]  0.1  0.2  0.3  0.4
[4,] -0.2  0.1  0.4  0.7
> t(m %*% ginv(m)) 
     [,1] [,2] [,3] [,4]
[1,]  0.7  0.4  0.1 -0.2
[2,]  0.4  0.3  0.2  0.1
[3,]  0.1  0.2  0.3  0.4
[4,] -0.2  0.1  0.4  0.7

# 广义逆矩阵%*%原矩阵=转置后的,广义逆矩阵%*%原矩阵
> ginv(m) %*% m   
     [,1] [,2] [,3] [,4]
[1,]  0.7  0.4  0.1 -0.2
[2,]  0.4  0.3  0.2  0.1
[3,]  0.1  0.2  0.3  0.4
[4,] -0.2  0.1  0.4  0.7
> t(ginv(m) %*% m)
     [,1] [,2] [,3] [,4]
[1,]  0.7  0.4  0.1 -0.2
[2,]  0.4  0.3  0.2  0.1
[3,]  0.1  0.2  0.3  0.4
[4,] -0.2  0.1  0.4  0.7

3.3 特征值和特征向量


# 创建一个方阵
> m0<-matrix(sample(100,16),4,4);m0
     [,1] [,2] [,3] [,4]
[1,]   97   93   11   70
[2,]   19   58   23   90
[3,]   17   94    6   35
[4,]   79   71   43   88

# 计算特征值和特征向量
> eigen(m0)
eigen() decomposition
$values
[1] 236.14449+ 0.00000i  40.51501+ 0.00000i -13.82975+32.81166i -13.82975-32.81166i

$vectors
             [,1]          [,2]                  [,3]                  [,4]
[1,] 0.6055193+0i -0.6428709+0i  0.2114869-0.0613776i  0.2114869+0.0613776i
[2,] 0.4115920+0i  0.3939259+0i -0.0781518+0.3993580i -0.0781518-0.3993580i
[3,] 0.3054253+0i  0.6482306+0i  0.7557355+0.0000000i  0.7557355+0.0000000i
[4,] 0.6088134+0i -0.1064728+0i -0.3210016-0.3342655i -0.3210016+0.3342655i

当symmetric=TRUE时,计算对称矩阵的特征值和特征向量,当m0不是对称矩阵时,则取下三角对称结构进行计算。


> eigen(m0,symmetric = TRUE)
eigen() decomposition
$values
[1] 231.824838  87.176816  -3.422899 -66.578756

$vectors
           [,1]       [,2]       [,3]        [,4]
[1,] -0.4801643  0.7155875  0.5040624 -0.05742733
[2,] -0.5024315 -0.5470752  0.2262475 -0.63014551
[3,] -0.3634220 -0.4138943  0.3287896  0.76714627
[4,] -0.6204267  0.1316618 -0.7659181  0.10538188

# 生成下三角矩阵的对称矩阵
> m0[upper.tri(m0)]<-t(m0)[upper.tri(m0)];m0
     [,1] [,2] [,3] [,4]
[1,]   97   19   17   79
[2,]   19   58   94   71
[3,]   17   94    6   43
[4,]   79   71   43   88

# 计算特征值,与eigen(m0,symmetric = TRUE) 一致
> eigen(m0)
eigen() decomposition
$values
[1] 231.824838  87.176816  -3.422899 -66.578756

$vectors
           [,1]       [,2]       [,3]        [,4]
[1,] -0.4801643  0.7155875  0.5040624 -0.05742733
[2,] -0.5024315 -0.5470752  0.2262475 -0.63014551
[3,] -0.3634220 -0.4138943  0.3287896  0.76714627
[4,] -0.6204267  0.1316618 -0.7659181  0.10538188

4. 矩阵分解

下面将介绍4种矩阵常用的分解的龙虎大战坐庄方法 ,包括三解分解LU,choleskey分解,QR分解,奇异值分解SVD。

4.1 三解分解LU
三角分解法是将原方阵分解成一个上三角形矩阵和一个下三角形矩阵,这样的分解法又称为LU分解法。它的用途主要在简化一个大矩阵的行列式值的计算过程,求逆矩阵,和求解联立方程组。这种分解法所得到的上下三角形矩阵不唯一,一对上三角形矩阵和下三角形矩阵,矩阵相乘会得到原矩阵。


创建一个矩阵
> m0<-matrix(sample(100,16),4);m0
     [,1] [,2] [,3] [,4]
[1,]   25   88   35   87
[2,]   26   75   73   15
[3,]   36   62   80   53
[4,]   70   44   21   47

# 三角分解
> lu<-lu.decomposition(m0)
> lu
$L                                        # 下三角矩阵
     [,1]      [,2]     [,3] [,4]
[1,] 1.00  0.000000 0.000000    0
[2,] 1.04  1.000000 0.000000    0
[3,] 1.44  3.917676 1.000000    0
[4,] 2.80 12.251816 4.617547    1

$U                                        # 上三角矩阵
     [,1]   [,2]      [,3]      [,4]
[1,]   25  88.00   35.0000   87.0000
[2,]    0 -16.52   36.6000  -75.4800
[3,]    0   0.00 -113.7869  223.4262
[4,]    0   0.00    0.0000 -303.5137

# 三角分解的下三角矩阵L %*% 三角分解的上三角矩阵U = 原矩阵
> lu$L %*% lu$U
     [,1] [,2] [,3] [,4]
[1,]   25   88   35   87
[2,]   26   75   73   15
[3,]   36   62   80   53
[4,]   70   44   21   47

4.2 choleskey分解
Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解。它要求矩阵的所有特征值必须大于零,故分解的下三角的对角元也是大于零的。Cholesky分解法又称平方根法,是当A为实对称正定矩阵时,LU三角分解法的变形。


创建对称方阵
> m1<-diag(5)+2;m1
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# choleskey分解
> chol(m1)
         [,1]     [,2]      [,3]      [,4]      [,5]
[1,] 1.732051 1.154701 1.1547005 1.1547005 1.1547005
[2,] 0.000000 1.290994 0.5163978 0.5163978 0.5163978
[3,] 0.000000 0.000000 1.1832160 0.3380617 0.3380617
[4,] 0.000000 0.000000 0.0000000 1.1338934 0.2519763
[5,] 0.000000 0.000000 0.0000000 0.0000000 1.1055416

# choleskey分解的性质
# chol分解的逆矩阵 %*% chol分解矩阵 = 原矩阵
> t(chol(m1))%*%chol(m1)
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# chol分解矩阵的对角线值的平方的积 = 行列式的模数
> prod(diag(chol(m1))^2)
[1] 11
> det(m1)
[1] 11

# chol分解的逆
> chol2inv(m1)
             [,1]         [,2]        [,3]         [,4]         [,5]
[1,]  0.166658199 -0.055580958 -0.01859473 -0.006401463 -0.002743484
[2,] -0.055580958  0.166590459 -0.05578418 -0.019204390 -0.008230453
[3,] -0.018594726 -0.055784179  0.16598080 -0.057613169 -0.024691358
[4,] -0.006401463 -0.019204390 -0.05761317  0.160493827 -0.074074074
[5,] -0.002743484 -0.008230453 -0.02469136 -0.074074074  0.111111111

> chol2inv(chol(m1))
           [,1]       [,2]       [,3]       [,4]       [,5]
[1,]  0.8181818 -0.1818182 -0.1818182 -0.1818182 -0.1818182
[2,] -0.1818182  0.8181818 -0.1818182 -0.1818182 -0.1818182
[3,] -0.1818182 -0.1818182  0.8181818 -0.1818182 -0.1818182
[4,] -0.1818182 -0.1818182 -0.1818182  0.8181818 -0.1818182
[5,] -0.1818182 -0.1818182 -0.1818182 -0.1818182  0.8181818

4.3 QR分解
QR分解法是将矩阵分解成一个正规正交矩阵与上三角形矩阵,所以称为QR分解法,与此正规正交矩阵的通用符号Q有关。


# 创建对称方阵
> m1<-diag(5)+2;m1
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# QR分解
> q<-qr(m1);q
$qr
     [,1]       [,2]       [,3]       [,4]       [,5]
[1,] -5.0 -4.8000000 -4.8000000 -4.8000000 -4.8000000
[2,]  0.4 -1.4000000 -0.6857143 -0.6857143 -0.6857143
[3,]  0.4  0.2142857 -1.2205720 -0.4012839 -0.4012839
[4,]  0.4  0.2142857  0.1560549 -1.1527216 -0.2852095
[5,]  0.4  0.2142857  0.1560549  0.1246843  1.1168808

$rank
[1] 5

$qraux
[1] 1.600000 1.928571 1.975343 1.992196 1.116881

$pivot
[1] 1 2 3 4 5

attr(,"class")
[1] "qr"

# 计算QR分解矩阵的R值
> qr.R(q)
     [,1] [,2]       [,3]       [,4]       [,5]
[1,]   -5 -4.8 -4.8000000 -4.8000000 -4.8000000
[2,]    0 -1.4 -0.6857143 -0.6857143 -0.6857143
[3,]    0  0.0 -1.2205720 -0.4012839 -0.4012839
[4,]    0  0.0  0.0000000 -1.1527216 -0.2852095
[5,]    0  0.0  0.0000000  0.0000000  1.1168808

# 计算QR分解矩阵的Q值
> qr.Q(q)
     [,1]        [,2]        [,3]        [,4]       [,5]
[1,] -0.6  0.62857143  0.36784361  0.26144202 -0.2030692
[2,] -0.4 -0.77142857  0.36784361  0.26144202 -0.2030692
[3,] -0.4 -0.05714286 -0.85272836  0.26144202 -0.2030692
[4,] -0.4 -0.05714286 -0.03344033 -0.89127960 -0.2030692
[5,] -0.4 -0.05714286 -0.03344033 -0.02376746  0.9138115

# QR分解的性质
# Q的值 %*% R值 = 原矩阵
> qr.Q(q) %*% qr.R(q) #=m1
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# X值 = 原矩阵,同Q的值 %*% R值 
> qr.X(q) #=m1
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# Q值的转置矩阵 * Q值 = 
> t(qr.Q(q)) %*% qr.Q(q)
              [,1]          [,2]          [,3]         [,4]          [,5]
[1,]  1.000000e+00 -3.469447e-17 -2.081668e-17 1.214306e-17  5.551115e-17
[2,] -3.469447e-17  1.000000e+00  7.199102e-17 5.204170e-17 -7.632783e-17
[3,] -2.081668e-17  7.199102e-17  1.000000e+00 3.382711e-17 -6.938894e-18
[4,]  1.214306e-17  5.204170e-17  3.382711e-17 1.000000e+00  3.122502e-17
[5,]  5.551115e-17 -7.632783e-17 -6.938894e-18 3.122502e-17  1.000000e+00

4.4 奇异值分解SVD
奇异值分解 (singular value decomposition, SVD) 是一种正交矩阵分解法。SVD是最可靠的分解法,但是它比QR 分解法要花上近十倍的计算时间。[U,S,V]=svd(A),其中U和V分别代表两个正交矩阵,而S代表一对角矩阵。 和QR分解法相同, 原矩阵A不必为正方矩阵。使用SVD分解法的用途是解最小平方误差法和数据压缩。


# 创建对称方阵
> m1<-diag(5)+2;m1
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# 奇异值分解
> s<-svd(m1);s
$d
[1] 11  1  1  1  1

$u
           [,1]          [,2]          [,3]          [,4]       [,5]
[1,] -0.4472136  4.440892e-17 -3.330669e-17 -3.108624e-16  0.8944272
[2,] -0.4472136 -3.219647e-16  2.414735e-16  8.660254e-01 -0.2236068
[3,] -0.4472136 -5.773503e-01  5.773503e-01 -2.886751e-01 -0.2236068
[4,] -0.4472136 -2.113249e-01 -7.886751e-01 -2.886751e-01 -0.2236068
[5,] -0.4472136  7.886751e-01  2.113249e-01 -2.886751e-01 -0.2236068

$v
           [,1]          [,2]          [,3]       [,4]       [,5]
[1,] -0.4472136  0.000000e+00  0.000000e+00  0.0000000  0.8944272
[2,] -0.4472136 -1.665335e-16  1.457168e-16  0.8660254 -0.2236068
[3,] -0.4472136 -5.773503e-01  5.773503e-01 -0.2886751 -0.2236068
[4,] -0.4472136 -2.113249e-01 -7.886751e-01 -0.2886751 -0.2236068
[5,] -0.4472136  7.886751e-01  2.113249e-01 -0.2886751 -0.2236068

# 奇异分解性质
# svd的u矩阵 %*% svd的d矩阵的对角性值 %*% svd的v的转置矩阵 = 原矩阵
> s$u %*% diag(s$d) %*% t(s$v) #=m1
[,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# svd的u矩阵的转置矩阵 %*% svd的u矩阵 = svd的s矩阵的转置矩阵 %*% svd的v矩阵
> t(s$u) %*% s$u
[,1]          [,2]          [,3]          [,4]          [,5]
[1,]  1.000000e+00  0.000000e+00 -5.551115e-17  0.000000e+00 -2.775558e-17
[2,]  0.000000e+00  1.000000e+00 -2.775558e-17 -1.387779e-16  8.326673e-17
[3,] -5.551115e-17 -2.775558e-17  1.000000e+00  6.245005e-17 -6.938894e-17
[4,]  0.000000e+00 -1.387779e-16  6.245005e-17  1.000000e+00 -1.387779e-16
[5,] -2.775558e-17  8.326673e-17 -6.938894e-17 -1.387779e-16  1.000000e+00
> t(s$v) %*% s$v
[,1]          [,2]         [,3]          [,4]          [,5]
[1,]  1.000000e+00  0.000000e+00 5.551115e-17 -1.665335e-16  2.775558e-17
[2,]  0.000000e+00  1.000000e+00 8.326673e-17 -8.326673e-17  0.000000e+00
[3,]  5.551115e-17  8.326673e-17 1.000000e+00  1.595946e-16  2.775558e-17
[4,] -1.665335e-16 -8.326673e-17 1.595946e-16  1.000000e+00 -8.326673e-17
[5,]  2.775558e-17  0.000000e+00 2.775558e-17 -8.326673e-17  1.000000e+00

5. 特殊矩阵

下面介绍的多种特殊矩阵,都是在matrixcalc库中提供的。

5.1 Hankel Matrix
汉克尔矩阵 (Hankel Matrix) 是具有恒定倾斜对角线的方形矩阵。 Hankel矩阵的行列式称为catalecticant。该函数根据n向量x的值构造n阶Hankel矩阵。 矩阵的每一行是前一行中值的循环移位。

矩阵定义:


> hankel.matrix(6, 1:50)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    4    5    6
[2,]    2    3    4    5    6    7
[3,]    3    4    5    6    7    8
[4,]    4    5    6    7    8    9
[5,]    5    6    7    8    9   10
[6,]    6    7    8    9   10   11

5.2 Hilbert Matrix
希尔伯特矩阵是一种数学变换矩阵,正定,且高度病态(即,任何一个元素发生一点变动,整个矩阵的行列式的值和逆矩阵都会发生巨大变化),病态程度和阶数相关。希尔伯特矩阵是一种特殊的汉克尔矩阵,该函数返回n乘n希尔伯特矩阵。

矩阵定义:


> hilbert.matrix(4)
          [,1]      [,2]      [,3]      [,4]
[1,] 1.0000000 0.5000000 0.3333333 0.2500000
[2,] 0.5000000 0.3333333 0.2500000 0.2000000
[3,] 0.3333333 0.2500000 0.2000000 0.1666667
[4,] 0.2500000 0.2000000 0.1666667 0.1428571

5.3 Creation Matrix
创造矩阵,n阶创建矩阵也称为推导矩阵,该函数返回阶数n创建矩阵,在主对角线下方的子对角线上具有序列1,2,…,n-1的方阵。

矩阵定义:


> creation.matrix(5)  
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    1    0    0    0    0
[3,]    0    2    0    0    0
[4,]    0    0    3    0    0
[5,]    0    0    0    4    0

5.4 Stirling Matrix
斯特林公式(Stirling’s approximation)是一条用来取n的阶乘的近似值的数学公式。一般来说,当n很大的时候,n阶乘的计算量十分大,所以斯特林公式十分好用,而且,即使在n很小的时候,斯特林公式的取值已经十分准确。

斯特林矩阵(Stirling Matrix),该函数构造并返回斯特林矩阵,该矩阵是包含第二类斯特林数的下三角矩阵。

矩阵定义:


> stirling.matrix(4)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    0    0    0
[3,]    0    1    1    0    0
[4,]    0    1    3    1    0
[5,]    0    1    7    6    1

5.5 Pascal matrix
帕斯卡矩阵:由杨辉三角形表组成的矩阵称为帕斯卡(Pascal)矩阵。此函数返回n乘以Pascal矩阵。在数学中,尤其是矩阵理论和组合学,Pascal矩阵是一个下三角矩阵,行中有二项式系数。 通过对相同顺序的对称Pascal矩阵执行LU分解并返回下三角矩阵,可以容易地获得它。

帕斯卡的三角形是由数字行组成的三角形。 第一行具有条目1.每个后续行通过添加前一行的相邻条目而形成,替换为0,其中不存在相邻条目。 pascal函数通过选择与指定矩阵维度相对应的Pascal三角形部分来形成Pascal矩阵,

矩阵定义:


> pascal.matrix(4)
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    1    1    0    0
[3,]    1    2    1    0
[4,]    1    3    3    1

5.6 Fibonacci matrix
斐波纳契矩阵,该函数构造了从Fibonacci序列导出的n + 1平方Fibonacci矩阵。

计算公式:


> fibonacci.matrix(4)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    0    0    0
[2,]    2    1    1    0    0
[3,]    3    2    1    1    0
[4,]    5    3    2    1    1
[5,]    8    5    3    2    1

5.7 Frobenius Matrix
Frobenius矩阵也称为伴随矩阵,它出现在线性一阶微分方程组的解中。
此函数返回一个在数值数学中有用的Fronenius矩阵。

矩阵定义:


> frobenius.matrix(4)
     [,1] [,2] [,3] [,4]
[1,]    0    0    0   -1
[2,]    1    0    0    4
[3,]    0    1    0   -6
[4,]    0    0    1    4

5.8 Duplication matrix
复制矩阵,当A是对称矩阵时,该函数构造将vech(A)映射到vec(A)的线性变换D.

计算公式:


> D.matrix(3) 
      [,1] [,2] [,3] [,4] [,5] [,6]
 [1,]    1    0    0    0    0    0
 [2,]    0    1    0    0    0    0
 [3,]    0    0    1    0    0    0
 [4,]    0    1    0    0    0    0
 [5,]    0    0    0    1    0    0
 [6,]    0    0    0    0    1    0
 [7,]    0    0    1    0    0    0
 [8,]    0    0    0    0    1    0
 [9,]    0    0    0    0    0    1

5.9 K matrix
k矩阵是由H.matrices()函数构造的,利用直积进行计算子列表的分量。K.matrix(r, c = r),返回阶数为p = r * c的方阵,对于r行c列的矩阵A,计算A和t(A)的直积。

计算公式:


> K.matrix(2,3)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    0    0    0    0    0
[2,]    0    0    1    0    0    0
[3,]    0    0    0    0    1    0
[4,]    0    1    0    0    0    0
[5,]    0    0    0    1    0    0
[6,]    0    0    0    0    0    1

5.10 E Matrices
E矩阵列表,E.matrices(n)使得每个子列表的分量,是从n阶单位矩阵计算向量的外积导出的方阵。

计算公式:


> E.matrices(3)
[[1]]
[[1]][[1]]
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    0    0
[3,]    0    0    0

[[1]][[2]]
     [,1] [,2] [,3]
[1,]    0    1    0
[2,]    0    0    0
[3,]    0    0    0

[[1]][[3]]
     [,1] [,2] [,3]
[1,]    0    0    1
[2,]    0    0    0
[3,]    0    0    0


[[2]]
[[2]][[1]]
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    1    0    0
[3,]    0    0    0

[[2]][[2]]
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    1    0
[3,]    0    0    0

[[2]][[3]]
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    1
[3,]    0    0    0


[[3]]
[[3]][[1]]
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    1    0    0

[[3]][[2]]
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    1    0

[[3]][[3]]
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    1

5.11 H Matrices
H矩阵列表,H.matrices(r, c = r)使得r阶c阶的子列表的分量,计算从r行和c列的单位矩阵的列向量的外积导出的方阵。


> H.matrices(2,3)
[[1]]
[[1]][[1]]
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    0    0

[[1]][[2]]
     [,1] [,2] [,3]
[1,]    0    1    0
[2,]    0    0    0

[[1]][[3]]
     [,1] [,2] [,3]
[1,]    0    0    1
[2,]    0    0    0


[[2]]
[[2]][[1]]
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    1    0    0

[[2]][[2]]
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    1    0

[[2]][[3]]
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    1

5.12 T Matrices
T矩阵列表,T.matrices(n) 高级别列表中的组件数为n。 n个组件中的每一个也是列表。 每个子列表具有n个分量,每个分量是n阶矩阵。

计算公式:


> T.matrices(3)
[[1]]
[[1]][[1]]
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    0    0
[3,]    0    0    0

[[1]][[2]]
     [,1] [,2] [,3]
[1,]    0    1    0
[2,]    1    0    0
[3,]    0    0    0

[[1]][[3]]
     [,1] [,2] [,3]
[1,]    0    0    1
[2,]    0    0    0
[3,]    1    0    0


[[2]]
[[2]][[1]]
     [,1] [,2] [,3]
[1,]    0    1    0
[2,]    1    0    0
[3,]    0    0    0

[[2]][[2]]
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    1    0
[3,]    0    0    0

[[2]][[3]]
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    1
[3,]    0    1    0


[[3]]
[[3]][[1]]
     [,1] [,2] [,3]
[1,]    0    0    1
[2,]    0    0    0
[3,]    1    0    0

[[3]][[2]]
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    1
[3,]    0    1    0

[[3]][[3]]
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    1

通过R语言,龙虎大战坐庄龙虎大战坐庄我 们 实现了对于矩阵的各种计算操作,非常方便!有了好的龙虎大战坐庄工具 ,不管是学习还是应用,都会事半功倍。本文只是列举了矩阵的操作龙虎大战坐庄方法 ,没有解释计算的推到过程,推到过程请参考线性代码的教科书。

转载请注明出处:
http://cnc77.com/r-matrix

打赏作者

信息熵的时间度量模型muti

R的极客理想系列文章,涵盖了R的思想,使用,龙虎大战坐庄工具 ,创新等的一系列要点,以龙虎大战坐庄我 个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,龙虎大战坐庄互联网 ….都在使用R语言。

要成为有理想的极客,龙虎大战坐庄龙虎大战坐庄我 们 不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让龙虎大战坐庄龙虎大战坐庄我 们 一起动起来吧,开始R的极客理想。

龙虎大战坐庄关于 作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://cnc77.com
  • email: bsspirit@gmail.com

转载请注明出处:
http://cnc77.com/r-entropy-muti

前言

信息熵作为信息度量的基础,被广泛的算法模型所使用,比如决策树,熵值法等,基于信息熵的扩展龙虎大战坐庄龙虎大战坐庄我 们 可以做很多的事情。

本文将给大家介绍一个新模型思路,用信息熵来衡量具有时间特征的两随机变量的信息差。由Mark D. Scheuerell教授开发,所属美国西雅图西雅图国家海洋和大气管理局,国家海洋渔业局西北渔业科学中心,鱼类生态学部。

目录

  1. muti包介绍
  2. muti包函数定义
  3. muti包使用
  4. 实战案例

1. muti包介绍

muti包是用于计算两个离散随机变量之间的互信息(MI)的R包。 muti包设计目标,是考虑时间序列分析的情况下开发的,但没有任何龙虎大战坐庄方法 可以将龙虎大战坐庄方法 与时间索引本身联系起来。

muti项目的龙虎大战坐庄主页 :https://mdscheuerell.github.io/muti/

熵是用来含量信息的单位,龙虎大战坐庄我 在文章 用R语言实现信息度量 中已有详细介绍。那么,通过计算两随机变量的互信息,来估计一个变量中包含的另一个变量的信息量,可以理解为是两个变量之间协方差的非参数度量。

muti包中用来计算互信息的计算过程,包括以下几步:

第一步,计算信息熵H(x),H(y):

第二步,计算联合熵H(x,y):

第三步,计算互信息MI(x;y):

第四步,标准化互信息,映射到[0,1]区间:

开发环境所使用的系统环境

  • 龙虎大战坐庄Win 10 64bit
  • R: 3.4.2 x86_64-w64-mingw32/x64 b4bit

muti目前还没有发布到CRAN上面,需要通过devtools包从github上进行安装。安装操作如下。


> if(!require("devtools")) {
>  install.packages("devtools")
>  library("devtools")
> }
> devtools::install_github("mdscheuerell/muti")

2. muti包函数定义

muti包,有很强的专业性,没有太多的函数定义,函数都是为了实现模型的算法。

  • muti(),主函数,用于计算两个向量的互信息并画图输出
  • mutual_info(),计算互信息
  • newZ(),基于原始数据的重新采样创建新矢量。
  • symbolize(),将数字向量转换为符号向量。
  • transM(),估计用于计算互信息上的临界阈值的转换矩阵。
  • plotMI(),画图输出

2.1 muti()函数:
主函数,用于计算两个向量的互信息并画图输出。

函数定义:


muti(x, y, sym = TRUE, n_bins = NULL, normal = FALSE, lags = seq(-4, 4), mc = 100, alpha = 0.05)

参数列表:

  • x, 向量
  • y, 向量
  • sym,是否转换成符号向量,默认TRUE,
  • n_bins, 离散化龙虎大战坐庄方法 ,当sym=FALSE时, n_bins = ceiling(2 * length(x)^(1/3))
  • normal, 是否标准化,默认FALSE
  • lags, 滞后周期,默认为seq(-4, 4)
  • mc, 用于估计互信息的临界阈值的蒙特卡洛数量的模拟,默认100,
  • alpha, 用于定义临界阈值的取值点0.05

2.2 mutual_info()函数:
计算互信息,源函数的定义过于复杂,龙虎大战坐庄我 重新进行的龙虎大战坐庄优化 ,使得结构简单易懂。


# 源函数
mutual_info <- function(su,normal) {
  ## function to calculate mutual information between 2 time series
  ## need to remove obs with NA's b/c they don't contribute to info
  su <- su[!apply(apply(su,1,is.na),2,any),]
  ## get number of bins
  bb <- max(su)
  ## get ts length
  TT <- dim(su)[1]
  ## init entropies
  H_s <- H_u <- H_su <- 0
  for(i in 1:bb) {
    ## calculate entropy for 1st series
    p_s <- sum(su[,1] == i)/TT
    if(p_s != 0) { H_s <- H_s - p_s*log(p_s, 2) }
    ## calculate entropy for 2nd series
    p_u <- sum(su[,2] == i)/TT
    if(p_u != 0) { H_u <- H_u - p_u*log(p_u, 2) }
    for(j in 1:bb) {
      ## calculate joint entropy
      p_su <- sum(su[,1]==i & su[,2]==j)/TT
      if(p_su != 0) { H_su <- H_su - p_su*log(p_su, 2) }
    } ## end j loop
  } ## end i loop
  ## calc mutual info
  MI <- H_s + H_u - H_su
  if(!normal) { return(MI) } else { return(MI/sqrt(H_s*H_u)) }
}

# 修改后的函数
library(philentropy)
mutual_info <- function(su,normal){
  n<-nrow(su)
  px<-table(su[,1])/n
  py<-table(su[,2])/n
  
  hX<-H(px)
  hY<-H(py)
  
  tf<-table(su[,1],su[,2])
  pXY<-prop.table(tf)
  hXY<-JE(pXY)
  
  MI<-hX+hY-hXY
  if(normal){
    return(MI/sqrt(hX * hY))
  }else{
    return(MI)
  }
}

2.3 symbolize()函数:
将数字向量转换为符号向量,同样源函数定义过于复杂,龙虎大战坐庄我 重写了一下。


# 源函数
symbolize <- function(xy) {
  ## of interest and converts them from numeric to symbolic.
  ## check input for errors
  if(dim(xy)[2] != 2) {
    stop("xy must be an [n x 2] matrix \n\n", call.=FALSE)
  }
  ## get ts length
  TT <- dim(xy)[1]
  ## init su matrix
  su <- matrix(NA, TT, 2)
  ## convert to symbols
  ## loop over 2 vars
  for(i in 1:2) {
    for(t in 2:(TT-1)) {
      ## if xy NA, also assign NA to su
      if(any(is.na(xy[(t-1):(t+1),i]))) { su[t,i] <- NA }
      ## else get correct symbol
      else {
        if(xy[t,i] == xy[t-1,i] & xy[t,i] == xy[t+1,i]) {
          ## same
          su[t,i] <- 3
          }
        if(xy[t,i] > xy[t-1,i]) {
          ## peak
          if(xy[t,i] > xy[t+1,i]) { su[t,i] <- 5 }
          ## increase
          else { su[t,i] <- 4 }
        }
        if(xy[t,i] < xy[t-1,i]) {
          ## trough
          if(xy[t,i] < xy[t+1,i]) { su[t,i] <- 1 }
          ## decrease
          else { su[t,i] <- 2 }
        }
      } ## end else
    } ## end t loop
  } ## end i loop
  ## return su matrix
  return(su)
}

# 龙虎大战坐庄优化
后的函数
symbolize<-function(xy){
  if(ncol(xy)!=2) {
    stop("xy must be an [n x 2] matrix \n\n", call. = FALSE)
  }
  
  idx<-2:(nrow(xy)-1)
  apply(xy,2,function(row){
    na<-which(is.na(row[idx]+row[idx-1]+row[idx+1]))
    n3<-which(row[idx]==row[idx-1] & row[idx]==row[idx+1])
    n5<-which(row[idx]>row[idx-1] & row[idx]>row[idx+1])
    n4<-which(row[idx]>row[idx-1] & row[idx]<=row[idx+1])
    n1<-which(row[idx]<row[idx-1] & row[idx]<row[idx+1])
    n2<-which(row[idx]<row[idx-1] & row[idx]>=row[idx+1])

    line<-rep(0,length(idx))
    line[na]<-NA
    line[n1]<-1
    line[n2]<-2
    line[n3]<-3
    line[n4]<-4
    line[n5]<-5
    c(NA,line,NA)
  })
}

2.4 transM()函数:
估计用于计算互信息上的临界阈值的转换矩阵,同样源函数定义过于复杂,龙虎大战坐庄我 重写了一下。


# 源函数
transM <- function(x,n_bins) {
  ## helper function for estimating transition matrix used in
  ## creating resampled ts for the CI on mutual info
  ## replace NA with runif()
  x[is.na(x)] <- runif(length(x[is.na(x)]),min(x,na.rm=TRUE),max(x,na.rm=TRUE))
  ## length of ts
  tt <- length(x)
  ## get bins via slightly extended range
  bins <- seq(min(x)-0.001,max(x),length.out=n_bins+1)
  ## discretize ts
  hin <- vector("numeric",tt)
  for(b in 1:n_bins) {
    hin[x > bins[b] & x <= bins[b+1]] <- b
  }
  ## matrix of raw-x & discrete-x
  xn <- cbind(x,hin)
  ## transition matrix from bin-i to bin-j
  MM <- matrix(0,n_bins,n_bins)
  for(i in 1:n_bins) {
    for(j in 1:n_bins) {
      MM[i,j] <- sum(hin[-tt]==i & hin[-1]==j)
    }
    if(sum(MM[i,])>0) { MM[i,] <- MM[i,]/sum(MM[i,]) }
    else { MM[i,] <- 1/n_bins }
  }
  return(list(xn=xn,MM=MM))
} ## end function


# 龙虎大战坐庄优化
后的函数
transM <- function(x,n_bins){
  x[is.na(x)] <- runif(length(x[is.na(x)]), min(x, na.rm = TRUE),max(x, na.rm = TRUE))
  n<-length(x)
  xn<-cbind(x,hin=cut(x,7))

  hin<-xn[,2]
  a<-data.frame(id=1:(n-1),a=hin[1:(n-1)])
  b<-data.frame(id=1:(n-1),b=hin[2:n])
  ab<-merge(a,b,by="id",all=TRUE)
  ab$v<-1
  abm<-aggregate(v~a+b,data=ab,sum)
  
  MM <- matrix(0, n_bins, n_bins)
  apply(abm,1,function(row){
    MM[row[1],row[2]]<<-row[3]
    invisible()
  })
  
  MN<-apply(MM,1,function(row){
    s<-sum(row)
    if(s>0){
      row/s
    }else{
      rep(1/n_bins,n_bins)
    }
  })
  MN<-t(MN) 
  return(list(xn = xn, MM = MN))
}

由于作者的代码,是以it的逻辑写的,其实并不太符合数据计算的逻辑,龙虎大战坐庄我 用向量化代码对源代码进行了龙虎大战坐庄优化 ,使用算法逻辑更清晰一下,同时龙虎大战坐庄我 也上传到了github:https://github.com/bsspirit/muti

3. muti包使用

接下来,就是对muti包的使用了,作者提供了一个实例给龙虎大战坐庄龙虎大战坐庄我 们 做参考。

3.1 数据的输入和输出
muti()函数,是muti包的主函数,调用muti()就完成了这个包的全部功能。

  • 输入,要求是2个数字型向量。
  • 输出,包括了2个部分,1.数据计算结果,2.可视化的结果。

新建一个程序

> set.seed(123)
> TT <- 30

# 创建两个数字型向量
> x1 <- rnorm(TT)
> y1 <- x1 + rnorm(TT)

# 计算互信息
> muti(x1, y1)
  lag MI_xy MI_tv
1  -4 0.238 0.595
2  -3 0.473 0.532
3  -2 0.420 0.577
4  -1 0.545 0.543
5   0 0.712 0.584
6   1 0.393 0.602
7   2 0.096 0.546
8   3 0.206 0.573
9   4 0.403 0.625

数据计算结果,包括3列:

  • lag: 滞后期,当lag=-1,length(x)=length(y)=26时,x取值[1:25],y取值[2:26];当lag=1,length(x)=length(y)=26时,x取值[2:26],y取值[1:25];
  • MI_xy: 互信息
  • MI_tv: 相应显着性阈值,用蒙特卡洛龙虎大战坐庄方法 生成互信息,按从大到小排序后取1-alpha位置的值。


可视化的结果,包括3个图:

  • 上图,原始数据的x,y两个变量
  • 中图,符号转型后的x,y两个变量,把原始数据离散化为5个值,差异化更大
  • 下图,滞后期的互信息(实线)和相应显着性阈值(虚线),同数据计算结果的输出。

3.2 离散化处理

在计算信息熵时,muti包是将源数据进行了离散化处理,通过离散化进行分组,计算不同分组的概率矩阵,从而得出互信息。

muti包提供了2种离散化。

  • 符号化: 把连续型的数据,映射到5个值的上,分别对应peak(峰值), increasing(增加), same(没变),decreasing(减少),trough(低谷)。比如,向量c(1.1,2.1,3.3,1.2,3.1),不计算第一个值和最后一个值,第二个值2.1>第一个值1.1,第二个值2.1<第三个值3.3时,第二个值为增加。第2,3,4值,分别为c(“增加”,“峰值”,“低谷”)。
  • 分箱法: 等宽分箱,如果没有指定箱数,则根据莱斯规则计算,其中n = ceiling(2*length(x)^(1/3))。

4. 实战案例

利用muti的模型特点,龙虎大战坐庄龙虎大战坐庄我 们 取一组金融市场的数据进行分析,满足x和y是两个连续型变量,同时x和y之间有一些在金融市场上有一定的关系,通过分析x对y的确定性的贡献度来解释互信息值的合理性。

龙虎大战坐庄龙虎大战坐庄我 们 取中国A股市场上证综指(000001.SH)和工商银行(601398.SH)的每日收盘价的向前复权的数据,从2014-10-21到2018-12-31。工商银行是上证综指的权重股,所以让x=工商银行,y=上证综指,通过计算互信息来验证工商银行能够给上证综指带来多少确定性的信息增益。

整个案例部分的代码,有多处省略,并不影响阅读和分析。如果需要代码的朋友,请扫文章下面的二维码,请作者喝杯咖啡。

整理来xts类型的时间序列数据,数据格式如下:

// 省略部分代码
> head(stock,20)
         date X601398 X000001
1  2014-10-21    2.92 2339.66
2  2014-10-22    2.91 2326.55
3  2014-10-23    2.92 2302.42
4  2014-10-24    2.92 2302.28
5  2014-10-27    2.90 2290.44
6  2014-10-28    2.92 2337.87
7  2014-10-29    2.94 2373.03
8  2014-10-30    2.98 2391.08
9  2014-10-31    3.04 2420.18
10 2014-11-03    3.01 2430.03
11 2014-11-04    3.02 2430.68
12 2014-11-05    2.97 2419.25
13 2014-11-06    2.97 2425.86
14 2014-11-07    2.97 2418.17
15 2014-11-10    3.05 2473.67
16 2014-11-11    3.16 2469.67
17 2014-11-12    3.12 2494.48
18 2014-11-13    3.12 2485.61
19 2014-11-14    3.14 2478.82
20 2014-11-17    3.11 2474.01

3列数据,第一列是时间列,第二列是X601398人收盘价,第三列是X000001收盘价。从价格上看,2个数据集并不在同一个度量内,龙虎大战坐庄龙虎大战坐庄我 们 转化度量单位为每日收益率。

数据变型后,每日收益率结果如下:

> head(ret,20)
                X601398       X000001
2014-10-21           NA            NA
2014-10-22 -0.003424658 -5.603378e-03
2014-10-23  0.003436426 -1.037158e-02
2014-10-24  0.000000000 -6.080559e-05
2014-10-27 -0.006849315 -5.142728e-03
2014-10-28  0.006896552  2.070781e-02
2014-10-29  0.006849315  1.503933e-02
2014-10-30  0.013605442  7.606309e-03
2014-10-31  0.020134228  1.217023e-02
2014-11-03 -0.009868421  4.069945e-03
2014-11-04  0.003322259  2.674864e-04
2014-11-05 -0.016556291 -4.702388e-03
2014-11-06  0.000000000  2.732252e-03
2014-11-07  0.000000000 -3.170010e-03
2014-11-10  0.026936027  2.295124e-02
2014-11-11  0.036065574 -1.617031e-03
2014-11-12 -0.012658228  1.004588e-02
2014-11-13  0.000000000 -3.555851e-03
2014-11-14  0.006410256 -2.731724e-03
2014-11-17 -0.009554140 -1.940439e-03

然后,龙虎大战坐庄龙虎大战坐庄我 们 把每日收益率转换为累积收益率,让收益曲线为连续值,这样就构建好了,muti模型需要的输入数据。

> head(rdf,20)
                X601398       X000001
2014-10-22 -0.003424658 -0.0056033783
2014-10-23  0.000000000 -0.0159168426
2014-10-24  0.000000000 -0.0159766804
2014-10-27 -0.006849315 -0.0210372447
2014-10-28  0.000000000 -0.0007650684
2014-10-29  0.006849315  0.0142627561
2014-10-30  0.020547945  0.0219775523
2014-10-31  0.041095890  0.0344152569
2014-11-03  0.030821918  0.0386252703
2014-11-04  0.034246575  0.0389030885
2014-11-05  0.017123288  0.0340177633
2014-11-06  0.017123288  0.0368429601
2014-11-07  0.017123288  0.0335561577
2014-11-10  0.044520548  0.0572775531
2014-11-11  0.082191781  0.0555679030
2014-11-12  0.068493151  0.0661720079
2014-11-13  0.068493151  0.0623808588
2014-11-14  0.075342466  0.0594787277
2014-11-17  0.065068493  0.0574228734
2014-11-18  0.058219178  0.0498833164

取前50个点,把x和y带入muti()函数。


> x<-rdf$X601398
> y<-rdf$X000001
> result<-muti(x,y)
> result
  lag MI_xy MI_tv
1  -4 0.148 0.534
2  -3 0.158 0.464
3  -2 0.190 0.385
4  -1 0.351 0.425
5   0 0.646 0.398
6   1 0.528 0.340
7   2 0.262 0.461
8   3 0.207 0.437
9   4 0.260 0.518

可视化输出:

  • 上图,为原始数据的x(上证综指)和y(工商银行)的累积收益率曲线。
  • 中图,为符号转型后的x(上证综指)和y(工商银行)的累积收益率曲线。
  • 下图,符号转型后的数据,滞后期从-4到4的9个的互信息(实线)和和相应显着性阈值(虚线)。

当滞后期为lag=0时,MI互信息最大为0.646,表示工商银行给上证综指带来的信息是最确定的。而当有不同的滞后期时,MI互信息都为小,表示工商银行给上证综指的信息确定性变小了。

这个结果是符合金融市场的业务规则的。工商银行是上证综指的权重股,工商银行的当天股价的走势是会影响当天的上证综指的走势的,但是并不会影响上证综指前一天或第二天的走势,所以得到的结果是符合逻辑的。

另外,龙虎大战坐庄龙虎大战坐庄我 们 再引入一个统计变量进行验证。假设越2组数据越确定,那么2组数据的差值的标准差应该是越小,差异越大则标准差越大,然后再进行验证。

输出数据集为4列,最后一列增加了标准差(sd)。

// 省略部分代码
> rsdf
  lag MI_xy MI_tv       sd
1  -4 0.148 0.534 1.991886
2  -3 0.158 0.464 2.400721
3  -2 0.190 0.385 2.107189
4  -1 0.351 0.425 2.378768
5   0 0.646 0.398 1.477030
6   1 0.528 0.340 2.342482
7   2 0.262 0.461 1.997676
8   3 0.207 0.437 2.267211
9   4 0.260 0.518 2.116934

从数据可以看出,第5行当lag=0时,sd=1.477是最小的,说明这个时候,2组数据的差异是最小的,而其他的时候sd都比较大。龙虎大战坐庄龙虎大战坐庄我 们 再有可视化进行输出。

图中有3条线,黑实线为互信息,灰色线为阈值,红色线为标准差。这样就更直观地看出,标准化的曲线和互信息的曲线刚好相反。当差异越小时,标准差越小,互信息越大,x对y的确定性越大。这样就证明了,muti包也是适合用在金融市场的模型。

虽然本案例没有使用模型构建策略生成交易信号,但论证了模型的可行性和正确性。大家可以基于文章的思路,找到适合的场景进行策略应用,比如跨期套利或跨品种套利时,2种不同合约之间的确定性的依赖程度,来形成交易策略。

龙虎大战坐庄龙虎大战坐庄我 们 完成对muti包的解读,使用,测试,验证的全过程。通过新的模型思路,获得新的知识,才能在金融市场上有龙虎大战坐庄更多 的机会,从而突破信息相同,龙虎大战坐庄技术 龙虎大战坐庄工具 相同等的白热化竞争环境。

整个案例部分的代码,有多处省略,并不影响阅读和分析。如果需要代码的朋友,请扫文章下面的二维码,请作者喝杯咖啡。

转载请注明出处:
http://cnc77.com/r-entropy-muti

打赏作者

R语言统计特征描述包descriptr

R的极客理想系列文章,涵盖了R的思想,使用,龙虎大战坐庄工具 ,创新等的一系列要点,以龙虎大战坐庄我 个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,龙虎大战坐庄互联网 ….都在使用R语言。

要成为有理想的极客,龙虎大战坐庄龙虎大战坐庄我 们 不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让龙虎大战坐庄龙虎大战坐庄我 们 一起动起来吧,开始R的极客理想。

龙虎大战坐庄关于 作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://cnc77.com
  • email: bsspirit@gmail.com

转载请注明出处:
http://cnc77.com/r-descriptr

前言

龙虎大战坐庄龙虎大战坐庄我 们 获得数据后需要了解数据,通常会用到统计特征来观察数据,比如字段类型,数据集长度,均值,方差,数据分布,概率密度等。

descriptr包,为了龙虎大战坐庄龙虎大战坐庄我 们 提供了一套用来观察数据统计特征的龙虎大战坐庄工具 集,主要特性包括统计特征计算,离散度,频率表,交叉表,分组摘要和多个单/双向表的度量,可以让龙虎大战坐庄龙虎大战坐庄我 们 非常方便的观察的数据特征。

目录

  1. descriptr包介绍
  2. descriptr包函数列表
  3. descriptr包函数使用

1. descriptr包介绍

descriptr包,主要用于生成描述性统计信息。它提供了3种数据处理视角,连续变量、类别变量(离散变量)和可视化。descriptr包的统计特征计算部分的源代码,结构非常工整,大量用到dplyr包来构建。

开发环境所使用的系统环境

  • 龙虎大战坐庄Win 10 64bit
  • R: 3.4.2 x86_64-w64-mingw32/x64 b4bit

descriptr包的安装比较简单,直接用install.pacakges()函数就行。


> install.packages("descriptr")
> library(descriptr)

2. descriptr包函数列表

descriptr包,提供了3种处理视角,连续变量、类别变量、和可视化。龙虎大战坐庄龙虎大战坐庄我 们 将分别介绍这3种处理视角的函数。

descriptr包,提供了2个数据集,龙虎大战坐庄龙虎大战坐庄我 们 可以基于这2个数据集进行学习和测试。文中的例子,都是基于mtcarz的数据集进行构建的。

数据集:

  • hsb, 高中数据集
  • mtcarz, 汽车数据集,复制系统的mtcars数据集

2.1 连续变量
2.1.1 统计概览

  • ds_summary_stats, 统计概率
  • ds_auto_summary_stats, 自动统计概率
  • ds_group_summary, 分组描述性统计
  • ds_auto_group_summary, 自动描述性统计
  • ds_tidy_stats, 多变精简统计概率
  • ds_multi_stats,已弃用函数,用ds_tidy_stats()替代

2.1.2 统计特征计算

  • ds_mode, 计算众数
  • ds_extreme_obs, 计算极端值
  • ds_freq_cont, 计算频数
  • ds_freq_table, 计算频率分布表
  • ds_percentiles,计算分位数
  • ds_range, 计算宽度, max(x)-min(x)
  • ds_kurtosis, 计算峰度
  • ds_skewness, 计算偏度
  • ds_gmean, 计算几何平均值, prod(x)^(1/length(x))
  • ds_hmean, 计算谐波均值, length(x)/sum(sapply(x, function(x) {1/x} ))
  • ds_css, 计算修正平方和, sum((x1-mean)^2+(x2-mean)^2+…)
  • ds_mdev, 计算平均绝对差, sum( abs(x1-mean) + abs(x2-mean) + …)
  • ds_cvar, 计算变异系数, sd(x)/mean(x) * 100%
  • ds_std_error, 计算标准误差, sd(x)/(length(x)^0.5)
  • ds_tailobs,计算最大最小的多个值

2.1.3 度量特征

  • ds_measures_location,位置的度量,包括均值,中位数和众数
  • ds_measures_symmetry, 对称性的度量,包括峰度和偏度
  • ds_measures_variation,变异的度量,包括宽度,方差,标准差

2.1.4 其他函数

  • ds_rindex, 计算值的索引,同which
  • ds_screener, 以表格展示数据

2.2 类别变量

  • ds_twoway_table,计算双向表
  • ds_cross_table, 展示双向表
  • ds_auto_freq_table, 展示多个单向表
  • ds_auto_cross_table, 展示多个双向表
  • ds_tway_tables, 已弃用函数,用ds_auto_cross_table()替换
  • ds_oway_tables,已弃用函数,用ds_auto_freq_table()替换

2.3 可视化
2.3.1 画图函数

  • ds_plot_bar Generate bar plots
  • ds_plot_bar_grouped Generate grouped bar plots
  • ds_plot_bar_stacked Generate stacked bar plots
  • ds_plot_box_group Compare distributions
  • ds_plot_box_single Generate box plots
  • ds_plot_density Generate density plots
  • ds_plot_histogram Generate histograms
  • ds_plot_scatter Generate scatter plots

2.3.2 已弃用函数,调用vistributions包

  • dist_binom_perc, 可视化二项分布
  • dist_binom_plot, 可视化二项分布
  • dist_binom_prob,可视化二项分布
  • dist_chi_perc, 可视化卡方分布
  • dist_chi_plot, 可视化卡方分布
  • dist_chi_prob, 可视化卡方分布
  • dist_f_perc, 可视化F分布
  • dist_f_plot, 可视化F分布
  • dist_f_prob, 可视化F分布
  • dist_norm_perc, 可视化正态分布
  • dist_norm_plot, 可视化正态分布
  • dist_norm_prob, 可视化正态分布
  • dist_t_perc, 可视化T分布
  • dist_t_plot, 可视化T分布
  • dist_t_prob, 可视化T分布

2.4 演示小程序
一个演示的小程序,可以快速看到功能界面,使用shiny来构建的。

  • ds_launch_shiny_app, Shiny演示小程序

3. descriptr包函数使用

接下来,龙虎大战坐庄龙虎大战坐庄我 们 找一些对于龙虎大战坐庄龙虎大战坐庄我 们 观察数据非常方便的函数进行列举。

首先,龙虎大战坐庄龙虎大战坐庄我 们 先了解一个龙虎大战坐庄龙虎大战坐庄我 们 要使用的数据集mtcarz


> mtcarz
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive      21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Valiant             18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
Merc 240D           24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
Merc 230            22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
Merc 280            19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
Merc 280C           17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
Merc 450SL          17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
Fiat 128            32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
Honda Civic         30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
Toyota Corolla      33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
Toyota Corona       21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
Fiat X1-9           27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
Porsche 914-2       26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
Lotus Europa        30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
Ford Pantera L      15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
Ferrari Dino        19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
Maserati Bora       15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8
Volvo 142E          21.4   4 121.0 109 4.11 2.780 18.60  1  1    4    2

3.1 数据展示
通过ds_screener()函数进行静态数据集展示,替代函数原系统的str()函数。


# 查看数据静态结构
> ds_screener(mtcarz)
-----------------------------------------------------------------------
|  Column Name  |  Data Type  |  Levels   |  Missing  |  Missing (%)  |
-----------------------------------------------------------------------
|      mpg      |   numeric   |    NA     |     0     |       0       |
|      cyl      |   factor    |   4 6 8   |     0     |       0       |
|     disp      |   numeric   |    NA     |     0     |       0       |
|      hp       |   numeric   |    NA     |     0     |       0       |
|     drat      |   numeric   |    NA     |     0     |       0       |
|      wt       |   numeric   |    NA     |     0     |       0       |
|     qsec      |   numeric   |    NA     |     0     |       0       |
|      vs       |   factor    |    0 1    |     0     |       0       |
|      am       |   factor    |    0 1    |     0     |       0       |
|     gear      |   factor    |   3 4 5   |     0     |       0       |
|     carb      |   factor    |1 2 3 4 6 8|     0     |       0       |
-----------------------------------------------------------------------

 Overall Missing Values           0 
 Percentage of Missing Values     0 %
 Rows with Missing Values         0 
 Columns With Missing Values      0 

# str()函数的静态结构
> str(mtcarz)
'data.frame':	32 obs. of  11 variables:
 $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
 $ cyl : Factor w/ 3 levels "4","6","8": 2 2 1 2 3 2 3 1 1 2 ...
 $ disp: num  160 160 108 258 360 ...
 $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
 $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
 $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
 $ qsec: num  16.5 17 18.6 19.4 17 ...
 $ vs  : Factor w/ 2 levels "0","1": 1 1 2 2 1 2 1 2 2 2 ...
 $ am  : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
 $ gear: Factor w/ 3 levels "3","4","5": 2 2 2 1 1 1 1 2 2 2 ...
 $ carb: Factor w/ 6 levels "1","2","3","4",..: 4 4 1 1 2 1 4 2 2 4 ... 

3.2 统计概览
通过ds_summary_stats()函数,查看数据集中某个连续型变量的所有统计特征值。


# 统计概览
> ds_summary_stats(mtcarz,mpg)
-------------------------------------------- Variable: mpg --------------------------------------------

                        Univariate Analysis                          
 N                       32.00      Variance                36.32 
 Missing                  0.00      Std Deviation            6.03 
 Mean                    20.09      Range                   23.50 
 Median                  19.20      Interquartile Range      7.38 
 Mode                    10.40      Uncorrected SS       14042.31 
 Trimmed Mean            19.95      Corrected SS          1126.05 
 Skewness                 0.67      Coeff Variation         30.00 
 Kurtosis                -0.02      Std Error Mean           1.07 

                              Quantiles                               
              Quantile                            Value                
             Max                                  33.90                
             99%                                  33.44                
             95%                                  31.30                
             90%                                  30.09                
             Q3                                   22.80                
             Median                               19.20                
             Q1                                   15.43                
             10%                                  14.34                
             5%                                   12.00                
             1%                                   10.40                
             Min                                  10.40                

                            Extreme Values                            
                Low                                High                
  Obs                        Value       Obs                        Value 
  15                         10.4        20                         33.9  
  16                         10.4        18                         32.4  
  24                         13.3        19                         30.4  
   7                         14.3        28                         30.4  
  17                         14.7        26                         27.3  

输出分成了3个部分:Univariate Analysis(单变量分析),Quantiles(分位数),Extreme Values(极值)。

  • Univariate Analysis(单变量分析),包括N(个数),Missing(缺失值),Mean(均值),Median(中位数),Mode(众数),Trimmed Mean(修正均值),Skewness(偏度),Kurtosis(峰度),Variance(方差),Std Deviation(标准差),Range(范围,最大-最小),Interquartile Range(四分位数范围),Uncorrected SS(未修正平方和),Corrected SS(修正平方和), Coeff Variation(变异系数,标准差/均值),Std Error Mean(标准误差均值)
  • Quantiles(分位数),从最小值到最小值,按顺序排列,对应的数值。
  • Extreme Values(极值),包括最小值前5个,最大值前5个。

3.3 统计特征快速查看
通过ds_tidy_stats()函数,查看数据集中各变量的统计特征,维度比较少。


# 多变量统计
> ds_tidy_stats(mtcarz, mpg, disp, hp)
# A tibble: 3 x 16
  vars    min   max  mean t_mean median  mode range variance  stdev  skew kurtosis coeff_var
  <chr> <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl>    <dbl>  <dbl> <dbl>    <dbl>     <dbl>
1 disp   71.1 472   231.   228    196.  276.  401.   15361.  124.   0.420  -1.07        53.7
2 hp     52   335   147.   144.   123   110   283     4701.   68.6  0.799   0.275       46.7
3 mpg    10.4  33.9  20.1   20.0   19.2  10.4  23.5     36.3   6.03 0.672  -0.0220      30.0
# ... with 3 more variables: q1 <dbl>, q3 <dbl>, iqrange <dbl>

3.4 频率表
通过ds_freq_table()函数,把数据集中某个连续型变量,进行等宽划分,形成频率表。


# 划分成5个等宽的频率
> ds_freq_table(mtcarz,mpg,5)
                              Variable: mpg                               
|-----------------------------------------------------------------------|
|    Bins     | Frequency | Cum Frequency |   Percent    | Cum Percent  |
|-----------------------------------------------------------------------|
| 10.4 - 15.1 |     6     |       6       |    18.75     |    18.75     |
|-----------------------------------------------------------------------|
| 15.1 - 19.8 |    12     |      18       |     37.5     |    56.25     |
|-----------------------------------------------------------------------|
| 19.8 - 24.5 |     8     |      26       |      25      |    81.25     |
|-----------------------------------------------------------------------|
| 24.5 - 29.2 |     2     |      28       |     6.25     |     87.5     |
|-----------------------------------------------------------------------|
| 29.2 - 33.9 |     4     |      32       |     12.5     |     100      |
|-----------------------------------------------------------------------|
|    Total    |    32     |       -       |    100.00    |      -       |
|-----------------------------------------------------------------------|

3.5 分组统计
通过ds_group_summary()函数,把数据集中变量进行分组,再分别计算统计特征。


> k<-ds_group_summary(mtcarz,cyl,mpg);k
                                       mpg by cyl                                         
-----------------------------------------------------------------------------------------
|     Statistic/Levels|                    4|                    6|                    8|
-----------------------------------------------------------------------------------------
|                  Obs|                   11|                    7|                   14|
|              Minimum|                 21.4|                 17.8|                 10.4|
|              Maximum|                 33.9|                 21.4|                 19.2|
|                 Mean|                26.66|                19.74|                 15.1|
|               Median|                   26|                 19.7|                 15.2|
|                 Mode|                 22.8|                   21|                 10.4|
|       Std. Deviation|                 4.51|                 1.45|                 2.56|
|             Variance|                20.34|                 2.11|                 6.55|
|             Skewness|                 0.35|                -0.26|                -0.46|
|             Kurtosis|                -1.43|                -1.83|                 0.33|
|       Uncorrected SS|              8023.83|              2741.14|              3277.34|
|         Corrected SS|               203.39|                12.68|                 85.2|
|      Coeff Variation|                16.91|                 7.36|                16.95|
|      Std. Error Mean|                 1.36|                 0.55|                 0.68|
|                Range|                 12.5|                  3.6|                  8.8|
|  Interquartile Range|                  7.6|                 2.35|                 1.85|
-----------------------------------------------------------------------------------------

3.6 分组分类统计
通过ds_auto_group_summary()函数,把数据集中变量进行分组,再分别两两计算统计特征。


# 分组分类
> ds_auto_group_summary(mtcarz, cyl, gear, mpg)
                                       mpg by cyl                                         
-----------------------------------------------------------------------------------------
|     Statistic/Levels|                    4|                    6|                    8|
-----------------------------------------------------------------------------------------
|                  Obs|                   11|                    7|                   14|
|              Minimum|                 21.4|                 17.8|                 10.4|
|              Maximum|                 33.9|                 21.4|                 19.2|
|                 Mean|                26.66|                19.74|                 15.1|
|               Median|                   26|                 19.7|                 15.2|
|                 Mode|                 22.8|                   21|                 10.4|
|       Std. Deviation|                 4.51|                 1.45|                 2.56|
|             Variance|                20.34|                 2.11|                 6.55|
|             Skewness|                 0.35|                -0.26|                -0.46|
|             Kurtosis|                -1.43|                -1.83|                 0.33|
|       Uncorrected SS|              8023.83|              2741.14|              3277.34|
|         Corrected SS|               203.39|                12.68|                 85.2|
|      Coeff Variation|                16.91|                 7.36|                16.95|
|      Std. Error Mean|                 1.36|                 0.55|                 0.68|
|                Range|                 12.5|                  3.6|                  8.8|
|  Interquartile Range|                  7.6|                 2.35|                 1.85|
-----------------------------------------------------------------------------------------

                                       mpg by gear                                        
-----------------------------------------------------------------------------------------
|     Statistic/Levels|                    3|                    4|                    5|
-----------------------------------------------------------------------------------------
|                  Obs|                   15|                   12|                    5|
|              Minimum|                 10.4|                 17.8|                   15|
|              Maximum|                 21.5|                 33.9|                 30.4|
|                 Mean|                16.11|                24.53|                21.38|
|               Median|                 15.5|                 22.8|                 19.7|
|                 Mode|                 10.4|                   21|                   15|
|       Std. Deviation|                 3.37|                 5.28|                 6.66|
|             Variance|                11.37|                27.84|                44.34|
|             Skewness|                -0.09|                  0.7|                 0.56|
|             Kurtosis|                -0.38|                -0.77|                -1.83|
|       Uncorrected SS|              4050.52|               7528.9|              2462.89|
|         Corrected SS|               159.15|               306.29|               177.37|
|      Coeff Variation|                20.93|                21.51|                31.15|
|      Std. Error Mean|                 0.87|                 1.52|                 2.98|
|                Range|                 11.1|                 16.1|                 15.4|
|  Interquartile Range|                  3.9|                 7.08|                 10.2|
-----------------------------------------------------------------------------------------

3.7 测量
通过ds_measures_xxx()的几个函数,把数据集中变量,分别进行不同维度的统计特征。如果您想要查看位置,变化,对称性,百分位数和极端观测值的度量,请使用以下函数。 除了ds_extreme_obs()之外,所有这些都将使用单个或多个变量。 如果未指定变量,则它们将返回数据集中所有连续变量的结果。

数据集变化分析:范围,四分位范围,方差,标准差,变异系数,标准误差


> ds_measures_variation(mtcarz)
# A tibble: 6 x 7
  var    range     iqr  variance      sd coeff_var std_error
  <chr>  <dbl>   <dbl>     <dbl>   <dbl>     <dbl>     <dbl>
1 disp  401.   205.    15361.    124.         53.7   21.9   
2 drat    2.17   0.840     0.286   0.535      14.9    0.0945
3 hp    283     83.5    4701.     68.6        46.7   12.1   
4 mpg    23.5    7.38     36.3     6.03       30.0    1.07  
5 qsec    8.40   2.01      3.19    1.79       10.0    0.316 
6 wt      3.91   1.03      0.957   0.978      30.4    0.173 

数据集数值分析:均值,修正均值,中位数,众数


> ds_measures_location(mtcarz)
# A tibble: 6 x 5
  var     mean trim_mean median   mode
  <chr>  <dbl>     <dbl>  <dbl>  <dbl>
1 disp  231.      228    196.   276.  
2 drat    3.60      3.58   3.70   3.07
3 hp    147.      144.   123    110   
4 mpg    20.1      20.0   19.2   10.4 
5 qsec   17.8      17.8   17.7   17.0 
6 wt      3.22      3.20   3.32   3.44

数据集分位数分析:从最小值到最大值排序


> ds_percentiles(mtcarz)
# A tibble: 6 x 12
  var     min  per1  per5 per10     q1 median     q3  per95  per90  per99    max
  <chr> <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
1 disp  71.1  72.5  77.4  80.6  121.   196.   326    449    396.   468.   472   
2 drat   2.76  2.76  2.85  3.01   3.08   3.70   3.92   4.31   4.21   4.78   4.93
3 hp    52    55.1  63.6  66     96.5  123    180    254.   244.   313.   335   
4 mpg   10.4  10.4  12.0  14.3   15.4   19.2   22.8   31.3   30.1   33.4   33.9 
5 qsec  14.5  14.5  15.0  15.5   16.9   17.7   18.9   20.1   20.0   22.1   22.9 
6 wt     1.51  1.54  1.74  1.96   2.58   3.32   3.61   5.29   4.05   5.40   5.42

极值分析


> ds_extreme_obs(mtcarz,mpg)
# A tibble: 10 x 3
   type  value index
   <chr> <dbl> <int>
 1 high   33.9    20
 2 high   32.4    18
 3 high   30.4    19
 4 high   30.4    28
 5 high   27.3    26
 6 low    10.4    15
 7 low    10.4    16
 8 low    13.3    24
 9 low    14.3     7
10 low    14.7    17

3.8 类别变量频率表
通过ds_cross_table()函数,查看数据集中类别变量的双向表。


> ds_cross_table(mtcarz, cyl, gear)
    Cell Contents
 |---------------|
 |     Frequency |
 |       Percent |
 |       Row Pct |
 |       Col Pct |
 |---------------|

 Total Observations:  32 

----------------------------------------------------------------------------
|              |                           gear                            |
----------------------------------------------------------------------------
|          cyl |            3 |            4 |            5 |    Row Total |
----------------------------------------------------------------------------
|            4 |            1 |            8 |            2 |           11 |
|              |        0.031 |         0.25 |        0.062 |              |
|              |         0.09 |         0.73 |         0.18 |         0.34 |
|              |         0.07 |         0.67 |          0.4 |              |
----------------------------------------------------------------------------
|            6 |            2 |            4 |            1 |            7 |
|              |        0.062 |        0.125 |        0.031 |              |
|              |         0.29 |         0.57 |         0.14 |         0.22 |
|              |         0.13 |         0.33 |          0.2 |              |
----------------------------------------------------------------------------
|            8 |           12 |            0 |            2 |           14 |
|              |        0.375 |            0 |        0.062 |              |
|              |         0.86 |            0 |         0.14 |         0.44 |
|              |          0.8 |            0 |          0.4 |              |
----------------------------------------------------------------------------
| Column Total |           15 |           12 |            5 |           32 |
|              |        0.468 |        0.375 |        0.155 |              |
----------------------------------------------------------------------------

3.9 类别变量的双向表
通过ds_twoway_table()函数,查看数据集中类别变量的分组后的情况。


> ds_twoway_table(mtcarz, cyl, gear)
Joining, by = c("cyl", "gear", "count")
# A tibble: 8 x 6
  cyl   gear  count percent row_percent col_percent
  <fct> <fct> <int>   <dbl>       <dbl>       <dbl>
1 4     3         1  0.0312      0.0909      0.0667
2 4     4         8  0.25        0.727       0.667 
3 4     5         2  0.0625      0.182       0.4   
4 6     3         2  0.0625      0.286       0.133 
5 6     4         4  0.125       0.571       0.333 
6 6     5         1  0.0312      0.143       0.2   
7 8     3        12  0.375       0.857       0.8   
8 8     5         2  0.0625      0.143       0.4   

3.10 可视化连续型数据
分别以柱状图,密度图,分箱图,散点图,对连续型数据进行可视化,从左到右的4个图。


> ds_plot_histogram(mtcarz, mpg, disp)
> ds_plot_density(mtcarz, mpg, disp)
> ds_plot_box_single(mtcarz, mpg, disp)
> ds_plot_scatter(mtcarz, mpg, disp)

3.11 可视化类别型数据
分别以bar图对类别型数据可视化,从左到右的4个图。


> ds_plot_bar(mtcarz,cyl, gear)
> ds_plot_bar_stacked(mtcarz, cyl, gear)
> ds_plot_bar_grouped(mtcarz, cyl, gear)
> ds_plot_box_group(mtcarz, cyl, gear, mpg)

3.12 可视化分布图

5种统计分布的可视化效果,由于使用时提示已弃用,改为调用vistributions包的对应函数,所以大家可以改用vistributions包。

二项分布


> dist_binom_prob(10, 0.3, 4, type = 'exact')

卡方分布


> dist_chi_perc(0.22, 13, 'upper')

F分布


> dist_f_perc(0.125, 9, 35, 'upper')

正态分布


> dist_norm_perc(0.95, mean = 2, sd = 1.36, type = 'both')

T分布


> dist_t_prob(1.445, 7, 'interval')

3.13 启动shiny小程序

提供了一个界面,方便小白进行操作,其实没什么用。>_<

本文对于descriptr包进行的完整的介绍,descriptr主要用于统计特征的快速查看,一个方便的龙虎大战坐庄工具 包,对于初识数据集是非常有龙虎大战坐庄帮助 的。

转载请注明出处:
http://cnc77.com/r-descriptr

打赏作者

2019中国R大会龙虎大战坐庄北京 : 凯利公式-用胜率和赔率量化龙虎大战坐庄你 的投资

跨界知识聚会系列文章,“知识是用来分享和传承的”,各种会议、龙虎大战坐庄论坛 、沙龙都是分享知识的绝佳场所。龙虎大战坐庄我 也有幸作为演讲嘉宾参加了一些国内的大型会议,向大家展示龙虎大战坐庄我 所做的一些成果。从听众到演讲感觉是不一样的,把知识分享出来,龙虎大战坐庄你 才能收获龙虎大战坐庄更多 。

龙虎大战坐庄关于 作者

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://cnc77.com/meeting-r-20190526

前言

本次R语言(龙虎大战坐庄北京 )大会由统计之都主办,主题围绕着“数据科学”的展开,1个主会场和19个分会场,涉及70多个分享主题。

本次大会结合了学术界、产业界,跨领域跨专业的交流,通过碰撞形成了新的火花,会议覆盖统计学、大数据、人工智能相关理论及其在各行各业的具体应用,包括医疗健康、生物信息、医疗大数据、心理学、量化金融、工业工程、智能制造、龙虎大战坐庄软件 龙虎大战坐庄工具 、计算平台、概率统计、统计理论、机器学习、人工智能、大数据应用、自然语言、新闻传播、社交网络、商务统计、人文科学等数据科学话题。

目录

  1. 龙虎大战坐庄我 分享的主题:凯利公式-用胜率和赔率量化龙虎大战坐庄你 的投资
  2. 会议体验和照片分享

1. 龙虎大战坐庄我 分享的主题:凯利公式-用胜率和赔率量化龙虎大战坐庄你 的投资

本次大会龙虎大战坐庄我 被安排在量化金融专场,从理论到落地,有公式有代码。以“故事”为切入点做为开场,用简单的话语把复杂的金融逻辑讲明白。龙虎大战坐庄我 的分享的PPT龙虎大战坐庄下载 ,本次活动的官方介绍的链接

龙虎大战坐庄我 的分享是最后一天,在人民龙虎大战坐庄大学 第三教学楼教室,100人的会场完全坐满。 龙虎大战坐庄我 分享主题:凯利公式-用胜率和赔率量化龙虎大战坐庄你 的投资,核心内容基于龙虎大战坐庄我 写的一篇文章,用R语言解读凯利公式

龙虎大战坐庄我 主要为分三个部分进行介绍:

  • 故事开始:通过一个故事,带个场景。
  • 凯利公式:原理,公式,变型。
  • 赌局最优解:通过公式结合真实市场,定义不同场景下最优解。
  • 让时间帮龙虎大战坐庄龙虎大战坐庄我 们 龙虎大战坐庄赚钱 :用R语言程序,按假设条件进行仿真。

由于时间有限,还请有问题的朋友,给龙虎大战坐庄我 留言再进行沟通。

2. 会议体验和照片分享

本次大会场次真多,涉及的内容真是保罗万象,不仅是R语言相关,还有很多行为应用,极大地拓宽的R语言会议的视野。

2.1 会议体验证和总结
龙虎大战坐庄我 的分享依然保持着一惯风格,有理论,有公式,有代码,有落地,会后总体的反馈还是不错的,大家纷纷表示都听懂了,金融模型很有意思。

龙虎大战坐庄我 的分享。

满座的会场。

龙虎大战坐庄我 要招人

2.2 相关照片

会场介绍

会后的大家交流

李孟育博士,南华期货

部分工作人员大合照:

新朋友,老朋友,跨界的朋友,龙虎大战坐庄台湾 的朋友,用过R语言让龙虎大战坐庄龙虎大战坐庄我 们 有共同的话题,感谢人民龙虎大战坐庄大学 和统计之都,所有主办方的小伙们,辛苦了!明年再见。

转载请注明出处:
http://cnc77.com/meeting-r-20190526

打赏作者