R语言基础编程技巧汇编 - 9

1. qq图的做法

qq图是正态分位数图，纵坐标是变量的取值，关键是横坐标。

qqplot全名应该是正态分位数图，横坐标的做法：

首先把变量按从小到大的顺序排列，计算变量的长度，即总共有多少个取值，再按顺序计算变量的所有取值的累积百分比，所谓的累积百分比，也就是可以看成是累积概率，比如有10个值，按照从小到大的顺序，第一个值的排序是1，那么他的所占的百分比就是10%，紧接着后一个值所占的百分比也会是10%，但是累积概率值为20%，依次往后计算，因为最后一个值的累积百分比是100%，即等于1，这个值如果计算它的正态分布概率的分位数的话，是无限大的，因此需要对这个值进行修正一下，就是因为这一个值无限大，所以对全体计算出来的累积百分比减去一个适当小的数，修正后的累积百分比与原百分比相差不多，但是回避了最后一个值是1而无法计算的问题。

有了累积百分比之后，相对应的就是累积的概率值。将累积概率值修正后，即得到累积概率，比如以10个值为例，第一个值的累积概率为0.05，查正态分布表，0.05的累积概率，对应的正态分布的Z值为-1.64，这样一次计算，所得的Z值，就是qqplot的横坐标数据。下面以10个数据和30个数据为例说明。

my.qqplot <-function(y){

op <- par(mfrow= c(1, 1))

N <- length(y)

n <- seq(1, N)

xais <-qnorm((n - (.5*N) /N)/ N)

#####中间三句可选，只是为了输出计算过程######

mid <-cbind(sort(y), n, n/N, (n-(.5*N)/N)/N , xais)

colnames(mid)<- c("y", "rank", "cumpercent","adj-cumper","xaix")

print(mid)

#####中间三句可选，只是为了输出计算过程######

par(mfrow =c(2,1))

qqnorm(y)

plot(sort(y) ~xais, main = 'my qqplot')

par(op)

}

y <- rnorm(10,mean = 20, s = 50)

my.qqplot(y)

y rankcumpercentadj-cumper xaix

[1,] 2.877321 1 0.1 0.05 -1.6448536

[2,] 6.930063 2 0.2 0.15 -1.0364334

[3,]16.461444 3 0.3 0.25-0.6744898

[4,]36.130825 4 0.4 0.35 -0.3853205

[5,]40.477883 5 0.5 0.45 -0.1256613

[6,]50.534636 6 0.6 0.55 0.1256613

[7,]53.425025 7 0.7 0.65 0.3853205

[8,]54.554269 8 0.8 0.75 0.6744898

[9,]120.496268 9 0.9 0.85 1.0364334

[10,]125.290253 10 1.0 0.95 1.6448536

2. read.csv读取一定数量的行

a=read.csv("C:/a.csv",nrows=500)

3. 数据高位用0补足

有1到100的一组数据，变成6位的数据，不足的用0来补，希望得到000001，........，000099，000100，解决方法：

a<- sprintf("%06d",1:100)

4. 用XML包抓取网页

library(XML)

dat = readLines("http://www.waaku.com/InfoView6792299.html")

doc = htmlTreeParse(dat, asText =TRUE, useInternalNodes =TRUE)

xpathSApply(xmlRoot(doc),
"//div[@id='new6']",xmlValue)

5. 处理命令行参数

optparse包和getopt包

6. 数据表中的变量访问方法

head(sleep)

BodyWgt BrainWgt NonD Dream Sleep Span GestPred Exp Danger

1 6654.000 5712.0 NA NA 3.3 38.6 645 3 5 3

2 1.000 6.6 6.3 2.0 8.3 4.5 42 3 1 3

3 3.385 44.5 NA NA 12.5 14.0 60 1 1 1

4 0.920 5.7 NA NA 16.5 NA 25 5 2 3

5 2547.000 4603.0 2.1 1.8 3.9 69.0 624 3 5 4

6 10.550 179.5 9.1 0.7 9.8 27.0 180 4 4 4

