R语言在概率与统计及线性代数领域中的应用

1、. R语言在概率与统计及线性代数领域中的应用 一、概率与统计：概率与统计在科学领域、工程领域与社会研究领域都有相当强的实用性。最近，机率与统计在人工智能、自然语言翻译、语音识别、机器人控制等领域，都有相当深入的进展，但是这些进展所使用到的数学，都需要相当深的机率统计基础。 在实务研究上，要重视“程序语言实作”的方式，因此采用R 这个专为机率统计而设计的语言来做研究，透过 R 语言，我们可以较为轻松的实作出许多用传统语言难以实作的程序范例，这些范例可以引导我们深入理解机率统计的理论，并且达到“以理论指导实务，以实务印证理论”的功能。R 语言当中默认就包含了各式各样的机率模型，以及各种统计工具。举

2、例而言，当看到常态分配的机率模型时，我们可以直接使用下列指令来画出常态分布，并且用程序产生符合常态分布模型的样本，以便进行某种交互式的学习，用实验体会常态分布的意义。> dnorm(0)1 0.3989423> dnorm(0.5)1 0.3520653> dnorm(2.5)1 0.0175283> curve(dnorm(x), from = -3.5, to = 3.5, ylab="pdf", main="N(0,1) probability density function (pdf)")>以上只是一个绘图的范例而

3、已，R 语言能轻易做出各式各样的机率统计程序，像是抽样、叙述报告、检定、推论、回归分析、变异数分析 (ANOVA)、质量管理、时间数列等等。范例：(1)模拟掷骰子：> sample(1:6, 1)1 2> sample(1:6, 5)1 4 1 3 2 5> sample(1:6, 10)错误在sample(1:6, 10) : cannot take a sample larger than the population when 'replace = FALSE'> sample(1:6, 10, replace=T) 1 4 3 2 4 1 6 4

4、 2 5 1> > sample(1:6, 10, replace=T) 1 4 3 2 4 1 6 4 2 5 1> x <- sample(1:6, 10, replace=T)> x 1 1 2 5 2 4 1 4 4 6 5> x=2> x1 2> x <- sample(1:6, 10, replace=T)> x 1 2 4 2 2 5 4 4 5 3 2> sum(x=2)1 2> x <- sample(1:6, 10, replace=T)> x 1 4 1 6 2 1 3 5 5 3 6&g

