R語言學(xué)習(xí)筆記三

10)求解線性方程組和逆矩陣

Solve函數(shù)求出a %*% x = b中的x向量值，即求解線性方程組，通常使用前2個(gè)參數(shù)，第一個(gè)是a，為系數(shù)矩陣  ，第二是b為常數(shù)項(xiàng)，當(dāng)b缺失時(shí)，默認(rèn)為單位矩陣。

ax=b=>ax=I=>a-1ax=a-1I=a-1 =>x=a-1

從上面推導(dǎo)可看出，不提供參數(shù)b就能求出a的逆。

求線性方程組的解

> b

[,1]

[1,]    8

[2,]    9

> a

[,1] [,2]

[1,]    1    3

[2,]    2    4

> solve(a,b)

[,1]

[1,] -2.5

[2,]  3.5

求矩陣的逆

> a

[,1] [,2]

[1,]    1    3

[2,]    2    4

> solve(a)

[,1] [,2]

[1,]   -2  1.5

[2,]    1 -0.5

11)求矩陣的特征值

λ是A的特征值等價(jià)于說線性系統(tǒng)（A– λI） v = 0 (其中I是單位矩陣)有非零解v (一個(gè)特征向量)，因此等價(jià)于說行列式:

 

函數(shù):是一個(gè)關(guān)于λ的多項(xiàng)式，稱為A的特征多項(xiàng)式。矩陣的特征值也就是其特征多項(xiàng)式的零點(diǎn)。求一個(gè)矩陣A的特征值可以通過求解方程來得到。

Eigen函數(shù)可以求特征值，其調(diào)用格式如下：

eigen(x, symmetric, only.values = FALSE)

x為需要求特征值的矩陣；symmetric是邏輯型表示是否是對(duì)稱矩陣。

對(duì)稱矩陣是一個(gè)方形矩陣，其轉(zhuǎn)置矩陣和自身相等。

對(duì)稱矩陣中的右上至左下方向元素以主對(duì)角線（左上至右下）為軸進(jìn)行對(duì)稱。

only.values 如果為TRUE，則只返回特征值，否則返回特征值和特征向量。

> a<-array(c(1:16),dim=c(4,4))

> eigen(a)

$values

[1]  3.620937e+01 -2.209373e+00  1.599839e-15  7.166935e-16

$vectors

[,1]        [,2]       [,3]       [,4]

[1,] 0.4140028  0.82289268 -0.5477226  0.1125155

[2,] 0.4688206  0.42193991  0.7302967  0.2495210

[3,] 0.5236384  0.02098714  0.1825742 -0.8365883

[4,] 0.5784562 -0.37996563 -0.3651484  0.4745519

> eigen(a,only.values=FALSE)

$values

[1]  5.3722813 -0.3722813

$vectors

[,1]       [,2]

[1,] -0.5657675 -0.9093767

[2,] -0.8245648  0.4159736

12)求矩陣的行列式值 

> a<-array(c(1:4),dim=c(2,2))

> det(x)

[1] -2

13)奇異分解

M的奇異值分解 :M=UDV',其中U'U=V'V=I。

V的列(columns)組成一套對(duì)的正交"輸入"或"分析"的基向量，是x的左奇異矩陣。U的列(columns)組成一套對(duì)的正交"輸出"的基向量，是x的右奇異矩陣。D返回一個(gè)向量，向量的元素是對(duì)角線上的元素。Svd函數(shù)完成奇異分解

> array(c(1:16),dim=c(4,4))->a

> a

[,1] [,2] [,3] [,4]

[1,]    1    5    9   13

[2,]    2    6   10   14

[3,]    3    7   11   15

[4,]    4    8   12   16

> svd(a)

$d

[1] 3.862266e+01 2.071323e+00 1.291897e-15 6.318048e-16

$u

[,1]       [,2]       [,3]       [,4]

[1,] -0.4284124 -0.7186535  0.5462756 -0.0397869

[2,] -0.4743725 -0.2738078 -0.6987120  0.4602190

[3,] -0.5203326  0.1710379 -0.2414027 -0.8010772

[4,] -0.5662928  0.6158835  0.3938391  0.3806452

