R Resampling (リサンプリング)

ゼミでリサンプリングの手法を学んだので、知識と各手法の善し悪しについて整理しておく。主にこちらのサイトを参考にした。

リサンプリングの定義

初期サンプルから、新しいサンプルを作りだすことに基づいて得られる統計的推定の方法

リサンプリングの目的

1. 推定量のバイアスを低減する(Jackknife, Bootstrap)

2. 正規分布を仮定せずに、信頼区間を推定する(Bootstrap)

3. 正確な検定を行う(Permutation test)

4. モデルの精度検証を行う(Cross Validation)

代表的な手法

1. Jackknife

 - n個のデータのうち、k個を抽出し、残りn-k個のデータで統計量の算出やモデルの作成を行う。抽出するk個の全パターンについて、この操作を繰り返す

 - 外れ値の検出に用いられる

 - 一般的に分散や共分散の推定にも使われる

 - 推定量のバイアスを低減するための方法

 - 不偏推定量に対しては無意味(??⇒要調査)

 - 推定量が観測値の滑らかな関数であるのみ有効

2. Bootstrap

 -1つの標本から復元抽出を繰り返して大量の標本を生成し、それらの標本から推定値 を計算し、母集団の性質やモデルの推測の誤差などを分析する方法(こちらより引用)

3. Permutation test

 - 2群の差を検定する際によく用いられる方法

 - 2個それぞれのデータ数を固定したまま、データをランダムに入れ替える。その結果から統計値を算出し、オリジナルの統計値と比較することで、2つの集団に差があるかどうかを検定する。(2群に差があるのならば、オリジナル統計値よりも高い値が出る可能性が低い)

 - 対応のないt検定の代替手法と考えられる

 - 非復元抽出である

4. Cross Validation

 - 再現性の確認に用いられる(モデルの精度評価など)

 - データをk個の部分集合に分割し、そのうちのk-1個をtraining dataとし、残り1個をtest dataとし、training dataからモデルを作成し、test dataで精度を測定する。test dataを順に各部分集合に割り当て、この操作をk回繰り返す。こちらの説明が分かりやすい。

Rで実際にやってみる

1. Jackknife

適当にデータを仮定し、相関係数について外れ値を検出してみる。

#内閣府よりGDPの成長率を取得
GDP <- read.csv("http://www.esri.cao.go.jp/jp/sna/data/data_list/sokuhou/files/2013/qe133_2/__icsFiles/afieldfile/2013/12/05/ritu-jfy1332.csv",skip=8,header=FALSE)
GDP <- GDP[,1:3]
GDP <- GDP[-nrow(GDP),]
colnames(x) <- c("Fiscal Year","GDP(Expenditure Approach)","PrivateConsumption")

#GDPと個人消費(Private Consumption)の相関係数を計算
COR <- cor(GDP[,2],GDP[,3])
plot(GDP[,2],GDP[,3],main="GDPと個人商品の相関関係",xlab="GDP前年度比成長率(%)",ylab="個人消費前年度比成長率(%)")
text(COR,1,paste("correlation coefficient =",round(COR,3)),cex=1)

#Jackknife法で相関係数を計算
COR_J <- rep(0,nrow(GDP))

for (i in 1:nrow(GDP)) {
  COR_J[i] <- cor(GDP[-i,2],GDP[-i,3])
}

#相関係数が著しく異なる年を特定
plot(c(1995:2012),COR_J,main="GDPと個人消費の相関係数のJackknife推定値",ylab="相関係数",xlab="Jackknife法において除いた年度")
abline(h=COR,col="red")

結果はこちら。リーマンショックが起こった2008年度,2009年度のデータを除くと相関係数が大きく変化することがわかるので、これらを"外れ値"と疑うことが出来る。

f:id:nakhirot:20131216215437p:plain

2. Bootstrap

正規分布でない母集団を仮定し、Bootstrapサンプルから推定した平均値・分散と、母集団の平均・分散を比較してみる。

#Making triangular distribution
N <- 10000
x <- runif(N)
y <- numeric(N)

for (i in 1:N) {
  if(x[i] > 0.5) {
{y[i] <- 1 - sqrt(0.5*(1-x[i]))}
  } else {
    y[i] <- sqrt(0.5*x[i])
  }
}

