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

用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-test-x2

前言

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

本文是第三篇卡方检验,卡方检验关注点在实际值与预测值的差距,对于预测类的模型给出了客观地评价龙虎大战坐庄方法 。

目录

  1. 卡方检验介绍
  2. 数据集
  3. 四格表法推导过程
  4. 程序实现

1. 卡方检验介绍

卡方检验,又称χ2检验,是一种非参数检验,主要是比较两个以及两个以上样本率以及两个分类变量之间是否具有显著的相关性,其根本思想是统计样本的实际观测值与理论推断值之间的偏离程度。卡方检验,是由英国统计学家Karl Pearson在1900年首次提出的,在《哲学杂志》上发表。

卡方检验,实际观测值与理论推断值之间的偏离程度就决定卡方值的大小,如果卡方值越大,两个变量偏差程度越大;反之,卡方值越小两个变量偏差越小;如果两个变量的完全相等,卡方值就为0,表明观测值与理论值完全符合。

通常卡方检验的主要应用方向为:卡方拟合优度检验,卡方独立性检验。卡方检验可以分为成组比较(不配对资料)和个别比较(配对,或同一对象两种处理的比较)两类。卡方检验试用条件包括,随机样本数据,卡方检验的理论频数不能太小。

卡方检验的计算公式:


公式解释:

  • χ2,为统计量,用于衡量实际值与理论值的差异程度。
  • A,为实际值
  • T,为理论值
  • 自由度,df = n-1

卡方检验有3种推导过程:四格表法的卡方检验,行列表法的卡方检验,列联表法的卡方检验,本文将介绍四格表法的卡方检验。

2. 数据集

卡方检验,对于数据有比较严格的要求,所以龙虎大战坐庄龙虎大战坐庄我 们 需要先找到一个合适的数据集,作为测试数据集。

龙虎大战坐庄龙虎大战坐庄我 们 以“抛硬币”来作为卡方检验的使用场景。对于1次抛硬币来说,只有2种结果,正面是字,反面是头像,两种结果的出现概率相同都是50%,且每次抛硬币都是独立事件。那么,根据这个抛硬币的结果做出假设,原假设(H0)抛硬币出现正面和反面的概率都是50%,为验证这个假设成立,就是卡方检验。

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

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

接下来,用随机过程模拟抛硬币的过程,设置抛硬币100次,1为正面,0为反面。


> N=100                 # 抛硬币100次
> set.seed(1)           # 随机种子
> coin<-sample(x=c(0,1), prob=c(0.5,0.5),size = N, replace = TRUE)  # 模拟抛硬币
> coin                  # 打印结果
  [1] 1 1 0 0 1 0 0 0 0 1 1 1 0 1 0 1 0 0 1 0 0 1 0 1 1 1 1 1 0 1 1 0 1 1 0 0 0 1 0 1 0 0 0 0
 [45] 0 0 1 1 0 0 1 0 1 1 1 1 1 0 0 1 0 1 1 1 0 1 1 0 1 0 1 0 1 1 1 0 0 1 0 0 1 0 1 1 0 1 0 1
 [89] 1 1 1 1 0 0 0 0 1 1 0 0

统计抛硬币的结果


> table(coin)  # 统计0和1的出现次数
coin
 0  1 
48 52 

用可视化展示累积的正面出现概率


> cnt <- cumsum(coin)  # 积累求和,正面1出现的次数
> n<-1:N               # 次数向量
> rate<- cnt/1:N       # 累积的正面出现概率
> plot(n,rate,type="o",log="x",   # 画图,x轴做log变型
+      xlim=c(1,N), ylim=c(0.0,1.0), cex.axis=1.5,
+      main="Running Proportion of Heads",cex.axis=1.5)

接下来,龙虎大战坐庄龙虎大战坐庄我 们 就以coin作为数据集,进行卡方检验的数据集。

3. 四格表法推导过程

