R語(yǔ)言在概率與統(tǒng)計(jì)及線性代數(shù)領(lǐng)域中的應(yīng)用

1、. R語(yǔ)言在概率與統(tǒng)計(jì)及線性代數(shù)領(lǐng)域中的應(yīng)用 一、概率與統(tǒng)計(jì)：概率與統(tǒng)計(jì)在科學(xué)領(lǐng)域、工程領(lǐng)域與社會(huì)研究領(lǐng)域都有相當(dāng)強(qiáng)的實(shí)用性。最近，機(jī)率與統(tǒng)計(jì)在人工智能、自然語(yǔ)言翻譯、語(yǔ)音識(shí)別、機(jī)器人控制等領(lǐng)域，都有相當(dāng)深入的進(jìn)展，但是這些進(jìn)展所使用到的數(shù)學(xué)，都需要相當(dāng)深的機(jī)率統(tǒng)計(jì)基礎(chǔ)。 在實(shí)務(wù)研究上，要重視“程序語(yǔ)言實(shí)作”的方式，因此采用R 這個(gè)專為機(jī)率統(tǒng)計(jì)而設(shè)計(jì)的語(yǔ)言來(lái)做研究，透過(guò) R 語(yǔ)言，我們可以較為輕松的實(shí)作出許多用傳統(tǒng)語(yǔ)言難以實(shí)作的程序范例，這些范例可以引導(dǎo)我們深入理解機(jī)率統(tǒng)計(jì)的理論，并且達(dá)到“以理論指導(dǎo)實(shí)務(wù)，以實(shí)務(wù)印證理論”的功能。R 語(yǔ)言當(dāng)中默認(rèn)就包含了各式各樣的機(jī)率模型，以及各種統(tǒng)計(jì)工具。舉

2、例而言，當(dāng)看到常態(tài)分配的機(jī)率模型時(shí)，我們可以直接使用下列指令來(lái)畫出常態(tài)分布，并且用程序產(chǎn)生符合常態(tài)分布模型的樣本，以便進(jìn)行某種交互式的學(xué)習(xí)，用實(shí)驗(yàn)體會(huì)常態(tài)分布的意義。> 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)")>以上只是一個(gè)繪圖的范例而

3、已，R 語(yǔ)言能輕易做出各式各樣的機(jī)率統(tǒng)計(jì)程序，像是抽樣、敘述報(bào)告、檢定、推論、回歸分析、變異數(shù)分析 (ANOVA)、質(zhì)量管理、時(shí)間數(shù)列等等。范例：(1)模擬擲骰子：> sample(1:6, 1)1 2> sample(1:6, 5)1 4 1 3 2 5> sample(1:6, 10)錯(cuò)誤在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)統(tǒng)計(jì)抽樣結(jié)果：> x <- sample(6,6000,replace=T)> for (i in 1:6) print(sum(x=i)1 10361 9781 10371 10161 9861 947繪制統(tǒng)計(jì)圖：> 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.線性代數(shù)：使用 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錯(cuò)誤于y %*% x : 非整合參數(shù)> 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)人工智能中的機(jī)率統(tǒng)計(jì)范例 - 牙醫(yī)判斷牙痛問(wèn)題問(wèn)題描述：當(dāng)病人來(lái)看牙醫(yī)時(shí)，該病人可能有蛀牙或沒蛀牙，也可能有牙痛或沒有牙痛，而牙醫(yī)可能會(huì)找到牙痛的原因或找不到。因此有下列三個(gè)隨機(jī)變量：X:(蛀) 蛀牙與否 (Cavity)Y:(痛) 牙痛與否 (Toothache)Z:(找) 是否找到痛的牙 (C

15、atch)假如這個(gè)問(wèn)題個(gè)統(tǒng)計(jì)機(jī)率都已經(jīng)知道了，如下表所示：牙痛 (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如果將這個(gè)表格寫成一整排，那么將會(huì)以下列機(jī)率表格顯示。蛀 X痛 Y找 ZP(X,Y,Z)0000.5760010.1440100.0640110.0161000.0081010.0721100.0121110.108問(wèn)題1：請(qǐng)計(jì)算 P(沒痛) = P(Y=0) = ?2：請(qǐng)計(jì)算 P(找到 | 牙痛) = ?3：請(qǐng)問(wèn)這是一個(gè)合理

16、的機(jī)率分布嗎？4：請(qǐng)計(jì)算 P(Z=1 | X=1) = ?5：請(qǐng)計(jì)算 P(Z=1, Y=1) = ?6：請(qǐng)計(jì)算 P(蛀 | 找到), P(蛀), P(找到), P(找到 | 蛀) ，然后驗(yàn)證下列貝氏定理是否成立。(1) P(找到|蛀)=P(蛀|找到)P(找到) P(蛀)解答：由于 R 的 數(shù)組是用以行為主的順序 (Column Major Order)，因此沒辦法直接與上列的表格對(duì)起來(lái)。所以我們必須先將真值表改為以行為主的方式，改寫后如下表所示：蛀 X痛 Y找 ZP(X,Y,Z)0000.5761000.0080100.0641100.0120010.1441010.0720110.01611

17、10.108由于 R 的數(shù)組是從 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計(jì)算過(guò)程：> 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：請(qǐng)問(wèn)這是一個(gè)合理的機(jī)率分布嗎？ (是的，因?yàn)榭偤蜑?1，而且每個(gè)機(jī)率直都介于 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：請(qǐng)計(jì)算 P(Z=1 | X=1) = 0.9> sum(p2,2)/sum(p2,)1 0.9解答5：請(qǐng)計(jì)算 P(Z=1, Y=1) = 0.124> p,2,2 沒蛀 蛀 0.016 0.108 > sum(p,2,2)1 0.124解答6：請(qǐng)計(jì)算 P(蛀 | 找到), P(蛀), P(找到), P(找到 | 蛀) ，然后驗(yàn)證下列貝氏定理是否成立。P(蛀 | 找到) = p(找到|蛀) * p(蛀)/p(找到)說(shuō)明：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.線性代數(shù)：使用 R 軟件計(jì)算最小二乘法： > 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.線性代數(shù)：使用 R 軟件計(jì)算指數(shù)矩陣矩陣指數(shù)（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