t_mean <- mean(y)
hist(y,main="母集団の分布")

母集団の分布はこちら。平均は0.5となる。

f:id:nakhirot:20131216215902p:plain

次に、母集団(個数10,000)から、ブートストラップサンプル(個数1000)を抽出し、その平均を計算する。これを100回繰り返し、平均値の分布を調べてみる。

#Repeat k resampling with Bootstrap method with replacement
#Repeat k resampling with Bootstrap method without replacement
n <- 1000
k <- 100
sample1 <- sample2 <- matrix(0,nrow=k,ncol=n)
s_mean1 <- s_mean2 <- numeric(k)

for (i in 1:k) {
  sampleNO1 <- sample(N,n,replace=TRUE)
  sample1[i,] <- y[sampleNO1]
  s_mean1[i] <- mean(sample1[i,])
  sampleNO2 <- sample(N,n,replace=FALSE)
  sample2[i,] <- y[sampleNO2]
  s_mean2[i] <- mean(sample2[i,])
}

#Comparing true mean with means derived from the resamplings
hist(s_mean1,col="#ff00ff40",border="#ff00ff",ylim=c(0,40),main="Distribution of means derived from Bootstrap method")
hist(s_mean2,col="#0000ff40",border="#0000ff",ylim=c(0,40),add=TRUE)
legend("topleft",legend=c("with replacement","without replacement"),fill=c("#ff00ff","#0000ff"))
abline(v = t_mean,col="red",lwd=2)
text(t_mean,5.5,paste("Mean =",round(t_mean,3)),cex=2)

ブートストラップサンプルの平均値の分布はこちら。ピンクはブートストラップサンプル抽出に重複を許す場合、青は重複を許さない場合。母集団の平均を中心に分布していることが確認出来る。

f:id:nakhirot:20131216220554p:plain

3. Permutation test

以前こちらで扱った、「高齢ラットと若齢ラットの膀胱筋片における高濃度ノルアドレナリン刺激に伴う最大弛緩」を用いてPermutation Testを行う。

old <- c(20.8,2.8,50.0,33.3,29.4,38.9,29.4,52.6,14.3)
young <- c(45.5,55.0,60.7,61.5,61.1,65.5,42.9,37.5)
stat_t <- mean(old) - mean(young)

対応の無い2群(old, young)の平均値の差が有意であるかどうかを検定する。2群を結合したデータ(下記ではdat)の順番をランダムに並び替え、新たな2群(old_s,young_s)を生成する。各群の平均値を計算し、両者の差を計算する。これを1000回繰り返す。

k <- 1000
mat <- matrix(0,nrow=k,ncol=length(c(old,young)))
old_s <- mat[,1:length(old)]
young_s <- mat[,(length(old)+1):ncol(mat)]
mean_o <- numeric(k)
mean_y <- numeric(k)
stat <- numeric(k)

for (i in 1:k) {
  dat <- c(old,young)
  mat[i,] <- dat[order(runif(length(dat)))]
  old_s[i,] <- mat[i,1:length(old)]
  young_s[i,] <- mat[i,(length(old)+1):ncol(mat)]
  mean_o[i] <- mean(old_s[i,])
  mean_y[i] <- mean(young_s[i,])
  stat[i] <- mean_o[i] - mean_y[i]
}

1000回の試行で生成した平均値の差と、本来の平均値の差を比較する。

f:id:nakhirot:20131216230956p:plain

本来の2群の平均値の差(-23.546)を下回った試行は1000回中2回しかないことがわかる。確率にして2/1000=0.2%。したがって、2群の平均値の差は有意水準5%で有意であることがわかる。 

4. Cross Validation

以前、Decision Treeモデルを作成する際、精度評価にCross Validationを使用した。

f:id:nakhirot:20131216231505p:plain

f:id:nakhirot:20131216231516p:plain

f:id:nakhirot:20131216231530p:plain

上記の"x-validation Relative Error"は、Cross validationを行って算出したモデルの精度である。

f:id:nakhirot:20131216231632p:plain

コードはこちらを参照。