5、t; x=2 1 FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE> sum(x=2)1 1> sum(x=6)1 2>(2)统计抽样结果：> x <- sample(6,6000,replace=T)> for (i in 1:6) print(sum(x=i)1 10361 9781 10371 10161 9861 947绘制统计图：> x = sample(1:6, 10000, replace=T)> hist(x)>(3)用 R 软件仿真扑克牌：> num

6、=1:13> type=c("a","b","c","d")> num*type> type=c("space", "heart", "diamond", "club")> interaction(num, type) 1 1.space 2.heart 3.diamond 4.club 5.space 6.heart 7 7.diamond 8.club 9.space 10.heart 11.diamond

7、12.club 13 13.space 52 Levels: 1.club 2.club 3.club 4.club 5.club 6.club 7.club . 13.spaceWarning message:In ans * length(l) + if1 : longer object length is not a multiple of shorter object length> expand.grid(num, type) Var1 Var21 1 space2 2 space3 3 space4 4 space5 5 space6 6 space7 7 space8 8

8、space9 9 space10 10 space11 11 space12 12 space13 13 space14 1 heart15 2 heart16 3 heart17 4 heart18 5 heart19 6 heart20 7 heart21 8 heart22 9 heart23 10 heart24 11 heart25 12 heart26 13 heart27 1 diamond28 2 diamond29 3 diamond30 4 diamond31 5 diamond32 6 diamond33 7 diamond34 8 diamond35 9 diamond

9、36 10 diamond37 11 diamond38 12 diamond39 13 diamond40 1 club41 2 club42 3 club43 4 club44 5 club45 6 club46 7 club47 8 club48 9 club49 10 club50 11 club51 12 club52 13 club> poker = expand.grid(num, type)二、1.线性代数：使用 R 软件求解特征值与特征向量 > matrix(rep(3,8), nrow=4) ,1 ,21, 3 32, 3 33, 3 34, 3 3> m

10、atrix(3, nrow=4, ncol=3) ,1 ,2 ,31, 3 3 32, 3 3 33, 3 3 34, 3 3 3> w=matrix(0, nrow=3, ncol=4)> w ,1 ,2 ,3 ,41, 0 0 0 02, 0 0 0 03, 0 0 0 0> w2,3=5> w ,1 ,2 ,3 ,41, 0 0 0 02, 0 0 5 03, 0 0 0 0> x=matrix(1:12, nrow=3, ncol=4)> x ,1 ,2 ,3 ,41, 1 4 7 102, 2 5 8 113, 3 6 9 12> dim(x

11、)1 3 4> y = 1:12> y %*% x错误于y %*% x : 非整合参数> y = 1:3> y %*% x ,1 ,2 ,3 ,41, 14 32 50 68> t(x) ,1 ,2 ,31, 1 2 32, 4 5 63, 7 8 94, 10 11 12> diag(c(3,7,2) ,1 ,2 ,31, 3 0 02, 0 7 03, 0 0 2> diag(x)1 1 5 9> z = matrix(1:9, nrow=3, ncol=3)> det(z)1 0> z ,1 ,2 ,31, 1 4 72, 2

12、5 83, 3 6 9> z2,1=10> z ,1 ,2 ,31, 1 4 72, 10 5 83, 3 6 9> det(z)1 48> eigen(z)$values1 17.936765+0.000000i -1.468383+0.721055i -1.468383-0.721055i$vectors ,1 ,2 ,31, 0.4100596+0i 0.2162169-0.1118877i 0.2162169+0.1118877i2, 0.6875542+0i -0.8668747+0.0000000i -0.8668747+0.0000000i3, 0.599

13、2665+0i 0.4306386+0.0617265i 0.4306386-0.0617265i> x ,1 ,2 ,3 ,41, 1 4 7 102, 2 5 8 113, 3 6 9 12> rownames(x) = paste("r", 1:3, sep="_")> x ,1 ,2 ,3 ,4r_1 1 4 7 10r_2 2 5 8 11r_3 3 6 9 12> colnames(x) = paste("c", 1:4, sep="")> x c1 c2 c3 c4r_1

14、 1 4 7 10r_2 2 5 8 11r_3 3 6 9 12> hist(x)> 范例：(1)R 软件与排列组合> prod(5:1)1 120> prod(5:3)1 60> choose(5,3)1 10> choose(52, 4)1 270725(2)人工智能中的机率统计范例 - 牙医判断牙痛问题问题描述：当病人来看牙医时，该病人可能有蛀牙或没蛀牙，也可能有牙痛或没有牙痛，而牙医可能会找到牙痛的原因或找不到。因此有下列三个随机变量：X:(蛀) 蛀牙与否 (Cavity)Y:(痛) 牙痛与否 (Toothache)Z:(找) 是否找到痛的牙 (C

15、atch)假如这个问题个统计机率都已经知道了，如下表所示：牙痛 (Y=1)不牙痛 (Y=0)找到 (Z=1)找不到 (Z=0)找到 (Z=1)找不到 (Z=0)蛀牙(X=1)0.1080.0120.0720.008沒蛀牙 (X=0)0.0160.0640.1440.576如果将这个表格写成一整排，那么将会以下列机率表格显示。蛀 X痛 Y找 ZP(X,Y,Z)0000.5760010.1440100.0640110.0161000.0081010.0721100.0121110.108问题1：请计算 P(没痛) = P(Y=0) = ?2：请计算 P(找到 | 牙痛) = ?3：请问这是一个合理

16、的机率分布吗？4：请计算 P(Z=1 | X=1) = ?5：请计算 P(Z=1, Y=1) = ?6：请计算 P(蛀 | 找到), P(蛀), P(找到), P(找到 | 蛀) ，然后验证下列贝氏定理是否成立。(1) P(找到|蛀)=P(蛀|找到)P(找到) P(蛀)解答：由于 R 的 数组是用以行为主的顺序 (Column Major Order)，因此没办法直接与上列的表格对起来。所以我们必须先将真值表改为以行为主的方式，改写后如下表所示：蛀 X痛 Y找 ZP(X,Y,Z)0000.5761000.0080100.0641100.0120010.1441010.0720110.01611

17、10.108由于 R 的数组是从 1 开始算的，因此还必须将上表修改如下：蛀 X痛 Y找 ZP(X,Y,Z)1110.5762110.0081210.0642210.0121120.1442120.0721220.0162220.108> p <- array(c(0.576, 0.008, 0.064, 0.012, 0.144, 0.072, 0.016, 0.108),c(2,2,2)> p, , 1 ,1 ,21, 0.576 0.0642, 0.008 0.012, , 2 ,1 ,21, 0.144 0.0162, 0.072 0.108> p1,1,11

18、0.576> p2,1,11 0.008> p1,2,11 0.064> p2,2,11 0.012> p1,1,21 0.144> p2,1,21 0.072> p1,2,21 0.016> p2,2,21 0.108> dimnames(p)1 = c("沒蛀", "蛀")> dimnames(p)2 = c("沒痛", "痛")> dimnames(p)3 = c("沒找", "找")> p, , 沒找

19、沒痛 痛沒蛀 0.576 0.064蛀 0.008 0.012, , 找 沒痛 痛沒蛀 0.144 0.016蛀 0.072 0.108解答1：P(没痛) = 0.8计算过程：> p,"沒痛", 沒找 找沒蛀 0.576 0.144蛀 0.008 0.072> sum(p,"沒痛",)1 0.8解答2：P(找到 | 牙痛) = 0.62> p,"找" 沒痛 痛沒蛀 0.144 0.016蛀 0.072 0.108> sum(p,"找")1 0.34> sum(p,"痛&quo

20、t;,"找")1 0.124> sum(p,"痛","找")/sum(sum(p,"痛",)1 0.62解答3：请问这是一个合理的机率分布吗？ (是的，因为总和为 1，而且每个机率直都介于 0 到1之间)> sum(p)1 1> 0<=p, , 沒找 沒痛 痛沒蛀 TRUE TRUE蛀 TRUE TRUE, , 找 沒痛 痛沒蛀 TRUE TRUE蛀 TRUE TRUE> p<=1, , 沒找 沒痛 痛沒蛀 TRUE TRUE蛀 TRUE TRUE, , 找 沒痛 痛沒蛀 TRU

21、E TRUE蛀 TRUE TRUE解答4：请计算 P(Z=1 | X=1) = 0.9> sum(p2,2)/sum(p2,)1 0.9解答5：请计算 P(Z=1, Y=1) = 0.124> p,2,2 沒蛀 蛀 0.016 0.108 > sum(p,2,2)1 0.124解答6：请计算 P(蛀 | 找到), P(蛀), P(找到), P(找到 | 蛀) ，然后验证下列贝氏定理是否成立。P(蛀 | 找到) = p(找到|蛀) * p(蛀)/p(找到)说明：P(蛀 | 找到) = 0.5294118, P(蛀)=0.2, P(找到)=0.34, P(找到 | 蛀)=0.9P

22、(蛀 | 找到) = 0.5294118 = 0.9 * 0.2/0.34 = = p(找到|蛀) * p(蛀)/p(找到)> pab = sum(p"蛀","找")/sum(p,"找") ; pab = P(蛀 | 找到)> pba = sum(p"蛀","找")/sum(p"蛀",) ; pba = P(找到 | 蛀)> pa = sum(p"蛀",) ; pa = P(蛀)> pb = sum(p,"找")

23、 ; pb = P(找到)> pab1 0.5294118> pba1 0.9> pa1 0.2> pb1 0.34> pba*pa/pb1 0.5294118> pab-pba*pa/pb1 0>所以：p(蛀|找) = sum(p"蛀","找")/sum(p,"找") = pab = pba * pa / pb = p(找|蛀) * p(蛀)/p(找) = sum(p"蛀","找")/sum(p,"蛀")* sum(p,"

24、蛀")/ sum(p"找",)2.线性代数：使用 R 软件计算最小二乘法： > x<-c(0.0, 0.2, 0.4, 0.6, 0.8)> y<-c(0.9, 1.9, 2.8, 3.3, 4.2)> l <- lsfit(x, y)> l$coefficientsIntercept X1.02 4.00$residuals1 -0.12 0.08 0.18 -0.12 -0.02$intercept1 TRUE$qr$qt1 -5.85849810 2.52982213 0.23749843 -0.02946714 0

25、.10356728$qrIntercept X1, -2.2360680 -0.89442722, 0.4472136 0.63245553, 0.4472136 -0.19543954, 0.4472136 -0.51166735, 0.4472136 -0.8278950$qraux1 1.447214 1.120788$rank1 2$pivot1 1 2$tol1 1e-07attr(,"class")1 "qr"最小二乘法也可用 QR 分解完成> X<-matrix(c(rep(1,5), x), ncol=2)> X,1 ,

26、21, 1 0.02, 1 0.23, 1 0.44, 1 0.65, 1 0.8> qr(X)$qr,1 ,21, -2.2360680 -0.89442722, 0.4472136 0.63245553, 0.4472136 -0.19543954, 0.4472136 -0.51166735, 0.4472136 -0.8278950$rank1 2$qraux1 1.447214 1.120788$pivot1 1 2attr(,"class")1 "qr"3.线性代数：使用 R 软件计算指数矩阵矩阵指数（1）eA =I+A+A2 /2!+A3 /3!+.> example(expm)expm> (m1 <- Matrix(c(1,0,1,1), nc = 2)2 x 2 Matrix of class "dtrMatrix",1 ,21, 1 12, . 1expm> (e1 <- expm(m1) ; e <- exp(1)2 x 2 Matrix of class "dtrMatrix",1 ,21, 2.718282 2.7182822, . 2.718282expm> stopifnot(all.eq