Upgrade to Pro — share decks privately, control downloads, hide ads and more …

Bootstrap法とその周辺

Daisuke Yoneoka
November 14, 2023

 Bootstrap法とその周辺

Daisuke Yoneoka

November 14, 2023
Tweet

More Decks by Daisuke Yoneoka

Other Decks in Research

Transcript

  1. 複雑な問題 • 統計量の正確な標本分布の導出は困難 • 漸近理論の発展 o 最尤法 • 一致性、漸近正規性 •

    漸近分散:Fisher情報量の逆数 • 統計量の関数の分布に関してはデルタ法 • 以下の様な統計量においては結構難しい o トリム平均、Median → 分位点に関する理論の発展が必要
  2. 抽出方法 • 復元抽出 sampling with replacement o 同じ要素の抽出を許す o Rコード:

    sample()関数 • Ex. sample(data, 100, replace =True) • 非復元抽出 sampling without replacement o 同じ要素の抽出を許さない o Rコード: sample()関数 • Ex. sample(data, 100, replace =False)
  3. ごちゃごちゃしたこたぁ いいんだよ! men.h <- c(26.6,37.2,37.9,36.6,35.6,37.1,40.1,37.4,37.8,36.6) mean.boot <- numeric(2000) set.seed(314) for

    (b in 1:2000){ i <- sample(1:10, replace=TRUE) # 1から10までの整数を10回無作為に抽出 men.boot <- men.h[i] # b回目のブートストラップ標本 mean.boot[b] <- mean(men.boot) # b回目のブートストラップ標本平均 } hist(mean.boot, freq=F, xlab="bootstrap mean", main="") # 平均のヒストグラム # 両側95%信頼区間 sort(mean.boot)[c(0.025*2000,0.975*2000)] CPPUTUSBQNFBO %FOTJUZ         
  4. Cramér-Raoの下限と効率 • Fisher情報量(I(θ))の復習 o f(Y;θ)は尤度関数 o 意味:1回微分の分散か2回微分の期待値 • いくつかの正則条件のもとで o

    クラメール・ラオの下限 (ただし、 ) • の分散がクラメール・ラオの下限を達成するときに、 を有効推定量(Efficient estimator) o 最尤推定量が有効推定量であることは稀。通常は漸近的に達成(漸近有効)
  5. ブートストラップ誤差 • 統計的誤差 o 差込原理より として近似したことからくる誤差 o どうしようもないから諦めよう!(提案) o でも、nは大きくしようね!

    • モンテカルロ誤差 o シミュレーションに基づく誤差 o 何回反復させるかに依存しているので、十分回数やろう! o で、結局何回くらいが適当なの? • nが大きい場合、反復回数を増やす • 中央値のような標本の滑らかなでない関数の場合反復回数を増やす o Efron and Tibshirani(1993) によると、分散や標準誤差のブートストラップ推定 の場合は25-300回程度十分らしい!
  6. Jackknife法 • もう一つのリサンプリング法 o 重複を許さないリサンプリング法 o 狭義にはこんなかんじで1つだけサンプルを抜いてリサンプリング o イメージ的にはCross validationによく似ている。

    • どうでもいいけど、語源は「キャンプ場ですげー便利」 • 利点 o Bootstrapよりちょっと早い • 欠点 o 統計量が平滑でない値の場合、失敗する場合がある。(ex. Median) o 平滑性=データの変化がどれくらい統計量を変化させるか
  7. Jackknifeの失敗 x<-sample(1:100, size=10) #標準誤差のジャックナイフ推定量 M<-numeric(10) for (i in 1:10){ y<-x[-i]

    M[i]<-median(y) } Mbar<-mean(y) print(sqrt((10-1)/10*sum((M-Mbar)^2))) [1] 38.54363 #標準誤差のブートストラップ推定量 Mb<-replicate(1000, expr={ y<-sample(x, size=10, replace=T) median(y)} ) print(sd(Mb)) [1] 11.94611 Jackknifeのmedianの標準誤差とBootstrap のmedianの標準誤差が大きく違う 何かおかしい! Jackknifeが推定誤差を起こす!
  8. Bootstrap信頼区間 • 標準正規Bootstrap CI • 基本Bootstrap CI • Percentile Bootstrap

    CI • Bootstrap T CI • BCa法 (Bias corrected and accelerated method) o 性能や特性など詳しくは、A.C. Davison et al(1997)
  9. 標準正規/基本 B-CI • 標準正規CI o まぁ想像通りです。 o 仮定が強い • の分布が正規分布

    or • が標本平均 and サンプルサイズが大きい(中心極限定理) • 基本CI o Bootstrap近似に基づく。 • と近似(この分位点を計算)→誤差が大きいかも o ただし、 の経験累積分布から標本のα分位点
  10. Bootstrap T CI • 基本Bootstrap CIの場合、 としているので、 分布のずれがある場合うまく行かない! o 一次の正確度しかないから

    • 一次の正確度: • Cは被覆誤差 • C→0 (n→∞)がであってほしい • それじゃ、二次のモーメント(分散)まで考えてみればいい じゃない!”t型”統計量の標本分布をリサンプリングで作成 • 信頼区間 o は、 のα/2番目に小さい値 上側信頼限界
  11. Bootstrap T CI • 信頼区間 o は、 のα/2番目に小さい値 • 長所

    o 二次の正確性を持つ: • 短所 o σの推定が不可欠→ブートストラップ標本ごとにσを計算しなけれならないので、 計算負荷が大きい(つまり、ブートストラップのなかにブートストラップの入 れ子構造)
  12. Percentile法の正確度 • Percentile CIの方が標準正規CIより良い被覆率 • 変換後に左右対称となる の単調増加関数 が存在す るか否かに正確度が依存する。 o

    多くの場合、そんな なんて存在しないよ! • Efron and Tibshirani;1993 • 汪、桜井;2011 • ちょっと改善しましょ!→ BCa法
  13. で、CI求めるのってどれがいいの & 何回反復すりゃいいの? • うーん。Bootstrap-TかBCaかな? Byung-Jin Ahn et al; 2009

    • CIの計算には分散の計算時よりも大きい反復回数が必要 o 90-95% CIの場合は反復回数1000-2000回は必要だよ! • Efron and Tibshirani;1993
  14. 回帰分析に応用 • こんなかんじのコンプライアンスとコレステロール値の散布図と3次の回帰直線 # データセットづくり library(bootstrap) # z の値を大きさの順でデータを並べ替える zz

    <- sort(cholost$z) yy <- cholost$y[order(cholost$z)] # 拡大データフレームを作る mydata <- data.frame(z1=zz, z2=zz^2, z3=zz^3, yy) # データの散布図 plot(mydata$z1, mydata$yy, xlab="compliance", ylab="decrease in cholesterol level") # 最小2乗法による3次関数のあてはめ cubic <- lm(yy~., data=mydata) # 推定された回帰曲線を描く lines(zz, predict(cubic), lty=2)              DPNQMJBODF EFDSFBTFJODIPMFTUFSPMMFWFM
  15. # データの散布図 plot(mydata$z1, mydata$yy, xlab=“compliance”, ylab=“decrease in cholesterol level”) set.seed(314159)

    # 乱数の種を固定する B <- 100 # ブートストラップ反復回数 n <- length(cholost$z) # 標本の大きさ r60 <- numeric(B) # ブートストラップ回帰平均値 r100 <- numeric(B) # ブートストラップ回帰平均値 for (b in 1:100){ # ブートストラップ反復開始 bt <- sample(1:n, replace=TRUE) # ブートストラップ標本番号 mydata <- cholost[bt,] # ブートストラップ標本 zz <- sort(mydata$z) # zの値の並べ替え yy <- mydata$y[order(mydata$z)] # zの大きさの順にyを並べ替える mydata <- data.frame(z1=zz, z2=zz^2, z3=zz^3, yy) # データフレームを作る cubic <- lm(yy~., data=mydata) # 最小2乗法による3次関数のあてはめ lines(zz, predict(cubic)) # 求めた最小2乗曲線を描く # z=60, 100 のときのブートストラップ回帰予測値 dumy <- predict(cubic, data.frame(z1=c(60, 100), z2=c(60, 100)^2, z3=c(60, 100)^3)) r60[b] <- dumy[1] r100[b] <- dumy[2] }              DPNQMJBODF EFDSFBTFJODIPMFTUFSPMMFWFM