廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

R による溜池サンプリング(Reservoir Sampling)の実験

Data Stream Management という本に出てくる Reservoir Sampling(溜池サンプリング)という手法をシミュレーションしてみたい.

これはサイズ N の母集団(N は未知でもよい)からサイズ nn \le N)のランダム標本を非復元抽出で取ってくるアルゴリズムで, 大きすぎるデータからリサンプリングしたいときや, センサーなどから逐次的に流れてくるようなデータからランダム標本を得たいときに使えそう.

Data Stream Management: Processing High-Speed Data Streams (Data-Centric Systems and Applications)

Data Stream Management: Processing High-Speed Data Streams (Data-Centric Systems and Applications)

Data Stream Management よりも元論文(McLeod and
Bellhouse (1983); http://www.stats.uwo.ca/faculty/aim/vita/pdf/McLeod83.pdf)のほうがわかりやすかった.

アルゴリズム

1. k = 0 とする
2. 母集団にまだ要素が残っているならば k =k+1 として次のステップへ. 残っていなければ終了.
3 (a).  k \le n ならば母集団の k 番目の要素をサンプルの k 番目の要素としてステップ 2 へ.
3 (b). k > n ならば, 範囲 1,\ldots, k の離散一様乱数を取り j とする. もし j \le n ならば, 現在のサンプルの j 番目の要素を, 母集団の k 番目で置き換え, ステップ 2 へ.

なぜこれでうまくいくのかは, 大量のテキストからランダムに少数の行を抽出したい - Reservoir Sampling - 唯物是真 @Scaled_Wurm による説明がわかりやすい.

R による実験

母集団としてサイズ N=100,000 の混合正規乱数を使う.

ここから Reservoir Sampling によってサイズ n=5,000 の標本を取ってきた.

次の図の上が母集団のヒストグラムで下が標本のヒストグラム. 元の分布の形状を保っていることがわかる.

f:id:abrahamcow:20161117011107p:plain

次の図は母集団と標本の要素の順番ごとの, 要素の値の折れ線グラフである. 上が母集団で下が標本. もとは要素の順番に依存する傾向があったが, Reservoir Sampling された系列は定常になっていることがわかる.

f:id:abrahamcow:20161117011239p:plain

以下 R のコードです.

set.seed(1)
rmixnorm3 <- function(n) {
  n1 <- round(n*0.5)
  n2 <- round(n*0.2)
  n3 <-n-(n1+n2)
  c( rnorm(n1,-3), rnorm(n2,3),rnorm(n3,9) )
}
N <-100000
n <-5000
allX <-rmixnorm3(N) #母集団
###
#Reservoir Sampling
e<-allX[1:n] #最初は n 番目までを標本にする
for(i in (n+1):N){
  newX <- allX[i]
  u <-sample.int(i,1) #uniform on {l,2,...,i}
  if(u < n){
    e[u] <- newX
  }
}
###
#以下プロット
par(mfrow=c(2,1))
hist(allX,breaks = "FD")
hist(e,breaks = "FD")
plot(allX,  type = "l")
plot(e,  type = "l")

abrahamcow.hatenablog.com