四格表法的卡方检验,用于进行两个样本率或两个样本构成比的比较。

龙虎大战坐庄龙虎大战坐庄我 们 可以把抛硬币的过程,分成观察组和期望组,观察组就是上面随机过程的结果,期望组是龙虎大战坐庄龙虎大战坐庄我 们 根据概率的公式给出来的,把数据填入下面表格。

正面 反面 合计
观察次数 52 48 100
期望次数 50 50 100
合计 102 98 200

通过简单的计算,龙虎大战坐庄龙虎大战坐庄我 们 得出观察组的正面概率为52%和期望组的正面概率为50%,两者的差别可能是误差导致,也有可能是受到抛硬币的次数影响的。

为了确定真实原因,龙虎大战坐庄龙虎大战坐庄我 们 先假设抛硬币是符合大数定理的,是不受抛硬币的次数影响的,即抛硬币的次数和概率无关,龙虎大战坐庄龙虎大战坐庄我 们 可以得出正面结果的实际是(52+50) / (52+48+50+50)= 0.51 = 51%。

可推到出:

正面 反面 合计
观察次数 100*0.51 =51 100*(1-0.51) = 49 100
期望次数 100*0.51 =51 100*(1-0.51) = 49 100
合计 102 98 200

如果抛硬币次数与正面结果的概率是独立无关的,那么四格表里的理论值和实际值差别应该会很小。

代入公式进行计算。


X^2 = (52- 51)^2 / 51 + (48 - 49)^2 / 49 + (50 - 51)^2 / 51 + (50 - 49)^2 / 49 = 0.08003201

当以四格表法进行卡方检验时,还有一些限制条件,当两个独立样本比较时,需要分成以下3种情况:

  1. 所有的理论数T≥5并且总样本量n≥40,用Pearson卡方进行检验.
  2. 如果理论数T<5但T≥1,并且n≥40,用连续性校正的卡方进行检验.
  3. 如果有理论数T<1或n<40,则用Fisher’s检验.

根据限制条件,由于样本数据量为200大于40,所以需要对原公式进行修正。

修正公式为:


X^2 = (abs(52- 51)-0.5)^2 / 51 + (abs(48 - 49)-0.5)^2 / 49 + (abs(50 - 51)-0.5)^2 / 51 + (abs(50 - 49)-0.5)^2 / 49 
= 0.020008

最后,算出卡方统计量为0.020008。下面龙虎大战坐庄龙虎大战坐庄我 们 再用R语言的程序来算一下,看看是否是相同的结果。

4. 程序实现

用程序实现卡方检验,可以直接使用chisq.test()函数。


# 合并观测值和期望值,形成矩阵
> s<-c(table(coin),c(50,50))
> x <- matrix(s, ncol = 2, byrow=TRUE)
> x
     [,1] [,2]
[1,]   48   52
[2,]   50   50

# 卡方检验
> chisq.test(x)

	Pearson's Chi-squared test with Yates' continuity correction

data:  x
X-squared = 0.020008, df = 1, p-value = 0.8875

指标解释:

  • H0:抛硬币出现正面和反面的概率都是50%,与次数无关
  • X-squared统计量:0.020008
  • df,自由度,1
  • p-value值:0.8875

结果解读,以0.05为显著性水平,自由度df=1,统计量X-squared = 0.020008 小于临界值0.455(查表),说明X-squared值不显著,不能拒绝原假设。以0.05为显著性水平,p-value=0.8875大于0.05,所以也拒绝原假设,说明抛硬币的次数与结果为正面的概率是没有关系的。这个结果与龙虎大战坐庄龙虎大战坐庄我 们 构造的数据是一致的,也是符合大数定理的。

查看观察值和期望值。


# 查看观察值
> chi$observed
     [,1] [,2]
[1,]   48   52
[2,]   50   50

# 查看期望值
> chi$expected
     [,1] [,2]
[1,]   49   51
[2,]   49   51