以R语言自带的数据框sleep为例，访问sleep中第一列BodyWgt,有几种可能的方法：

（1）中括号+列名带引号，即sleep['BodyWgt']，将得到一个类型为数据框的访问结果，如下所示：

sleep['BodyWgt']

BodyWgt

1 6654.000

2 1.000

3 3.385

4 0.920

5 2547.000

6 10.550

str(sleep['BodyWgt'])

'data.frame': 62 obs.of 1 variable:

$ BodyWgt: num 6654 1 3.38 0.92 2547 ...

方法1的变形是sleep[1]，即访问sleep的第一列，得到的类型也是数据框。

（2）$+列名，即sleep$BodyWgt，列名不加引号，得到的数据类型是向量。

sleep$BodyWgt

[1] 6654.000 1.000 3.385 0.920 2547.000 10.550

str(sleep$BodyWgt)

num [1:62] 6654 1 3.38 0.92 2547 ...

（3）中括号+双下标，即sleep[,1]，意义为：访问第一列的所有的行，得到一个向量。

sleep[,1]

[1] 6654.000 1.000 3.385 0.920 2547.000 10.550

str(sleep[,1])

num [1:62] 6654 1 3.38 0.92 2547...

7. 时间格式

options(scipen=200)

t <- c(90000,93000,100000,102000)

tt=prettyNum(t,":",big.interval =2L)

strptime(tt,'%H:%M:%S')

8. scale函数做数据中心化和标准化

scale()

9. barplot上标出频数

x <- rep(1:4)

bp <- barplot(x)

text(bp,x/2,x)

10. rvest包网页抓取示例

下面用抓取豆瓣网电影数据为例，演示一下rvest包。

library(rvest)

url = 'http://movie.douban.com/top250?format=text'

web = html(url,encoding="UTF-8")

score = html_text(web %>%html_nodes("span") %>% html_nodes("em"))

year = web %>%html_nodes("p")

year= html_text(year[seq(1,50,2)])

gy = gregexpr('[0-9]{4}',year)

gd = gregexpr("导演",year)

gz = gregexpr("主演",year)

time = sapply(1:length(gy),function(i)substr(year[i],gy[[i]],gy[[i]]+attr(gy[[i]],'match.length')-1))

direct = sapply(1:length(gy),function(i)substr(year[i],gd[[i]] + 4,gz[[i]] - 4))

actor = sapply(1:length(gy),function(i)substr(year[i],gz[[i]] + 4,gy[[i]] - 4))

ga = gregexpr('\n',actor)

actor = sapply(1:length(gy),function(i)substr(actor[i],1,ga[[i]] - 1))

names = web %>%html_nodes("img") %>% html_attr("alt")

names =iconv(names,"utf-8","gbk")

rates = html_text(web %>%html_nodes("span"))

rates = rates[grep("人评价",rates)]

rates = gsub("人评价","",rates)

head(names)

head(time)

head(score)

head(rates)

head(direct)

head(actor)

11. 显示数学公式

plot(0, 0, type = "n", axes = F,xlab="", ylab = "")

text(0, 0,bquote(italic(mu[tilde(s)
^2](s
^2))), cex = 12)

12. 图例加在图形框以外

plot(c(1,1),
c(0.8,1.2), pch =20:21)

legend(0.6,1.27,
c("A","B"), lty =1:2,
pch =20:21, ncol =2,
bty ="n", xpd =TRUE)

13. 开n次方

X^(1/n)

14. 用rgl包绘制三维球体

library("rgl")

lat <- matrix(seq(90,-90, len=50)*pi/180,50, 50, byrow=TRUE)

long <- matrix(seq(-180, 180,len=50)*pi/180, 50, 50)

r <- 1

x <- -r*cos(lat)*cos(long)

y <- r*sin(lat)

z <- r*cos(lat)*sin(long)

open3d()

persp3d(x, y, z, col="white" )

15. 查看当前工作空间中有哪些变量

objects()