做最具实战的
数据科学职业社区

R语言与数据的预处理

在面对大规模数据时,对数据预处理,获取基本信息是十分必要的。今天分享的就是数据预处理的一些东西。

uFR9f

一、获取重要数据

在导入大规模数据时,我们通常需要知道数据中的关键内容:最值,均值,离差,分位数,原点矩,离差,方差等。在R中常用的函数与作用整理如下:

统计函数 作用
Max 返回数据的最大值
Min 返回数据的最小值
Which.max 返回最大值的下标
Which.min 返回最小值的下标
Mean 求均值
Median 求中位数
mad 求离差
Var 求方差(总体方差)
Sd 求标准差
Range 返回【最小值,最大值】
Quantile 求分位数
Summary 返回五数概括与均值
Finenum 五数概括(最值,上下四分位数,中位数)
Sort 排序(默认升序,decreasing=T时为降序)
Order 排序(默认升序,decreasing=T时为降序)
Sum 求和
length 求数据个数
emm Actuar包中求k阶原点矩
skewness Fbasic包中求偏度
kurtosis Fbasics包中求峰度

注:对象为分组数据,矩阵时返回的不是整体的方差,均值,而是每一列(组)的方差均值其余变量类似。

二、直方图与频数统计

对于数据分布的认识,在大规模时有必要使用直方图。在R语言中,直方图的函数调用为:

hist(x, breaks = "Sturges",
freq = NULL, probability = !freq,
include.lowest = TRUE, right = TRUE,
density = NULL, angle = 45, col = NULL, border = NULL,
main= paste("Histogram of" , xname),
xlim = range(breaks), ylim = NULL,
xlab = xname, ylab,
axes = TRUE, plot = TRUE, labels = FALSE,
nclass = NULL, warn.unused = TRUE, ...)

这里值得一提的是,分组参数breaks默认使用史特吉斯(Sturges)公式,根据测定数n 来计算组距数k,公式为:k=1+3.32 logn。当然也可以自己设定一个数组来决定分组。(举例参见《R语言绘图学习笔记》)

说完频率分布直方图,我们还有频率分布直方表。对于数据的统计,函数table可以统计出数据中完全相同的数据个数。例如对《全宋词》中暴力拆解(两个相邻字算一词)词语使用数目的统计程序如下:

  1. <span style=”font-size:18px;”>l=scan(“Ci.txt”,”character”,sep=”\n”);
  2. l.len=nchar(l);
  3. ci=l;
  4. sentences=strsplit(ci,”,|。|!|?|、”);# 句子用标点符号分割。
  5. sentences=unlist(sentences);
  6. sentences=sentences[sentences!=””];
  7. s.len=nchar(sentences);#单句太长了说明有可能是错误的字符,去除掉。
  8. sentences=sentences[s.len<=10];
  9. s.len=nchar(sentences);
  10. splitwords=function(x,x.len)substring(x,1:(x.len-1),2:x.len);
  11. words=mapply(splitwords,sentences,s.len,SIMPLIFY=TRUE,USE.NAMES=FALSE);
  12. words=unlist(words);
  13. words.freq=table(words);#词频统计
  14. words.freq=sort(words.freq,decreasing=TRUE);
  15. data.frame(Word=names(words.freq[1:100]),Freq=as.integer(words.freq[1:100]));</span>

而对于一堆数,我们按区间做的时候,就还需要函数cut.调用格式如下:

cut(x, breaks, labels = NULL,
include.lowest = FALSE, right = TRUE, dig.lab = 3,
ordered_result = FALSE, ...)

举一个具体例子,某一款保险产品,假设保单到达的速率为10张/天,理赔发生的速率为 1次/天。假设每张保单价格c=120,理赔额服从参数为v=1/1000 (以c*lambda1=1.2*lambda2/v设定)的指数分布。设定初始u=3000时,计算到第1000天为止发生破产的概率。(案例摘自《复 合泊松过程模型的推广和在R语言环境下的随机模拟》 )

破产过程的R代码如下:

  1. <span style=”font-size:18px;”>pois.proc= function(T, lambda){
  2. S = 0
  3. I =rpois(1, lambda*T) #产生t 泊松分布,这里调用R 内置的泊松函数避免循环。
  4. U =runif(I)
  5. S =sort(T * U) #排序产生顺序统计量的思想
  6. list(I= I, S = S)
  7. }
  8. broken.proc= function(k, u= 3000, c= 120){
  9. n =1000  #模拟到时刻 1000 为止的破产情况
  10. M =pois.proc(n, 10)
  11. N =pois.proc(n, 1)
  12. U = u   #初始盈余
  13. X = 0
  14. result=0
  15. A =sort(c(M$S, N$S))  #M$S和 N$S 是保单和理赔达到时刻
  16. for(iin 1:length(A)){
  17. if(any(A[i]==N$S)== 0)
  18. U=U+c
  19. else {
  20. X[i] =rexp(1, rate=1/1000)
  21. U = U -X[i]   #减去这个随机值
  22. if(U< 0){    #判断盈余是否小于0(保单到达的时候不需要判断)
  23. result<-A[i]   #盈余小于 0时,记录这个理赔到达(破产)的刻
  24. break}
  25. }
  26. }
  27. if(U>= 0){   #如果 for循环没有中断,判最终的盈余其实肯定非负
  28. result= n + 200
  29. }  #给结果赋值一个明显比模拟时刻大的数据,表示未破产
  30. return(result)   #返回最终结果
  31. }
  32. #根据这个破产过程可以模拟保险人的频数和频率:
  33. simulation= function(n=100){ #定义一个重复模拟破产过程的函数
  34. t =numeric(n)
  35. for(iin 1:n){
  36. t[i] =broken.proc(i)}   #产生 n次破产或者代表未破产的时刻
  37. return(t)}
  38. time=simulation(n= 1200)
  39. rangetime= time[time!=1200]
  40. breakratio= length(rangetime)/length(time);
  41. breakratio
  42. break.points<-c(0,10,20,30,40,50,100,200,300,400,500,1000,1200)
  43. table(cut(time,breaks=break.points))
  44. hist(rangetime,breaks = 50, xlab=’broken time’,xlim = c(0, max(rangetime)),main = ‘Histogramof Broken time’)</span>

