Blog Archives

用R语言解读牛顿冷却定律

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

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

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

环球彩票关于 作者:

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

转载请注明出处:
http://xiyueit.com/r-newton-cooling

前言

牛顿冷却定律是经典的热力学模型,描述了温度与时间之间的指数衰减函数。环球彩票环球彩票我 们 通过对基础科学的研究,可以极大地开阔思路,并形成跨学科的解决方案,引用新的思考维度,解决现有场景的难题。站在物理学的制高点,推进数据科学的发展。

目录

  1. 牛顿冷却定律
  2. 数据集
  3. 程序实现

1. 牛顿冷却定律

牛顿冷却定律,是一个热力学的基本定律,指物体所损失的热的速率与物体和其周围环境间的温度差是成比例的。当物体表面与周围存在温度差时,单位时间从单位面积散失的热量与温度差成正比,比例系数称为热传递系数。牛顿冷却定律,是由英国物理学家艾萨克·牛顿爵士所提出的,是一个经验发现的规律。

当物体处于与周围环境不同的温度时,它将逐渐冷却或加热直至温度相等。环球彩票环球彩票我 们 每个人都在环球彩票生活中遇到过,比如煮水沏茶的过程,把凉水烧热到100度的高温,然后沏茶,再等一会儿,到达可以喝的温度。水从100度的高温,自然到达可以喝的温度,就是冷却的过程,即物体温度与周围环境温度热交换的过程。

公式:

公式解释:

  • T: 温度
  • t: 时间
  • T(t):t时刻的温度
  • Ta:代表室温,T(t)-Ta就是当前温度与室温之间的温差。当前温度高于室温时,T(t)-Ta > 0。
  • T0:初始温度
  • k:传热系数,常数,表示室温与降温速率之间的比例关系
  • dT/dt:温度和时间的微分形式,表示物体的温度随时间下降的速度

把上面微分方程进行数学变型,产生可以下面等式。

经过数学变型后,就容易进行计算了。由于k是一个常数,所以环球彩票环球彩票我 们 需要先计算出k值,再用k值计算出每一时刻的模型值,最后用模型值和测试值进行对比。

进行k值的公式变型:

2. 数据集

用数据说话,环球彩票环球彩票我 们 通过一个实验来验证牛顿冷去定律,把牛顿冷却定律计算的数值与真实数据的测量值进行拟合。实验数据来源:Newton’s Law of Cooling An Experimental Investigation

实验开始,准备3个烧杯加热至沸腾,让3个烧杯中的水自然冷却,记录60分钟每分钟3个烧杯中水的温度。3个烧杯水中水量不同,第一个烧杯放入100毫升水,第二个烧杯放入300毫升水,第三个烧杯放入800毫升水。每个烧杯都有自己的温度计,温度计在测量之间保持在烧杯中,确保不会出现温度滞后。每分钟测量每个烧杯中水的温度,环境温度为23°C。

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

  • 环球彩票Win 10 64bit
  • R: 3.6.1 x86_64-w64-mingw32/x64 b4bit

记录3个烧杯温度的数据集,共4列,保存文件为Temperature.csv。

  • time_min:分钟,冷却时间,从0分钟到60分钟
  • temp_100ml:温度(摄氏度),100ml水随时间变化的温度
  • temp_300ml:温度(摄氏度),300ml水随时间变化的温度
  • temp_800ml:温度(摄氏度),800ml水随时间变化的温度

> dat
   time_min temp_100ml temp_300ml temp_800ml
1         0        100        100        100
2         1         95         95         96
3         2         82         91         95
4         3         79         87         92
5         4         74         84         90
6         5         70         81         88
7         6         67         78         85
8         7         65         76         83
9         8         61         73         80
10        9         59         71         78
11       10         57         70         76
12       11         56         68         75
13       12         54         66         74
14       13         52         64         73
15       14         51         63         71
16       15         50         61         70
17       16         49         60         68
18       17         48         58         66
19       18         47         58         66
20       19         45         56         65
21       20         45         55         63
22       21         44         55         62
23       22         43         54         61
24       23         42         53         60
25       24         42         52         60
26       25         41         51         59
27       26         41         50         58
28       27         40         49         56
29       28         39         48         56
30       29         38         48         55
31       30         38         47         54
32       31         38         46         53
33       32         38         46         52
34       33         37         45         52
35       34         36         45         51
36       35         36         45         50
37       40         34         42         47
38       45         33         40         45
39       50         31         38         43
40       55         30         37         41
41       60         29         36         40

观察数据集,用散点图画出的测试值的3种情况。


# 加载数据集
> dat<-read.csv("Temperature.csv")

# 以散点图观察数据
> plot(dat$time_min,dat$temp_100ml)
> points(dat$time_min,dat$temp_300ml,col="red")
> points(dat$time_min,dat$temp_800ml,col="blue")


