読者です 読者をやめる 読者になる 読者になる

Rを通じて統計学を学ぶ備忘録ブログ

SPSSからRに移行したい私のような人向けのR解説ブログ兼学習用備忘録。

第78回 ブートストラップ法

今回はブートストラップ法について見ていきましょう。ブートストラップ法はモンテカルロ法の一種で、標本から標本を再抽出することで、母集団の性質を推測する方法のことです。観測されたサンプルデータから母集団の性質を推測するとき、必ず誤差が生じています。その誤差を評価するために、1つの標本から復元抽出を繰り返すことで、たくさんのサンプルデータを作り、そのサンプルデータから母集団の性質やモデルの推測の誤差などを評価する方法がブートストラップ法です。

■ブートストラップ法~平均と信頼区間

平均値をブートストラップ法を使って評価していきましょう。その後、ブートストラップ分布から信頼区間を求めます。以下のような手順になります。

1.標本から無作為抽出を行い、1番目の標本を作り、1番目の平均値を求めます。

2.手順1をn回繰り返します。

3.n個の標本をもとにブートストラップ分布を作成し、信頼区間を求めます

f:id:teruaki-sugiura:20160313155420p:plain

------------------------手順1~2のイメージ----------------------
# 母集団からのサンプリングA
sample <- rnorm(10, 0, 1)
mean(sample)
var(sample)
# サンプリングAから1回目のリサンプリング
x <- sample(1:10, replace = TRUE)
x1 <- sample[x]
mean(x1)
var(x1)
# サンプリングAから2回目のリサンプリング
x <- sample(1:10, replace = TRUE)
x2 <- sample[x]
mean(x2)
var(x2)
----------------------------------------------------------------------------------  

サンプルからリサンプリングして、サンプルを沢山作って、サンプルの母集団の性質を評価することがブートストラップ法ですが、とりあえず、百聞は一見に如かずということなのでRで実践していきます。

標準正規分布からまずサンプルサイズ100のデータ”sample”を生成します。次にその"sample"からリサンプリングを実施し、平均値を求めます。それをmean.bootstrap[1](="sample"からリサンプリングされたサンプルサイズ100の平均値で、1つ目のブートストラップ標本)とします。

これをmean.bootstrap[1],mean.bootstrap[2],...,mean.bootstrap[100]まで繰り返して、そのヒストグラムを書きます。このヒストグラムのことをブートストラップ分布と呼びます。

--------------100サンプルサイズ平均値のリサンプリングを100回行う------
#イメージとしてmean.bootstrap[1],...,mean.bootstrap[100]を作る
sample <- rnorm(100, 0, 1)
mean.bootstrap <- numeric(100)
for(i in 1:100){
        t <- sample(1:100, replace = TRUE)
        x <- sample[t]
        mean.bootstrap[i]<- mean(x)
      }
head(round(mean.bootstrap, 3), 10)
# 95%信頼区間
quantile(mean.bootstrap, p=c(0.025,0.975))
hist(mean.bootstrap, freq = FALSE, xlab = "mean.bootstrap", main = "bootstrap method")
-------------------同様のことは以下のコードでもできます------------
library(simpleboot)
library(boot)
#mean.bootstrap[1],...,mean.bootstrap[100]を作る
mean.bootstrap <- one.boot(sample, mean, 100)  
#mean.bootstrap[1],...,mean.bootstrap[100]の中身を確認
mean.bootstrap$t
hist(mean.bootstrap, freq = FALSE, xlab = "mean.bootstrap", main = "bootstrap method")
boot.ci(mean.bootstrap, conf = 0.99, type = "all") 
f:id:teruaki-sugiura:20160313202544p:plain
---------------------------------------------------------------------------- 

 

boot.ci()の結果”Normal”は通常の区間推定、”Percentile”は推定値を小→大の順に並べた列の100α%点。区間推定値の精度を高めるには、ブートストラップの標本数を増やせば解決できます。

■ブートストラップ法で相関係数

観測されたデータx, yをもとに相関係数を求めて、それは真の相関係数の1つの推定値でしかありません。この値こそ、真の値です!!と主張しても、相手は納得し難いと思います。そこで、ある程度の散らばりがわかれば、区間を持って相関係数の推測について主張できます。そんな時、ブートストラップ法を使って、相関係数区間推定を行います。