$v

[,1]        [,2]       [,3]       [,4]

[1,] -0.1347221  0.82574206 -0.4654637 -0.2886928

[2,] -0.3407577  0.42881720  0.4054394  0.7318599

[3,] -0.5467933  0.03189234  0.5855124 -0.5976414

[4,] -0.7528288 -0.36503251 -0.5254881  0.1544743

14)構(gòu)造矩陣

使用matrix函數(shù)構(gòu)造矩陣，主要參數(shù)為：Data表示構(gòu)造所需數(shù)據(jù)，nrow為行數(shù),ncol為列數(shù)，byrow表示是否按行順序分配元素(默認(rèn)為FALSE)。

> matrix(c(1:10),2,5,TRUE)

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

[1,]    1    2    3    4    5

[2,]    6    7    8    9   10

> matrix(c(1:10),2,5)

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

[1,]    1    3    5    7    9

[2,]    2    4    6    8   10

15)最小二乘法擬合

最小二乘法（generalized least squares）是一種數(shù)學(xué)優(yōu)化技術(shù)，它通過最小化誤差的平方和找到一組數(shù)據(jù)的最佳函數(shù)匹配。 

設(shè)存在(x,y)這2個(gè)變量，對(duì)于一系列的x變量值，有一系列的y值與其對(duì)應(yīng)，可以找到這2個(gè)變量之間的相互關(guān)系。我們以一次函數(shù)為例，通過將這些(x,y)值標(biāo)注在直角坐標(biāo)系統(tǒng)，可以得到一條直線，讓這些點(diǎn)在這條直線附近，設(shè)直線方程為

y=kx+b 

其中：k、b 是任意實(shí)數(shù)

目標(biāo)是這些點(diǎn)到這條直線的距離的平方和最小，可運(yùn)用最小二乘法，最小二乘法擬合的過程就是回歸，這條直線就是回歸線。

Lsfit()函數(shù)完成最小二乘法擬合，其主要參數(shù)為：

X：一個(gè)矩陣的行對(duì)應(yīng)的情況和其列對(duì)應(yīng)為變量。

Y：結(jié)果，可以是一個(gè)矩陣，如果你想，以適應(yīng)多種左手側(cè)。

Wt:可選參數(shù)，加權(quán)最小二乘法的執(zhí)行權(quán)重向量。

Intercept:是否應(yīng)使用截距項(xiàng)。

Tolerance:公差將用于在矩陣分解

Yname:用于響應(yīng)變量的名稱。

我們以x=(1,2,3,4),y=(2,4,6,8)，可得到回歸線方程為Y=2x

> y<-c(2,4,6,8)

> x<-c(1,2,3,4)

> lsfit(x,y)

$coefficients

Intercept         X 

0         2 

上述結(jié)果中,Intercept項(xiàng)表示截距，x項(xiàng)表示方程的x的常數(shù)項(xiàng)。我們先假設(shè)回歸線為Y=2x+3，然后，根據(jù)回歸線構(gòu)造x和y值。

> y<-c(5,7,9,11)

> x<-c(1,2,3,4)

執(zhí)行l(wèi)sfit()函數(shù)

> lsfit(x,y)

$coefficients

Intercept         X 

3         2 

要正確得出方程的截距為3，x的常數(shù)項(xiàng)為2?，F(xiàn)實(shí)生活中，很難有如此精確的模型，我們?cè)俣鄻?gòu)造一些點(diǎn)：

> y<-c(5,7,9,11,16,20)

> x<-c(1,2,3,4,7,9)

> lsfit(x,y)

> x<-c(1,2,3,4,7,9)

> y<-c(5,7,9,11,16,20)

我們通過plot(x,y)來繪制這些點(diǎn)在直角坐標(biāo)系中的位置，這個(gè)圖也被稱為散點(diǎn)圖。

> plot(x,y)

> lsfit(x,y)

$coefficients

Intercept         X 

3.338028  1.845070 

$residuals

[1] -0.18309859 -0.02816901  0.12676056  0.28169014 -0.25352113  0.05633803