图解释:

  • x轴为时间
  • y轴为维度
  • 黑色散点,100毫升水的温度变化
  • 红色散点,300毫升水的温度变化
  • 蓝色散点,800毫升水的温度变化

初步观察,3个烧杯水开始温度重合,都是100度。100毫升水随时间变化快,800毫升水随时间变化慢。

3. 程序实现

接下来,分别计算出3种情况下的水温冷却的理论值,通过牛顿冷却定律,并验证理论值与观测值之差的差别,是否能通过理论值来描述测量值。

3.1 用牛顿冷却定律计算100毫升水的温度

计算步骤:

  1. 计算1分钟内的k值
  2. 计算每个统计时间的k值,并求平均k值
  3. 根据平均k值,计算每个统计时间温度
  4. 对比理论值和测量值,进行可视化
  5. 计算理论值和测量值的均方根误差
  6. 进行理论值和测量值的统计检验
  7. 得出结论

第一步,计算1分钟内的k值。


# 室温 23度
> Ta <- 23

# 根据公式计算k值, k = ln((T(1)-Ta)/(T(0)-Ta))/t1
> k = log((dat$temp_100ml[2]-Ta)/(dat$temp_100ml[1]-Ta)) / dat$time_min[2]
> k
[1] -0.0671393

当t=1时(1分钟),传热系数k值为-0.0671393。

第二步,计算每个统计时间的k值,并求平均k值


# 每个统计时候的k值
> kv = log((dat$temp_100ml[-1]-Ta)/(dat$temp_100ml[1]-Ta)) / dat$time_min[-1]
> kv
 [1] -0.06713930 -0.13313399 -0.10615124 -0.10299495 -0.09873156 -0.09326930 -0.08659083
 [8] -0.08827741 -0.08447628 -0.08174449 -0.07702708 -0.07581818 -0.07511612 -0.07225721
[15] -0.06986457 -0.06785681 -0.06617233 -0.06476398 -0.06593489 -0.06263815 -0.06187062
[22] -0.06127605 -0.06084202 -0.05830694 -0.05813735 -0.05590129 -0.05594785 -0.05611488
[29] -0.05640535 -0.05452517 -0.05276630 -0.05111735 -0.05165903 -0.05231930 -0.05082446
[36] -0.04864775 -0.04536045 -0.04528728 -0.04359810 -0.04253410

# k平均值
> k<-mean(kv);k
[1] -0.06758501

第三步,根据平均k值,计算每个统计时间温度


# 计算理论温度
> Tt<- Ta + (dat$temp_100ml[1] - Ta)*exp(k*dat$time_min[-1])
> Tt
 [1] 94.96792 90.26469 85.86882 81.76024 77.92015 74.33103 70.97645 67.84111 64.91067 62.17173
[11] 59.61179 57.21915 54.98287 52.89273 50.93919 49.11331 47.40676 45.81174 44.32095 42.92759
[21] 41.62529 40.40809 39.27044 38.20714 37.21333 36.28446 35.41630 34.60487 33.84648 33.13764
[31] 32.47513 31.85591 31.27716 30.73624 30.23066 28.15726 26.67841 25.62362 24.87129 24.33470

第四步,对比理论值和测量值,进行可视化


# 把100ml的理论值和测量值合并数据框
> df100<-data.frame(dat[,1:2],Tt=c(100,Tt))

# 进行可视化,测试值为黑色点,理论值为红色线
> plot(df100$time_min,df100$temp_100ml,ylim=c(0,100))
> lines(df100$time_min,df100$Tt,col="red")

第五步,计算理论值和测量值的均方根误差


> sqrt(sum((df100$Tt - df100$temp_100ml)^2)/nrow(df100))
[1] 4.800323

第六步,进行理论值和测量值的统计检验。

对观测值和理论分别进行正常分步的检验,判断是否属于正态分布。


# 观测值检验
> shapiro.test(df100[,2])

	Shapiro-Wilk normality test

data:  df100[, 2]
W = 0.88026, p-value = 0.0004557

# 理论值检验
> shapiro.test(df100[,3])

	Shapiro-Wilk normality test

data:  df100[, 3]
W = 0.90219, p-value = 0.001925

通过shapiro-wilk检验,以0.05为显著性水平,p-value都小于0.05,所以拒绝原假设,观测值和预测值都不符合正态分布,也就不能使用T检验和F检验,对数据集是分析了。

卡方检验,用于检验2个独立样本的观测值与期望值的差距。


# 卡方检验
> x <- matrix(c(df100[,2],df100[,3]), ncol = 2)
> chisq.test(x)

	Pearson's Chi-squared test

data:  x
X-squared = 10.266, df = 40, p-value = 1

以0.05为显著性水平,p-value=1都大于0.05,所以不能拒绝原假设,观测值和理论值没有明显的差异性,用理论值来近似观测值的环球彩票方法 是合理的。