# 残差
> chi$residuals
           [,1]      [,2]
[1,] -0.1428571  0.140028
[2,]  0.1428571 -0.140028

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


# 手动计算X^2值
> sr <- rowSums(x)
> sc <- colSums(x)
> n <- sum(x)
> E <- outer(sr, sc, "*")/n;E
     [,1] [,2]
[1,]   49   51
[2,]   49   51
> x2 <- sum((abs(x - E) - 0.5)^2/E);x2
[1] 0.020008

# 手动计算P值
> pchisq(x2, 1, lower.tail = FALSE)
[1] 0.8875147

本文用R语言详细介绍了卡方检验的计算过程,x^2越小,越能验证假设是正确的,越能说明观测值与期望值是一致的,x^2越大,越证明假设是错误的,预测越不准。为预测类的模型提供了评价标准。

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

打赏作者

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

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-f

前言

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

本文是第二篇F检验,T检验关注点在均值,而F检验关注点在方差。

目录

  1. F检验介绍
  2. 数据集
  3. F检验实现

1. F检验介绍

F检验法(F-test),初期叫方差比率检验(Variance Ratio),又叫联合假设检验(Joint Hypotheses Test),是英国统计学家Fisher提出的,主要通过比较两组数据的方差,以确定他们的密度是否有显著性差异。至于两组数据之间是否存在系统误差,则在进行F检验并确定它们的密度没有显著性差异之后,再进行T检验

F检验,是一种在零假设(H0)之下,统计值服从F-分布的检验。

样本标准偏差的平方公式:

F统计量计算公式:


公式解释

  • F:统计量,根据自由度查表,当F值小于查表值时没有显著差异,当F值大于等于查表值时有显著差异
  • S1:样本1的标准差
  • S2:样本2的标准差
  • 分子自由度: df=分子的数量-1
  • 分母自由度: df=分母的数量-1

T检验和F检验对比
T检验用来检测数据的准确度(系统误差),F检验用来检测数据的精密度(偶然误差)。在定量分析过程中,常遇到两种情况:一种是样本测量的平均值与真值不一致;另一种是两组测量的平均值不一致。

上述不一致是由于定量分析中的系统误差和偶然误差引起的,因此必须对两组分析结果的准确度或精密度是否存在显著性差异做出判断,两组数据的显著性检验顺序是先F检验后T检验。

T检验是检查两组均值的差异,而F检验是检查多组均值之间的差异。

对于多元线性回归模型,t检验是对于单个变量进行显著性,检验该变量独自对被解释变量的影响。f检验是检验回归模型的显著意义,即所有解释变量联合起来对被解释变量的影响,龙虎大战坐庄关于 线性回归请参考文章,R语言解读一元线性回归模型R语言解读多元线性回归模型

2. 数据集

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

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

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

数据集ToothGrowth,记录了60只豚鼠的牙齿生长速度,使用2种不同的龙虎大战坐庄方法 (OJ和VC),每天按3种不同的注射剂量进行注射,对牙齿的生长速度的对比数据,共3列,60条记录。

  • len列,为牙齿长度
  • supp列,为注射龙虎大战坐庄方法
  • dose列,为注射剂量

查看数据集,打印前10行


> head(ToothGrowth,10)
    len supp dose
1   4.2   VC  0.5
2  11.5   VC  0.5
3   7.3   VC  0.5
4   5.8   VC  0.5
5   6.4   VC  0.5
6  10.0   VC  0.5
7  11.2   VC  0.5
8  11.2   VC  0.5
9   5.2   VC  0.5
10  7.0   VC  0.5

F检验对于数据的正态性非常敏感,龙虎大战坐庄龙虎大战坐庄我 们 需要先对选定数据集进行进行正态分布检验。使用Shapiro-Will作为正态分布检验的龙虎大战坐庄方法 ,原假设H0:样本符合正态分布。