Coefficients為系數(shù)，包括截距和x的系數(shù)，residuals表示殘差，殘差分別反應(yīng)了這些點(diǎn)與直線的差異，殘差越小越好，我們將回歸線也畫上

> abline(lsfit(x,y))

可以看到擬合效果還是不錯(cuò)的，我們也可以使用lm（）函數(shù)，來建立線性模型進(jìn)行回歸分析：

畫x,y的散點(diǎn)圖： plot(x,y)

做相關(guān)回歸分析，結(jié)果存放在xy中： lm(y~x)->xy

顯示xy的相關(guān)回歸分析結(jié)果:summary(xy)

畫回歸線:>  abline( lm(y~x))

16)矩陣與向量連接

對(duì)于向量，Cbind將行轉(zhuǎn)變?yōu)榱羞M(jìn)行連接，而rbind將列轉(zhuǎn)變?yōu)樾羞M(jìn)行連接。

> x2<-c(101:105)

> x1<-c(1:10)

> cbind(x1,x2)

x1  x2

[1,]  1 101

[2,]  2 102

[3,]  3 103

[4,]  4 104

[5,]  5 105

[6,]  6 101

[7,]  7 102

[8,]  8 103

[9,]  9 104

[10,] 10 105

> rbind(x1,x2)

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]

x1    1    2    3    4    5    6    7    8    9    10

x2  101  102  103  104  105  101  102  103  104   105

對(duì)于矩陣，cbind完成橫向拼接，rbind完成縱向拼接，下面是關(guān)于矩陣的操作

> x3<-matrix(c(1:10),2,5)

> x4<-matrix(c(101:105),2,5)

> x3

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

[1,]    1    3    5    7    9

[2,]    2    4    6    8   10

> x4

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

[1,]  101  103  105  102  104

[2,]  102  104  101  103  105

> cbind(x3,x4)

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]

[1,]    1    3    5    7    9  101  103  105  102   104

[2,]    2    4    6    8   10  102  104  101  103   105

>  rbind(x3,x4)

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

[1,]    1    3    5    7    9

[2,]    2    4    6    8   10

[3,]  101  103  105  102  104

[4,]  102  104  101  103  105

As.vector將數(shù)組轉(zhuǎn)換為向量。

> x<-array(c(1:10),dim=c(2,5))

> x

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

[1,]    1    3    5    7    9

[2,]    2    4    6    8   10

> as.vector(x)

[1]  1  2  3  4  5  6  7  8  9 10

交叉因子頻率表

> y<-c(11,22,13,14,11,22,31,31,31,14)

> table(y)

y

11 13 14 22 31 

2  1  2  2  3 

cut把數(shù)值型對(duì)象分區(qū)間轉(zhuǎn)換為因子table

> cut(c(1:10),5)

[1] (0.991,2.79] (0.991,2.79] (2.79,4.6]   (2.79,4.6]   (4.6,6.4]   

[6] (4.6,6.4]    (6.4,8.21]   (6.4,8.21]   (8.21,10]    (8.21,10]   

Levels: (0.991,2.79] (2.79,4.6] (4.6,6.4] (6.4,8.21] (8.21,10]

我們來看一個(gè)綜合的例子，求出下面樣本的數(shù)字在某區(qū)間內(nèi)的分布數(shù)量，即求因子頻率。下面是美國(guó)地震臺(tái)網(wǎng)公布的全球2013年5月20日22點(diǎn)到24點(diǎn)的所有發(fā)生的地震的震級(jí)。

我們計(jì)算一下地震震級(jí)的區(qū)間頻率分布：首先，將地震震級(jí)數(shù)據(jù)放入一個(gè)向量中。

> mag<-c(1.6,0.9,2.1,2.2,2.3,1.7,1.3,1.6,4.7,1.2,0.9,4.7,0.6,5.3,1.1,4.8,4,4.2,4.6,1.3,2.1,1.5,3)

> mag

[1] 1.6 0.9 2.1 2.2 2.3 1.7 1.3 1.6 4.7 1.2 0.9 4.7 0.6 5.3 1.1 4.8 4.0 4.2

[19] 4.6 1.3 2.1 1.5 3.0

然后，使用cut函數(shù)將震級(jí)分成5個(gè)區(qū)間，并建立因子