第七步,得出结论,使用牛顿冷却定律进行理论值计算与实际观测值没有明显的差异性,用理论值来近似观测值的环球彩票方法 是合理的,根绝牛顿冷却定律得出的理论值可以作为有效的热传递的计算环球彩票方法 。

再分别计算300毫升水和800毫升水的理论值和观测值之间的关系。

3.2 用牛顿冷却定律计算300毫升水的温度
300毫升水,计算理论值。


# 300毫升水
> kv = log((dat$temp_300ml[-1]-Ta)/(dat$temp_300ml[1]-Ta)) / dat$time_min[-1]

# k均值
> k<-mean(kv);k
[1] -0.0447057

# 计算每个统计时间的温度
> Tt<- Ta + (dat$temp_300ml[1] - Ta)*exp(k*dat$time_min[-1]);Tt
 [1] 96.63347 93.41413 90.33555 87.39156 84.57629 81.88411 79.30963 76.84771 74.49342 72.24207
[11] 70.08915 68.03036 66.06159 64.17888 62.37850 60.65682 59.01043 57.43601 55.93043 54.49067
[21] 53.11386 51.79725 50.53820 49.33420 48.18284 47.08182 46.02894 45.02208 44.05925 43.13852
[31] 42.25804 41.41606 40.61089 39.84092 39.10461 35.87873 33.29902 31.23605 29.58630 28.26701

# 合并理论值和观测值
> df300<-data.frame(dat[,c('time_min',"temp_300ml")],Tt=c(100,Tt))

# 均方根误差
> sqrt(sum((df300$Tt - df300$temp_300ml)^2)/nrow(df300))
[1] 3.711398

# 可视化
> plot(df300$time_min,df300$temp_300ml,ylim=c(0,100))
> lines(df300$time_min,df300$Tt,col="red")

3.3 用牛顿冷却定律计算800毫升水的温度
800毫升水,计算理论值


> kv = log((dat$temp_800ml[-1]-Ta)/(dat$temp_800ml[1]-Ta)) / dat$time_min[-1]

# k均值
> k<-mean(kv);k
[1] -0.03266545

# 计算每个统计时间的温度
> Tt<- Ta + (dat$temp_800ml[1] - Ta)*exp(k*dat$time_min[-1]);Tt
 [1] 97.52540 95.13032 92.81222 90.56862 88.39712 86.29540 84.26123 82.29244 80.38692 78.54263
[11] 76.75762 75.02998 73.35785 71.73947 70.17309 68.65706 67.18974 65.76959 64.39507 63.06473
[21] 61.77714 60.53093 59.32478 58.15738 57.02750 55.93394 54.87552 53.85111 52.85963 51.90001
[31] 50.97123 50.07230 49.20226 48.36018 47.54516 43.84653 40.70523 38.03729 35.77137 33.84689

# 合并理论值和观测值
> df800<-data.frame(dat[,c('time_min',"temp_800ml")],Tt=c(100,Tt))

# 均方根误差
> sqrt(sum((df800$Tt - df800$temp_800ml)^2)/nrow(df800))
[1] 2.237821

# 可视化
plot(df800$time_min,df800$temp_800ml,ylim=c(0,100))
lines(df800$time_min,df800$Tt,col="red")


3.4 总结
最后,把3种情况的评价指标整理到一起,形成一个表格。

水量 k平均值 均方根误差
100ml -0.06758501 4.800323
300ml -0.0447057 3.711398
800ml -0.03266545 2.237821

800ml的烧瓶的水,均方根误差最小,理论值与观测值拟合最好;同时,传热系数k值最小,说明水量越多传热越慢。100ml的烧杯的水,均方根误差最大,理论值与观测值拟合最差;同时,传热系数k值最大,说明水量越少传热越快。

牛顿冷却定律是一种热力学模型,通过温度与时间之间的函数关系,构建出了一个指数衰减的过程,这种衰减过程不仅能表达热传导的关系,还被应用到了如新闻热度排名领域。对基础科学的公理和定理的研究,可以极大地开阔研究数据的思路,形成跨学科的解决方案,下篇文章环球彩票我 将介绍牛顿冷却定律在基于用户投票的排名算法场景中的应用。

转载请注明出处:
http://xiyueit.com/r-newton-cooling

打赏作者

用R语言解读统计检验-卡方检验

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

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

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

环球彩票关于 作者:

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

转载请注明出处:
http://xiyueit.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://xiyueit.com/r-test-x2

打赏作者

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

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

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

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

环球彩票关于 作者:

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

转载请注明出处:
http://xiyueit.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://xiyueit.com/r-test-f

打赏作者

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

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

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

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

环球彩票关于 作者:

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

转载请注明出处:
http://xiyueit.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://xiyueit.com/r-test-t

打赏作者

R语言中的矩阵计算

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

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

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

环球彩票关于 作者:

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

转载请注明出处:
http://xiyueit.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://xiyueit.com/r-matrix

打赏作者