R語言與數(shù)據(jù)的預(yù)處理

在面對大規(guī)模數(shù)據(jù)時(shí)，對數(shù)據(jù)預(yù)處理，獲取基本信息是十分必要的。今天分享的就是數(shù)據(jù)預(yù)處理的一些東西。

一、獲取重要數(shù)據(jù)

在導(dǎo)入大規(guī)模數(shù)據(jù)時(shí)，我們通常需要知道數(shù)據(jù)中的關(guān)鍵內(nèi)容：最值，均值，離差，分位數(shù)，原點(diǎn)矩，離差，方差等。在R中常用的函數(shù)與作用整理如下：

統(tǒng)計(jì)函數(shù)
作用
Max
返回?cái)?shù)據(jù)的最大值
Min
返回?cái)?shù)據(jù)的最小值
Which.max
返回最大值的下標(biāo)
Which.min
返回最小值的下標(biāo)
Mean
求均值
Median
求中位數(shù)
mad
求離差
Var
求方差（總體方差）
Sd
求標(biāo)準(zhǔn)差
Range
返回【最小值，最大值】
Quantile
求分位數(shù)
Summary
返回五數(shù)概括與均值
Finenum
五數(shù)概括（最值，上下四分位數(shù)，中位數(shù)）
Sort
排序（默認(rèn)升序，decreasing=T時(shí)為降序）
Order
排序（默認(rèn)升序，decreasing=T時(shí)為降序）
Sum
求和
length
求數(shù)據(jù)個數(shù)
emm
Actuar包中求k階原點(diǎn)矩
skewness
Fbasic包中求偏度
kurtosis
Fbasics包中求峰度

注：對象為分組數(shù)據(jù)，矩陣時(shí)返回的不是整體的方差，均值，而是每一列（組）的方差均值其余變量類似。

二、直方圖與頻數(shù)統(tǒng)計(jì)

對于數(shù)據(jù)分布的認(rèn)識，在大規(guī)模時(shí)有必要使用直方圖。在R語言中，直方圖的函數(shù)調(diào)用為：

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, ...)

這里值得一提的是，分組參數(shù)breaks默認(rèn)使用史特吉斯（Sturges）公式，根據(jù)測定數(shù)n 來計(jì)算組距數(shù)k，公式為：k＝1+3.32 logn。當(dāng)然也可以自己設(shè)定一個數(shù)組來決定分組。（舉例參見《R語言繪圖學(xué)習(xí)筆記》）

說完頻率分布直方圖，我們還有頻率分布直方表。對于數(shù)據(jù)的統(tǒng)計(jì)，函數(shù)table可以統(tǒng)計(jì)出數(shù)據(jù)中完全相同的數(shù)據(jù)個數(shù)。例如對《全宋詞》中暴力拆解（兩個相鄰字算一詞）詞語使用數(shù)目的統(tǒng)計(jì)程序如下：

[plain] view plaincopyprint?

    <span style=”font-size:18px;”>l=scan(“Ci.txt”,”character”,sep=”\n”);

    l.len=nchar(l);

    ci=l;

    sentences=strsplit(ci,”，|。|！|？|、”);# 句子用標(biāo)點(diǎn)符號分割。

    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);#詞頻統(tǒng)計(jì)

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

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

而對于一堆數(shù)，我們按區(qū)間做的時(shí)候，就還需要函數(shù)cut.調(diào)用格式如下：

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

舉一個具體例子，某一款保險(xiǎn)產(chǎn)品，假設(shè)保單到達(dá)的速率為10張/天，理賠發(fā)生的速率為 1次/天。假設(shè)每張保單價(jià)格c=120，理賠額服從參數(shù)為v=1/1000 (以c*lambda1=1.2*lambda2/v設(shè)定)的指數(shù)分布。設(shè)定初始u=3000時(shí)，計(jì)算到第1000天為止發(fā)生破產(chǎn)的概率。（案例摘自《復(fù) 合泊松過程模型的推廣和在R語言環(huán)境下的隨機(jī)模擬》 ）

破產(chǎn)過程的R代碼如下：

