【數據分析】R語言中用自助法求統計量置信區間

當樣本不符合理論分布假設時,求樣本統計量的置信區間就成為一個難題。而自助法(Bootstrap)的思路是對原始樣本重復抽樣產生多個新樣本,針對每個樣本求取統計量,然後得到它的經驗分布,再通過求經驗分布的分位數來得到統計量的置信區間,這種方法不需要對統計量有任何理論分布的假設。一般認為,只要樣本具有代表性,採用自助法需要的原始樣本只要20-30個,重復抽樣1000次就能達到滿意的結果。

在R中進行自助法是利用boot擴展包,其流程如下:

  • 編寫一個求取統計量的自定義函數

  • 將上面的函數放入boot()函數中進行運算,得到自助法的結果

  • 用boot.ci()函數求取置信區間

讓我們用mtcars數據集來作為例子,我們可以將wt和disp作為自變量,mpg 作為因變量,進行回歸後能得到一系列回歸統計量。其中我們感興趣的是判定系數R-square,希望用自助法求它的95%置信區間。

首先定義求R-square的函數,注意其中的indices是必不可少的參數,另外一個參數代表樣本數據
————————
rsq=function(data,indices){
 d=data[indices,]
 fit=lm(formula=mpg~wt+disp,data=d)
 return(summary(fit)$r.square)
 }
————————
載入boot擴展包,將隨機種子設為1234,以方便得到相同的結果,再利用boot函數得到結果results,其中R表示重復抽樣得到1000個樣本
————————
library(boot)
set.seed(1234)
results=boot(data=mtcars,statistic=rsq,R=1000)
print(results)
plot(results)
————————

【數據分析】R語言中用自助法求統計量置信區間

results這個數據結構中包括了原始樣本的統計量(results$t0)和再抽樣樣本的統計量(results$t0),上圖左側的直方圖表示了再抽樣樣本的統計量的經驗分布,其中的虛線表示了原始樣本的統計量,從中可以觀察到偏差。右側QQ圖有助於判斷經驗分布是否正態。下面我們用boot.ci函數從結果中提取置信區間。

————————
boot.ci(results,conf=0.95,type=c('perc','bca'))
————————
其中conf表示置信水平,type表示了用何種算法來求區間,perc即使用百分位方法,bca表示adjusted bootstrap percentile,即對偏差進行了調整。結果如下:

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = results, conf = 0.95, type = c("perc", "bca"))

Intervals : 
Level     Percentile            BCa          
95%   ( 0.6838,  0.8833 )   ( 0.6344,  0.8549 )

來源:http://www.r-bloggers.com/lang/chinese/540

【數據分析】R語言中用自助法求統計量置信區間

1、回復「數據分析師」查看數據分析師系列文章
2、回復「案例」查看大數據案例系列文章
3、回復「徵信」查看相關徵信的系列文章
4、回復「可視化」查看可視化專題系列文章
5、回復「SPPS」查看SPSS系列文章
6、回復「答案」查看hadoop面試題題目及答案
7、回復「愛情」查看大數據與愛情的故事
8、回復「笑話」查看大數據系列笑話
9、回復「大數據1、大數據2、大數據3、大數據4」查看大數據歷史機遇連載

PPV課大數據ID: ppvke123 (長按可復制)

大數據人才的搖籃!專注大數據行業人才的培養。每日一課,大數據(EXCEL、SAS、SPSS、Hadoop、CDA)影片課程。大數據資訊,每日分享!數據咖—PPV課數據愛好者俱樂部!
【數據分析】R語言中用自助法求統計量置信區間

閱讀原文


關於作者:
「PPV課大數據」是一個大數據學習平台,數據分析、數據挖掘交流和分享的圈子。我們用數據說話,傳播正能量,執著探索大數據應用價值!

微信號:ppvke123