> factor(cut(mag,5))

[1] (1.54,2.48]  (0.595,1.54] (1.54,2.48]  (1.54,2.48]  (1.54,2.48] 

[6] (1.54,2.48]  (0.595,1.54] (1.54,2.48]  (4.36,5.3]   (0.595,1.54]

[11] (0.595,1.54] (4.36,5.3]   (0.595,1.54] (4.36,5.3]   (0.595,1.54]

[16] (4.36,5.3]   (3.42,4.36]  (3.42,4.36]  (4.36,5.3]   (0.595,1.54]

[21] (1.54,2.48]  (0.595,1.54] (2.48,3.42] 

Levels: (0.595,1.54] (1.54,2.48] (2.48,3.42] (3.42,4.36] (4.36,5.3]

最后，統(tǒng)計(jì)因子頻率

factor(cut(mag,5))->magfactor

> table(magfactor)

magfactor

(0.595,1.54]  (1.54,2.48]  (2.48,3.42]  (3.42,4.36]   (4.36,5.3] 

8            7            1            2            5 

可以看出2013年5月20日22點(diǎn)到24點(diǎn)期間，全球發(fā)生的地震在(0.595,1.54]內(nèi)有8起，在(1.54,2.48]有7起等。

hist函數(shù)可用來繪制直方圖

> hist(mag,breaks=5)

 

18)list函數(shù)生成一個(gè)對(duì)象，這個(gè)對(duì)象可擁有自定義的組件，組件也可擁有自定義的屬性

> list(name="students",class="101",stdt.ages=c(22,25,20),stdt.name=c("zhangsang","lisi","wangwu"))->mystudents

> mystudents

$name

[1] "students"

$class

[1] "101"

$stdt.ages

[1] 22 25 20

$stdt.name

[1] "zhangsang" "lisi"      "wangwu"   

length返回組件的數(shù)量

> length(mystudents)

[1] 4

可以直接提取組件的內(nèi)容完成計(jì)算 

> c(mystudents$stdt.name,mystudents$stdt.ages)

[1] "zhangsang" "lisi"      "wangwu"    "22"        "25"        "20"       

19)data.frame

data.frame是一個(gè)list類型，內(nèi)部可以擁有很多組件，我們接著上例構(gòu)造一個(gè)學(xué)生集的data.frame

> data.frame(name=mystudents$stdt.name,age=mystudents$stdt.ages)->mysts

> mysts

name age

1 zhangsang  22

2      lisi  25

3    wangwu  20

attach將data.frame內(nèi)的組件復(fù)制一份后將變量放到搜索路徑上 ，我們用分離出來的變量將每個(gè)人的年齡加2歲

> age+2->mysts$age

> mysts

name age

1 zhangsang  26

2      lisi  29

3    wangwu  24

> age

[1] 24 27 22

> name

[1] zhangsang lisi      wangwu   

Levels: lisi wangwu zhangsang

使用detach將組件從搜索路徑上刪除。

> detach(mysts)

> age

錯(cuò)誤: 找不到對(duì)象'age'

> name

錯(cuò)誤: 找不到對(duì)象'name'

可使用search函數(shù)顯示當(dāng)前搜索路徑 

> attach(mysts)

> search()

[1] ".GlobalEnv"        "mysts"             "package:stats"    

[4] "package:graphics"  "package:grDevices" "package:utils"    

[7] "package:datasets"  "package:methods"   "Autoloads"        

[10] "package:base"   

使用ls()函數(shù)顯示搜索路徑上的組件

> ls(2)

[1] "age"  "name"

> ls(3)

[1] "acf"                  "acf2AR"               "add.scope"           

[4] "add1"                 "addmargins"           "aggregate"           

[7] "aggregate.data.frame" "aggregate.default"    "aggregate.ts"        

[10] "AIC"                  "alias"                "anova"     

19)read.table和scan讀取文件

read.table比scan更強(qiáng)大，在文件有文件頭的情況下，指定header=TRUE可以將文件頭做為變量名。

> read.table("h:/my_docs/eqweek.csv",header=TRUE,sep=",")->earthquake

