當樣本不符合理論分布假設時,求樣本統計量的置信區間就成為一個難題。而自助法(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)
————————
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
1、回復「數據分析師」查看數據分析師系列文章
2、回復「案例」查看大數據案例系列文章
3、回復「徵信」查看相關徵信的系列文章
4、回復「可視化」查看可視化專題系列文章
5、回復「SPPS」查看SPSS系列文章
6、回復「答案」查看hadoop面試題題目及答案
7、回復「愛情」查看大數據與愛情的故事
8、回復「笑話」查看大數據系列笑話
9、回復「大數據1、大數據2、大數據3、大數據4」查看大數據歷史機遇連載
關於作者:
「PPV課大數據」是一個大數據學習平台,數據分析、數據挖掘交流和分享的圈子。我們用數據說話,傳播正能量,執著探索大數據應用價值!
微信號:ppvke123