1. 標本からx, yを無作為抽出を行い、1番目の標本組を作り、1番目の相関係数を求めます。

2. 手順1をn回繰り返し、n組作ります。

3. n組の標本をもとにブートストラップ分布を作成し、信頼区間を求めます

-----------リサンプリングを100回行い、相関係数と信頼区間を求める----- 
x <- rpois(100, 3) 
y <- rnorm(100, 0) 
bootstrap.cor <- pairs.boot(x, y, cor, 100)
boot.ci(bootstrap.cor, conf = 0.95, type = "all")
f:id:teruaki-sugiura:20160313203416p:plain
hist(bootstrap.cor)
f:id:teruaki-sugiura:20160313202415p:plain
----------------------------------------------------------------------------
a <- c(71, 68, 66, 67, 70, 71, 70, 73, 72, 65)
b <- c(69, 64, 65, 63, 65, 62, 65, 64, 66, 59)
bootstrap.cor <- numeric(200) # 200回分の格納用 
for ( i in 1:200 ) {
          t <- sample(1:10, 10, replace = TRUE)
          x <- a[t]
          y <- b[t]
          bootstrap.cor[i] <- cor(x, y)
        } 
hist(n)
----------------------------------------------------------------------------  

■ブートストラップ法の仮説検定~ブートストラップ検定~

2標本問題の統計検定に対するブートストラップ検定を見ていきます。なんらかの分布に従う標本xとyが得られた時に、ブートストラップ法を利用し平均値の差を検討する問題を考えます。以下の手順で検定で行っていきます。

1. 新しい変数zをxとyを混合することで作ります。この時、x, y, zに対する確率変数をX, Y ,Z、検定統計量をt(Z)=t(X, Y)とすると、p値は以下の式で求められます。

f:id:teruaki-sugiura:20160313190112p:plain

t(z)は検定統計量t(Z)の実現値です。z=(x, y)からサンプルサイズm, nの標本x^*b, y^*bを無作為復元抽出し、これをz^*b={x^*b, y^*b}とします。

2. 手順1をB回繰り返し、ブートストラップ検定統計量を計算します。そして、ブートストラップp値のモンテカルロ近似により計算します。

f:id:teruaki-sugiura:20160313195005p:plain

3. 得られたp値をもとに帰無仮説を棄却するかどうかを決定する。

----------------------------------------------------------------------------
x <- c(90, 200, 20, 40, 100, 150, 200, 300, 10)
y <- c(50, 10, 10, 10, 50, 30, 40, 50, 10, 30)
z <- c(x, y)
Z <- numeric(1000)
for(i in 1:1000){
  X <- sample(z, 9, replace = TRUE)
  Y <- sample(z, 10, replace = TRUE)
  Z[i] <- mean(X) - mean(Y)
  }
sum(Z >= mean(x) - mean(y))/1000
table(ifelse(Z >= mean(x) - mean(y), 1, 0))
hist(Z, main = "bootstarp test", xlab = "θ*", freq = TRUE)
abline(v = mean(x) - mean(y), lty = 2) 
f:id:teruaki-sugiura:20160313202640p:plain
-------------------------モンテカルロ検定------------------------------
library(boot)
fir <- fir$count
mean(fir)
var(fir)
N <- 1999
stat.s <- numeric(N)
s <- sum(fir)
stat <- function(d){ 49*var(d)/mean(d)}
y.s <- rmultinom(B, size = s, prob = rep(1/50, 50))
for(i in 1:B){
       stat.s[i] <- stat(y.s[,i])
     }
(1 + sum(stat.s >= stat(fir)))/(N + 1)
hist(stat.s, freq = FALSE)
abline(v = mean(stat.s), lty = 2)
---------------------------------多項分布の補足------------------------------
6面のサイコロで各出目の確率は1/6とする。
①rmultinom(n=1,size=1,prob=c(1/6,1/6,1/6,1/6,1/6,1/6))
これを1回実行すると、さいころを1回振る確率シミュレーションを1セット行う。 
②rmultinom(n=1,size=10,prob=c(1/6,1/6,1/6,1/6,1/6,1/6))
この場合、10 回振って出目の頻度を得る確率シミュレーションを1セット行う。
③rmultinom(n=10,size=5,prob=c(1/6,1/6,1/6,1/6,1/6,1/6))
これだと、5回振って出た目を得る確率シミュレーションを10セット行う。
----------------------------------------------------------------------------

 

