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)[n]^2](s[n]^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()