# 按不同的处理龙虎大战坐庄方法
,进行分组
> len_VC<-ToothGrowth$len[which(ToothGrowth$supp=='VC')]
> len_OJ<-ToothGrowth$len[which(ToothGrowth$supp=='OJ')]

# 正态分布检验
> shapiro.test(len_VC)

	Shapiro-Wilk normality test

data:  len_VC
W = 0.96567, p-value = 0.4284

# 正态分布检验
> shapiro.test(len_OJ)

	Shapiro-Wilk normality test

data:  len_OJ
W = 0.91784, p-value = 0.02359

两个样本的W统计量都接近1,且p-value都大于0.05,不能拒绝原假设,两组样本数据为正态分布。

查看数据的相关性。


> coplot(len ~ dose | supp, data = ToothGrowth, panel = panel.smooth,
       xlab = "ToothGrowth data: length vs dose, given type of supplement")

3. F检验实现

3.1 随机数进行F检验
龙虎大战坐庄龙虎大战坐庄我 们 先用一种随机数,来做一下F检验。以正态分布生成2组数据,数量,均值,方差都不同,进行F检验。


# 生成随机数
> set.seed(1)
> x <- rnorm(50, mean = 0, sd = 2)
> y <- rnorm(30, mean = 1, sd = 1)

# 进行F检验
> var.test(x, y)

	F test to compare two variances

data:  x and y
F = 2.6522, num df = 49, denom df = 29, p-value
= 0.006232
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 1.332510 4.989832
sample estimates:
ratio of variances 
          2.652168 

指标解释:

  • H0:原假设2组样本的方差,无显著差异
  • F统计量:2.6522
  • num df,分子自由度,50-1=49
  • denom df,分每自由度,30-1=29
  • p-value值:0.006232
  • 95 percent confidence interval:95%的置信区间
  • ratio of variances:方差比率2.652168

结果解读,以0.05为显著性水平,F = 2.6522大于临界值1.81(查表),F值显著,拒绝原假设。以0.05为显著性水平,p-value=0.006232小于0.05,拒绝原假设,两样本方差有显著性差异。这个结果与龙虎大战坐庄龙虎大战坐庄我 们 构造的数据是一致的,样本的方差就是不同的。

3.2 ToothGrowth进行F检验
使用ToothGrowth数据集进行F检验,原假设HO,用VC和OJ两种龙虎大战坐庄方法 按3种剂量进行注射,对于60只豚鼠的牙齿生长速度的方差,没有显著性差异。


> var.test(len_VC,len_OJ)

	F test to compare two variances

data:  len_VC and len_OJ
F = 1.5659, num df = 29, denom df = 29, p-value
= 0.2331
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.745331 3.290028
sample estimates:
ratio of variances 
          1.565937 

结果解读,以0.05为显著性水平,F=1.5659小于临界值1.90(查表),F值不显著,不能拒绝原假设。以0.05为显著性水平,p-value=0.2331大于0.05,不能拒绝原假设,所以两种龙虎大战坐庄方法 的3种剂量实验的方差,没有显著性的差异。

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

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


# 手动计算T值
> Xn<-length(len_VC)
> Yn<-length(len_OJ)
> Xm<-mean(len_VC)
> Ym<-mean(len_OJ)

# 计算两组样本的偏方差
> fx<-sum((len_VC-Xm)^2)/(Xn-1)
> fy<-sum((len_OJ-Ym)^2)/(Yn-1)

# 计算F值
> fx/fy
[1] 1.565937

# 手动计算P值,双边检验
> p_value<-pf(f_stat,Yn-1,Xn-1)
> p_value<-2*min(p_value, 1 - p_value);p_value
[1] 0.2331433

用F检验测试样本数据的偶然误差,对数据集进行方差齐性检验,从而判断数据是否有显著性差异,为方差分析提供了基本的判别龙虎大战坐庄方法 ,对于研究数据的波动性是非常有用的。

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

打赏作者

用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

打赏作者