■ブートストラップ法で回帰分析

例として、p個の説明変数(x_1, x_2,...,x_p)を持つ線形回帰モデルを考えます。つまり、x=(x_1, x_2,...,x_p)^TでE(Y|x)=β^T*xというモデルです。通常、このようなモデルでは、誤差項eには平均0,分散σ^2という制約がつきます。この条件があるからこそ、区間推定や仮説検定の対して厳密な理論を展開することができるのですが、誤差項が正規分布していないような場合に、ブートストラップは役に立ちます。ブートストラップの手順は以下の通りです。

※簡略したイメージです。

1. 中心化調整済み残差r_i-r^-で構成される経験分布関数から復元抽出し、ブートストラップ残差ε*_1,ε*_2,...,ε*_nを得る。

2. ブートストラップ標本Y*_i=β^_0+β^_1*x_i+ε*_iを計算する。

3. (x_1, Y*_1), (x_2, Y*_2),...,(x_n, Y*_n)からβ^_0*とβ^_1*、(σ^*)^2を計算する

ブートストラップ法では、誤差の従うなんらかの分布を、互いに独立な残差(e_i=y_i-y^_i)を利用することで、分布を検討します。なんらかの分布は適当な条件のもとで分布収束します。注意として、残差e_iは互いに独立で、iごとに異なるなんらかの分布に従う可能性があります。つまり、E(e_i)=0にはなるが、var(e_i)≠σ^2のようなことが起こりうります。

このような問題に対処するために以下のように標準化された残差の分散を考えます。

f:id:teruaki-sugiura:20160313154959p:plain

これは、nが大きくなるとσ^2に近似できます。この式を用いることで、以下のような平均0、分散がiに依存せずσ^2となる、修正済み残差r_iとなります。

f:id:teruaki-sugiura:20160313155009p:plain

ブートストラップ法で計算されるβは以下のように計算されます。

f:id:teruaki-sugiura:20160313155018p:plain

以上の手順で、信頼区間を改良すると、n回ブートストラップを適用したとすると、以下のようになります。

f:id:teruaki-sugiura:20160313155029p:plain

-------残差からのリサンプリングブートストラップ回帰分析-----------------
y <- c(21, 16, 23, 20, 22, 17, 17, 25, 23, 22, 21, 23, 17, 23, 17, 20, 24, 22, 26, 17)
x <- c(15, 23, 18, 17, 50, 19, 24, 25, 75, 11, 13, 35, 25, 75, 22, 18, 60, 12, 20, 21)
# least squares method of β1
β1 <- function(x, y){                
          Sxx <- sum(x - mean(x))*(x - mean(x))
          Sxy <- sum(x - mean(x))*(y - mean(y))
          Sxy/Sxx
        }
# least squares method of β0 
β0 <- function(x, y){                 
          mean(y) - β1(x, y)*mean(x)
        }
# 中心化修正済み残差
r <- function(x, y){                     
          n <- length(x)
          Sxx <- sum(x - mean(x)*(x - mean(x)))
          mu <- β0(x, y) + β1(x, y)*x
          rr <- (y - mu)/sqrt(1 - 1/n - (x - mean(x))^2/Sxx)
          rr - mean(rr)
        }
# 標準偏差
σ <- function(x, y){                 
          n <- length(x)
          mu <- β0(x, y) + β1(x, y)*x
          sqrt(sum(y - mu)*(y - mu)/(n - 2))
        }
# ブートストラップ開始
B <- 2000                                       
t.s <- numeric(B)                                                       
residual <- r(x, y)                    
Sxx <- sum(x - mean(x))^2)
# 1000回の反復
for (b in 1:B){                                             
        e.s <- sample(residual, replace=TRUE)    # 残差のリサンプリング
        y.s <- β0(x, y) + β1(x, y)*x + e.s      # y* の値
        t.s[b] <- β1(x, y.s) - β1(x, y) 
        t.s[b] <- t.s[b]/σ(x, y.s)*sqrt(Sxx)       # t* の値
      }