用R 语言模拟了1200 次,最终结果 1200 次中破产 628 次,破产率大概 52.3% 。输出各阶段破产时刻 频数和率结果如下:

区间

频数

(0,10]

389

(10,20]

89

(20,30]

45

(30,40]

28

(40,50]

16

(50,100]

36

(100,200]

17

(200,300]

6

(300,400]

2

(400,500]

0

(500,1000]

0

(1000,1200]

572

对于一些数据我们可能直接录入的是频率分布直方表,那么actuar包中提供了一个有用的数据结构grouped.data。调用格式:

grouped.data(..., right = TRUE, row.names = NULL, check.rows = FALSE,
             check.names = TRUE)

运用举例:

  1. <span style=”font-size:18px;”>library(actuar)
  2. z=rnorm(10000)
  3. break.points<-c(-Inf,-3,-2,-1,0,1,2,3,Inf)
  4. tz<-table(cut(z,breaks=break.points))
  5. tz
  6. zz<-grouped.data(Group=break.points,freq=as.matrix(tz))
  7. zz</span>

对比一下下面的输出结果,我们发现分组数据的均值计算与总体数据计算方法是不一样的。

  1. <span style=”font-size:18px;”>mean(zz)
  2. mean(zz[c(2:7),])
  3. mean(z)</span>

注:函数elev()可以计算有限期望值,可以避免mean(zz)不存在的尴尬。

当然对于数据的直观分析R提供的函数有许多,我们将常见的函数汇总如下:

  1. EDA <- function (x)
  2.   {
  3.     par(mfrow=c(2,2))              # 同时做4个图
  4.     hist(x)                        # 直方图
  5.     dotchart(x)                    # 点图
  6.     boxplot(x,horizontal=T)        # 箱式图
  7.     qqnorm(x);qqline(x)            # 正态概率图
  8.     par(mfrow=c(1,1))              # 恢复单图
  9.   }

三、正态检验与经验分布

对于数据的分布估计经验分布是一个非常好的估计。在actuar包中函数ogive给出的实现:

ogive(x, y = NULL, …)

## S3 method for class ‘ogive’

print(x, digits = getOption(“digits”) – 2, …)

## S3 method for class ‘ogive’

summary(object, …)

## S3 method for class ‘ogive’

knots(Fn, …)

## S3 method for class ‘ogive’

plot(x, main = NULL, xlab = “x”, ylab = “F(x)”, …)

还是以上面的例子数据zz为例:

ogive(zz)

plot(ogive(zz))

输出结果:

Ogive forgrouped data

Call:ogive(zz)

x =  -Inf,     -3,     -2, …,      3,    Inf

F(x) =     0, 0.0011, 0.0229,  …,0.9985,      1

13Y513Z9420-12643

由于大数定律的存在,很多情况下,正态性检验是十分有必要的一个分布检验,在R中提供的正态性检验可以汇总为下面的一个正态检验函数:

  1. <span style=”font-size:18px;”>NormTest<-function(data){
  2. library(fBasics)
  3. library(nortest)
  4. udata<-unique(data)
  5. result<-list()
  6. result$D<-dagoTest(data)
  7. result$jB<-jarqueberaTest(data)
  8. result$SW<-shapiroTest(data)
  9. result$lillie<-lillie.test(data)
  10. result$ad<-ad.test(data)
  11. result$cvm<-cvm.test(data)
  12. result$sf<-sf.test(data)
  13. return(result)
  14. }</span>

对于分布的检验还有卡方检验,柯尔莫哥洛夫检验等,在R中也有实现函数chisq.test()等。我们同样以一个例子来说明:

13Y513910140-21293

例子摘自:王兆军《数理统计讲义》习题6.3

解答如下:(结果以注释形式标明)

  1. <span style=”font-size:18px;”>v<-c(57,203,383,525,532,408,273,139,45,27,16)
  2. chisq.test(v)#p<0.05,认为检验总体是否与给定的p相同,p缺省表示等可能性检验
  3. #验证V的分布是否为poission分布
  4. x<-0:10
  5. options(digits=3)
  6. likely<-function(lamda=3){
  7. -sum(y*dpois(x,lamda=lamda,log=T))
  8. }
  9. library(stats4)
  10. mle(likely)
  11. chisq.fit<-function(x,y,r){
  12. options(digits=4)
  13. result<-list()
  14. n<-sum(y)
  15. prob<-dpois(x,3.87,log=F)
  16. y<-c(y,0)
  17. m<-length(y)
  18. prob<-c(prob,1-sum(prob))
  19. result$chisq<-sum((y-n*prob)^2/(n*prob))
  20. result$p.value<-pchisq(result$chisq,m-r-1,lower.tail=F)
  21. result
  22. }
  23. chisq.fit(x,v,1)#p<0.05 拒绝假设分布</span>
注:数据科学网遵循行业规范,任何转载的稿件都会明确标注作者和来源,若标注有误或遗漏,请联系主编邮箱:admin@idatacamp.com
数据科学 » R语言与数据的预处理
分享到: 更多 (0)

吐槽 抢沙发

评论前必须登录!