R语言与数据的预处理-CDA数据分析师官网

热线电话：13121318867

R语言与数据的预处理

2016-10-08

R语言与数据的预处理

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

一、获取重要数据

在导入大规模数据时，我们通常需要知道数据中的关键内容：最值，均值，离差，分位数，原点矩，离差，方差等。在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可以统计出数据中完全相同的数据个数。例如对《全宋词》中暴力拆解（两个相邻字算一词）词语使用数目的统计程序如下：

[plain] view plaincopyprint?

 l=scan(“Ci.txt”,”character”,sep=”\n”);

 l.len=nchar(l);

 ci=l;

 sentences=strsplit(ci,”，|。|！|？|、”);# 句子用标点符号分割。

 sentences=unlist(sentences);

 sentences=sentences[sentences!=””];

 s.len=nchar(sentences);#单句太长了说明有可能是错误的字符，去除掉。

 sentences=sentences[s.len<=10];

 s.len=nchar(sentences);

 splitwords=function(x,x.len)substring(x,1:(x.len-1),2:x.len);

 words=mapply(splitwords,sentences,s.len,SIMPLIFY=TRUE,USE.NAMES=FALSE);

 words=unlist(words);

 words.freq=table(words);#词频统计

 words.freq=sort(words.freq,decreasing=TRUE);

 data.frame(Word=names(words.freq[1:100]),Freq=as.integer(words.freq[1:100]));

而对于一堆数，我们按区间做的时候，就还需要函数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代码如下：

[plain] view plaincopyprint?

 pois.proc= function(T, lambda){

 S = 0

 I =rpois(1, lambda*T) #产生t 泊松分布，这里调用R 内置的泊松函数避免循环。

 U =runif(I)

 S =sort(T * U) #排序产生顺序统计量的思想

 list(I= I, S = S)

 }

 broken.proc= function(k, u= 3000, c= 120){

 n =1000 #模拟到时刻 1000 为止的破产情况

 M =pois.proc(n, 10)

 N =pois.proc(n, 1)

 U = u #初始盈余

 X = 0

 result=0

 A =sort(c(M$S, N$S)) #M$S和 N$S 是保单和理赔达到时刻

 for(iin 1:length(A)){

 if(any(A[i]==N$S)== 0)

 U=U+c

 else {

 X[i] =rexp(1, rate=1/1000)

 U = U -X[i] #减去这个随机值

 if(U< 0){ #判断盈余是否小于0（保单到达的时候不需要判断）

 result<-A[i] #盈余小于 0时，记录这个理赔到达（破产）的刻

 break}

 }

 }

 if(U>= 0){ #如果 for循环没有中断，判最终的盈余其实肯定非负

 result= n + 200

 } #给结果赋值一个明显比模拟时刻大的数据，表示未破产

 return(result) #返回最终结果

 }

 #根据这个破产过程可以模拟保险人的频数和频率：

 simulation= function(n=100){ #定义一个重复模拟破产过程的函数

 t =numeric(n)

 for(iin 1:n){

 t[i] =broken.proc(i)} #产生 n次破产或者代表未破产的时刻

 return(t)}

 time=simulation(n= 1200)

 rangetime= time[time!=1200]

 breakratio= length(rangetime)/length(time);

 breakratio

 break.points<-c(0,10,20,30,40,50,100,200,300,400,500,1000,1200)

 table(cut(time,breaks=break.points))

 hist(rangetime,breaks = 50, xlab=’broken time’,xlim = c(0, max(rangetime)),main = ‘Histogramof Broken time’)

用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)

运用举例：

[plain] view plaincopyprint?

 library(actuar)

 z=rnorm(10000)

 break.points<-c(-Inf,-3,-2,-1,0,1,2,3,Inf)

 tz<-table(cut(z,breaks=break.points))

 tz

 zz<-grouped.data(Group=break.points,freq=as.matrix(tz))

 zz

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

[plain] view plaincopyprint?

 mean(zz)

 mean(zz[c(2:7),])

 mean(z)

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

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

[plain] view plaincopyprint?

 EDA <- function (x)

 {

 par(mfrow=c(2,2)) # 同时做4个图

 hist(x) # 直方图

 dotchart(x) # 点图

 boxplot(x,horizontal=T) # 箱式图

 qqnorm(x);qqline(x) # 正态概率图

 par(mfrow=c(1,1)) # 恢复单图

 }

三、正态检验与经验分布

对于数据的分布估计经验分布是一个非常好的估计。在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

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

[plain] view plaincopyprint?

 NormTest<-function(data){

 library(fBasics)

 library(nortest)

 udata<-unique(data)

 result<-list()

 result$D<-dagoTest(data)

 result$jB<-jarqueberaTest(data)

 result$SW<-shapiroTest(data)

 result$lillie<-lillie.test(data)

 result$ad<-ad.test(data)

 result$cvm<-cvm.test(data)

 result$sf<-sf.test(data)

 return(result)

 }

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

解答如下：（结果以注释形式标明）

[plain] view plaincopyprint?

 v<-c(57,203,383,525,532,408,273,139,45,27,16)

 chisq.test(v)#p<0.05,认为检验总体是否与给定的p相同，p缺省表示等可能性检验

 #验证V的分布是否为poission分布

 x<-0:10

 options(digits=3)

 likely<-function(lamda=3){

 -sum(y*dpois(x,lamda=lamda,log=T))

 }

 library(stats4)

 mle(likely)

 chisq.fit<-function(x,y,r){

 options(digits=4)

 result<-list()

 n<-sum(y)

 prob<-dpois(x,3.87,log=F)

 y<-c(y,0)

 m<-length(y)

 prob<-c(prob,1-sum(prob))

 result$chisq<-sum((y-n*prob)^2/(n*prob))

 result$p.value<-pchisq(result$chisq,m-r-1,lower.tail=F)

 result

 }

 chisq.fit(x,v,1)#p<0.05 拒绝假设分布

CDA数据分析师考试相关入口一览（建议收藏）：

▷ 想报名CDA认证考试，点击>>> “CDA报名” 了解CDA考试详情；