> earthquake     DateTime.Latitude.Longitude.Depth.Magnitude.MagType.NbStations.Gap.Distance.RMS.Source.EventID.Version

1                   2013-05-20T23:57:12.000+00:00,63.45,-148.291,5.5,1.6,Ml,,,,0.8,ak,ak10720946,1.3691E+12

2               2013-05-20T23:52:59.000+00:00,61.337,-152.069,81.4,2.1,Ml,,,,1.15,ak,ak10720941,1.36909E+12

3           2013-05-20T23:49:15.100+00:00,19.99,-155.426,38.2,2.2,Md,,133,0.1,0.11,hv,hv60501711,1.3691E+12

4                2013-05-20T23:46:36.000+00:00,60.498,-142.974,4.2,2.3,Ml,,,,0.43,ak,ak10720934,1.36909E+12

..........................

然后我們提取震級(jí)數(shù)據(jù)，繪制直方圖，看一下截止2013.5.20為止的最近一周全球地震震級(jí)的分布情況，從圖中可以看出，大部分的震級(jí)都在1-2級(jí)以內(nèi)，地球還是比較安全的~

> hist(earthquake$Magnitude,5)

Scan按更詳細(xì)的設(shè)置讀取文件，更加接近于底層。其調(diào)用格式為：

scan(file = "", what = double(), nmax = -1, n = -1, sep = "",

quote = if(identical(sep, "\n")) "" else "'\"", dec = ".",

skip = 0, nlines = 0, na.strings = "NA",

flush = FALSE, fill = FALSE, strip.white = FALSE,

quiet = FALSE, blank.lines.skip = TRUE, multi.line = TRUE,

comment.char = "", allowEscapes = FALSE,

fileEncoding = "", encoding = "unknown", text)

我們讀取近一個(gè)星期的全球地震數(shù)據(jù)

> scan(file="h:/my_docs/eqweek.csv",skip=1,sep=",",what=list("",0,0,0,mag=0,"",0,0,0,0,"","",0))->eq

Read 1075 records

分5個(gè)區(qū)間統(tǒng)計(jì)分布頻率

> table(factor(cut(eq$mag,5)))

(0.995,2.1]   (2.1,3.2]   (3.2,4.3]   (4.3,5.4]  (5.4,6.51] 

693         200          46         126          10 

可使用edit()函數(shù)編輯數(shù)據(jù)集后存為另一個(gè)數(shù)據(jù)集

t(eq)->eq2

19)分析數(shù)據(jù)集

繼續(xù)以全球近一周地震數(shù)據(jù)為例。我們先將變量放到搜索路徑上

> attach(earthquake)

先分析一下地震震深：

> summary(Depth)

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 

0.10    5.80   12.15   30.82   38.00  630.70      39 

Min表示地震震深最小值，Max表示最大值，Median為中位數(shù)，Mean為平均值。

我們?cè)囍鴱南旅娴?a href='/map/sandiantu/' style='color:#000;font-size:inherit;'>散點(diǎn)圖中觀察一下地震震深與震級(jí)的關(guān)系：

Depth是震深，Magnitude是震級(jí)，很有意思的是表面上看過去一周中Depth和Magnitude之間沒有關(guān)系，仔細(xì)觀察后這個(gè)圖，發(fā)現(xiàn)一個(gè)有趣的結(jié)果：當(dāng)震深超過300后，震級(jí)都接近5或在5以上，而在300以內(nèi)時(shí)，震級(jí)并不確定?？梢宰鲫P(guān)于震深的直方圖

hist(Depth)

這些只是根據(jù)一個(gè)星期的數(shù)據(jù)分析的結(jié)果，不一定就代表真正的答案。

lines函數(shù)可完成畫線

比如說我們繪制一個(gè)(10,40)、(20,50)、(30,60)的散點(diǎn)圖，并將點(diǎn)連成線

> plot(c(10,20,30),c(40,50,60))

> lines(c(10,20,30),c(40,50,60))

Fivenum函數(shù)返回以下數(shù)據(jù)：minimum, lower-hinge, median, upper-hinge, maximum

> fivenum(Magnitude)