# t*のブートストラップ分布
hist(t.s, main="Bootstrap dist of t*", xlab="t*")
f:id:teruaki-sugiura:20160313202947p:plain
# β1の95%信頼上限下限区間
lower <- sort(t.s)[c(0.975)*B]               
lower <- β1(x, y) - lower*σ(x, y)/sqrt(Sxx)
upper <- sort(t.s)[c(0.025)*B]               
upper <- β1(x, y) - upper*σ(x, y)/sqrt(Sxx)
----------------------------------------------------------------------------------
B <- 2000       
t.s <- numeric(B)                                                       
residual <- r(x, y)                    
Sxx <- sum(x - mean(x))^2)
# 1000回の反復
for (b in 1:B){                                 
        e.s <- sample(residual, replace=TRUE)         # 残差のリサンプリング
        y.s <- β0(x, y) + β1(x, y)*x + e.s       # y* の値 
        t.s[b] <- β0(x, y.s) - β0(x, y)
        t.s[b] <- t.s[b]/σ(x, y.s)*sqrt(Sxx)       # t* の値
        }
# t*のブートストラップ分布 
hist(t.s, main="Bootstrap dist of t*", xlab="t*") 
f:id:teruaki-sugiura:20160313202955p:plain
# β0の95%信頼上限下限区間
lower <- sort(t.s)[c(0.975)*B]               
lower <- β0(x, y) - lower*σ(x, y)/sqrt(Sxx)
upper <- sort(t.s)[c(0.025)*B]               
upper <- β0(x, y) - upper*σ(x, y)/sqrt(Sxx) 
--------------------------単回帰分析&重回帰分析----------------------------
installed.packages("simpleboot")
library(simpleboot) 
lmboot.model <- lm(y ~ x)
lmboot.1000 <- lm.boot(lmboot.model, R = 1000)
summary(lmboot.1000)
plot(lmboot.1000, xlab = "x", ylab = "y")
f:id:teruaki-sugiura:20160313203035p:plain
lmboot.model2 <- lm(Ozone ~ Wind + Solar.R、data = airquality)
lmboot.2000 <- lm.boot(lmboot.model2, R = 1000)
summary(lmboot.2000)
 
-----サンプルデータからのリサンプリングブートストラップ回帰分析------
y <- c(21, 16, 23, 20, 22, 17, 17, 25, 23, 22, 21, 23, 17, 23, 17, 20, 24, 22, 26, 17)
x <- c(15, 23, 18, 17, 50, 19, 24, 25, 75, 11, 13, 35, 25, 75, 22, 18, 60, 12, 20, 21)
# least squares method of β1
β1 <- function(x, y){                 
            Sxx <- sum(x - mean(x))*(x - mean(x))) 
            Sxy <- sum(x - mean(x))*(y - mean(y)))
            Sxy/Sxx
        }
# least squares method of β0
β0 <- function(x, y){                
          mean(y) - β1(x, y)*mean(x)
          }
# ブートストラップ開始
B <- 5000                                
β0.s <- numeric(B)                  # β0*を入れる箱
β1.s <- numeric(B)                  # β1*を入れる箱
n <- length(x)                      # サンプルサイズ
# ブートストラップ開始
for (b in 1:B){                             
          bt <- sample(1:n, replace=TRUE)   # 残差のリサンプリング
          β0.s[b] <- β0(x[bt], y[bt])      # β0*の値
          β1.s[b] <- β1(x[bt], y[bt])      # β1*の値
        }    
# β1*, β0*のブートストラップ分布
hist(β1.s, main = "Bootstrap dist of β1*", xlab = "β1*")  
hist(β0.s, main = "Bootstrap dist of β0*", xlab = "β0*")  
f:id:teruaki-sugiura:20160313203140p:plain
f:id:teruaki-sugiura:20160313203142p:plain
# データからのリサンプリングに基づくβ1, β0の信頼区間
sort(β1.s)[c(B*0.025, B*0.975)]
sort(β0.s)[c(B*0.025, B*0.975)]
f:id:teruaki-sugiura:20160313203236p:plain 
------------------------------------------------------------------------------------- 

以上でおしまい。

広告を非表示にする