[plain] view plaincopyprint?

    <span style=”font-size:18px;”>pois.proc= function(T, lambda){

    S = 0

    I =rpois(1, lambda*T) #產(chǎn)生t 泊松分布，這里調(diào)用R 內(nèi)置的泊松函數(shù)避免循環(huán)。

    U =runif(I)

    S =sort(T * U) #排序產(chǎn)生順序統(tǒng)計(jì)量的思想

    list(I= I, S = S)

    }


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

    n =1000  #模擬到時(shí)刻 1000 為止的破產(chǎn)情況

    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 是保單和理賠達(dá)到時(shí)刻

    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]   #減去這個隨機(jī)值

    if(U< 0){    #判斷盈余是否小于0（保單到達(dá)的時(shí)候不需要判斷）

    result<-A[i]   #盈余小于 0時(shí)，記錄這個理賠到達(dá)（破產(chǎn)）的刻

    break}

    }

    }

    if(U>= 0){   #如果 for循環(huán)沒有中斷，判最終的盈余其實(shí)肯定非負(fù)

    result= n + 200

    }  #給結(jié)果賦值一個明顯比模擬時(shí)刻大的數(shù)據(jù)，表示未破產(chǎn)

    return(result)   #返回最終結(jié)果

    }

    #根據(jù)這個破產(chǎn)過程可以模擬保險(xiǎn)人的頻數(shù)和頻率：

    simulation= function(n=100){ #定義一個重復(fù)模擬破產(chǎn)過程的函數(shù)

    t =numeric(n)

    for(iin 1:n){

    t[i] =broken.proc(i)}   #產(chǎn)生 n次破產(chǎn)或者代表未破產(chǎn)的時(shí)刻

    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’)</span>

用R 語言模擬了1200 次，最終結(jié)果 1200 次中破產(chǎn) 628 次，破產(chǎn)率大概 52.3% 。輸出各階段破產(chǎn)時(shí)刻 頻數(shù)和率結(jié)果如下：


區(qū)間
    

頻數(shù)
    
    
    

(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
    
    
    

對于一些數(shù)據(jù)我們可能直接錄入的是頻率分布直方表，那么actuar包中提供了一個有用的數(shù)據(jù)結(jié)構(gòu)grouped.data。調(diào)用格式：

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

check.names = TRUE)

運(yùn)用舉例：

[plain] view plaincopyprint?

    <span style=”font-size:18px;”>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</span>

對比一下下面的輸出結(jié)果，我們發(fā)現(xiàn)分組數(shù)據(jù)的均值計(jì)算與總體數(shù)據(jù)計(jì)算方法是不一樣的。

[plain] view plaincopyprint?

    <span style=”font-size:18px;”>mean(zz)

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

    mean(z)</span>

注：函數(shù)elev()可以計(jì)算有限期望值，可以避免mean（zz）不存在的尷尬。

當(dāng)然對于數(shù)據(jù)的直觀分析R提供的函數(shù)有許多，我們將常見的函數(shù)匯總?cè)缦拢?br />
[plain] view plaincopyprint?

    EDA <- function (x)

      {

        par(mfrow=c(2,2))              # 同時(shí)做4個圖

        hist(x)                        # 直方圖

        dotchart(x)                    # 點(diǎn)圖

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

        qqnorm(x);qqline(x)            # 正態(tài)概率圖

        par(mfrow=c(1,1))              # 恢復(fù)單圖

      }

三、正態(tài)檢驗(yàn)與經(jīng)驗(yàn)分布

對于數(shù)據(jù)的分布估計(jì)經(jīng)驗(yàn)分布是一個非常好的估計(jì)。在actuar包中函數(shù)ogive給出的實(shí)現(xiàn)：

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)”, …)

還是以上面的例子數(shù)據(jù)zz為例：

ogive(zz)

plot(ogive(zz))

輸出結(jié)果：

Ogive forgrouped data

Call:ogive(zz)

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

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

由于大數(shù)定律的存在，很多情況下，正態(tài)性檢驗(yàn)是十分有必要的一個分布檢驗(yàn)，在R中提供的正態(tài)性檢驗(yàn)可以匯總為下面的一個正態(tài)檢驗(yàn)函數(shù)：

對于分布的檢驗(yàn)還有卡方檢驗(yàn)，柯爾莫哥洛夫檢驗(yàn)等，在R中也有實(shí)現(xiàn)函數(shù)chisq.test()等。我們同樣以一個例子來說明：

解答如下：（結(jié)果以注釋形式標(biāo)明）

[plain] view plaincopyprint?

    <span style=”font-size:18px;”>v<-c(57,203,383,525,532,408,273,139,45,27,16)

    chisq.test(v)#p<0.05,認(rèn)為檢驗(yàn)總體是否與給定的p相同，p缺省表示等可能性檢驗(yàn)

    #驗(yàn)證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 拒絕假設(shè)分布</span>