[1] 1.0 1.3 1.7 2.5 6.5

表示震級(jí)最小為1.0，最大為6.5，中位數(shù)為1.3，通過1.3將一組數(shù)據(jù)分為上下兩組，然后再計(jì)算上下兩組的中位數(shù)1.3與2.5。

rug函數(shù)顯示實(shí)際的數(shù)據(jù)點(diǎn)

>  hist(Magnitude)

> rug(Magnitude)

利用直方圖估計(jì)密度函數(shù)存在密度函數(shù)是不平滑的，密度函數(shù)受子區(qū)間寬度影響很大、當(dāng)數(shù)據(jù)維數(shù)超過2維時(shí)有局限性等問題，因此基于核密度估計(jì)的方法可解決這些問題。

核密度估計(jì)又叫核函數(shù)估計(jì)(kernel density estimation)，是在概率論中用來估計(jì)未知的密度函數(shù)，屬于非參數(shù)檢驗(yàn)方法之一，由Rosenblatt (1955)和Emanuel Parzen(1962)提出，又名Parzen窗（Parzen window）。Ruppert和Cline基于數(shù)據(jù)集密度函數(shù)聚類算法提出修訂的核密度估計(jì)方法。

我們使用density函數(shù)進(jìn)行核密度估計(jì)：

> hist(Magnitude,prob=TRUE)

> lines(density(Magnitude))

累積分布函數(shù)能完整描述一個(gè)實(shí)數(shù)隨機(jī)變量X的概率分布，是概率密度函數(shù)的積分 ，與概率密度函數(shù)相對(duì)，定義為隨機(jī)變量小于或者等于某個(gè)數(shù)值的概率P（X<=x），即：F(x) = P(X<=x)

Ecdf函數(shù)完成累積分布函數(shù)的計(jì)算，我們計(jì)算一下震級(jí)的累積分布

> plot(ecdf(Magnitude),do.points=FALSE,verticals=TRUE)

19)R條件、循環(huán)語句

可以使用for做為循環(huán)：

> z<-c()

> x<-(1:10)

> y<-(11:20)

> for (i in 1:length(x)){

+ z[i]=x[i]^2+y[i]^2

+ }

> z

[1] 122 148 178 212 250 292 338 388 442 500

使用while做為循環(huán)語句

> x<-c(1:10)

> i=1

> while (x[i]^2<10)

+ {

+ i=i+1

+ x[i]=x[i]^2

+ }

> x

[1]  1  4  3  4  5  6  7  8  9 10

使用if條件語句:

> z<-c()

> x<-(1:10)

> y<-(11:20)

> for (i in 1:length(x)){

+ if ((x[i]^3>y[i]^2)) 

+ z[i]=x[i]^3

+ else

+ z[i]=y[i]^2

+ }

> z

[1]  121  144  169  196  225  256  343  512  729 1000

19)定義函數(shù)與操作符

有限維空間中，已知向量的坐標(biāo)，就可以知道它的模長(zhǎng)。

我們定義一個(gè)求3維向量模長(zhǎng)的函數(shù)：

> vector_length<-function(x1,x2,x3){

+ vlength<-sqrt (x1^2+x2^2+x3^2)

+ vlength

+ }

> vector_length(12,33,19)

[1] 39.92493

我們?cè)俣x一個(gè)n維的求向量模長(zhǎng)的函數(shù)

> vectorn_length<-function(x){

+ temp<-0

+ for (i in 1:length(x)){

+ temp<-temp+x[i]^2

+ }

+ vlength<-sqrt(temp)

+ vlength

+ }

> vectorn_length(c(11,22,33,44,55))

[1] 81.57818

> vectorn_length(c(11,22,33,55))

[1] 68.69498

操作符的定義使用%符號(hào)%的方式定義，實(shí)際使用時(shí)%也屬于操作符的一部分。我們定義一個(gè)操作符求歐氏距離。

歐氏距離（Euclid Distance）是在n維空間中兩個(gè)點(diǎn)之間的真實(shí)距離，n維歐氏空間的每個(gè)點(diǎn)可以表示為 (x[1]，x[2]，…，x[n]) ，兩個(gè)點(diǎn) X = (x[1]，x[2]，…，x[n]) 和 Y= (y[1]，y[2]，…，y[n]) 之間的距離 d(X，Y) 定義為下面的公式:

d(X，Y) =sqrt   (∑( ( x[i] - y[i] )^2 ))  其中 i = 1，2，…，n以%~%為操作符：

> "%~%"<-function(x1,x2){

+ temp<-0

+ for (i in 1:length(x1)){

+ temp<-temp+(x1[i]-x2[i])^2

+ }

+ edis<-sqrt(temp)

+ edis

+ }

> c(1,2,3) %~% c(5,6,7)

[1] 6.928203

指定參數(shù)的默認(rèn)值，在函數(shù)調(diào)用時(shí)可以省略該參數(shù)，通過在function隨后的參數(shù)表中直接指定

> vectorn_length<-function(x=c(1,1,1)){

+ temp<-0

+ for (i in 1:length(x)){

+ temp<-temp+x[i]^2

+ }

+ vlength<-sqrt(temp)

+ vlength

+ }

> vectorn_length()

[1] 1.732051

> vectorn_length(c(11,11,11))

[1] 19.05256

(3)不定數(shù)量的函數(shù)參數(shù)

> mycount<-function(...){

+ temp=0

+ for (i in c(...)){

+ temp=temp+1

+ }

+ temp

+ }

> mycount(11,22,33)

[1] 3

> mycount(11,22,33,66)

[1] 4

> mycount(11,22,66)

[1] 3

(4)內(nèi)嵌函數(shù)

允許在函數(shù)內(nèi)定義函數(shù)，使用方法和定義函數(shù)一樣，只不過整個(gè)函數(shù)體要定義在另一個(gè)函數(shù)體內(nèi) 。

5、統(tǒng)計(jì)模型基本應(yīng)用

(1)分析數(shù)據(jù)分布情況

我們分析一下某老年常見病數(shù)據(jù)，下面是數(shù)據(jù)格式

加載數(shù)據(jù)后，查看老年病的老人年齡分布情況及概率密度分布

> read.table("h:/my_docs/aged_patients.csv",header=TRUE,sep=",")->agedpatients

> hist(agedpatients$年齡)

> hist(agedpatients$年齡,prob=TRUE)

> lines(density(agedpatients$年齡))

> rug(agedpatients$年齡)

> min(agedpatients$年齡)

[1] 60

> max(agedpatients$年齡)

[1] 90

> mean(agedpatients$年齡)

[1] 69.09839

> var(agedpatients$年齡)

[1] 38.60801

63~72這個(gè)階段的老人注意堅(jiān)持運(yùn)動(dòng)，保持健康，這個(gè)年齡段是這2類老年病的高發(fā)時(shí)間段， 發(fā)生概率較大。從60到90的老人都有可能患上這2類老年病，年紀(jì)的平均值是69歲，標(biāo)準(zhǔn)差為38.60801比較大，看來這2個(gè)老年病在老人中發(fā)生的比較普遍。

我們?cè)侔茨挲g分類匯總腫瘤和急腹癥的數(shù)量

attach(agedpatients)

> tapply(腫瘤,年齡,sum)

60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 87 88 90 

19 26 26 28 21 35 37 22 21 23 29 26 31 16 19 12 13 19  8  5  4  3  7  2  1  1  0  3  1 

> tapply(急腹癥,年齡,sum)

60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 87 88 90 

2  0  0  2  0  0  1  0  1  2  4  1  1  0  3  1  2  2  0  2  0  1  0  1  0  1  2  0  0 

分年齡統(tǒng)計(jì)腫瘤患者的數(shù)量

> table(factor(cut(agedpatients$年齡[agedpatients$腫瘤==1],5)))

(60,66] (66,72] (72,78] (78,84] (84,90] 

155     158     118      22       5 

分年齡統(tǒng)計(jì)急腹癥患者的數(shù)量

> table(factor(cut(agedpatients$年齡[agedpatients$急腹癥==1],5)))

(60,65.4] (65.4,70.8] (70.8,76.2] (76.2,81.6]   (81.6,87] 

4           8           8           5           4