resampling2.utf8.md

# −−−−−−−−−−−−−−−−−
# 【データからの無作為サンプリング】
# −−−−−−−−−−−−−−−−−

data <- c(0,1,2,3,4,5,6,7,8,9)
  #10個の整数0〜9をデータとする

sample(data, 10, replace=T)

##  [1] 5 7 9 0 5 8 6 6 1 1

  # data から重複を許して10個を再抽出（ブーツストラップ）
  # 試行を反復するたびに数列にばらつきが生じる

sample(data, 9, replace=F)

## [1] 8 9 5 0 1 2 3 6 4

  # data から重複を許さず9個を再抽出（ジャックナイフ）
  # 試行を反復するたびに数列にばらつきが生じる


# −−−−−−−−−−−−−−−−−−
# 【正規分布からの無作為サンプリング】
# −−−−−−−−−−−−−−−−−−

normal.data<- rnorm(10000, mean=0, sd=1); 
hist(normal.data, freq=F);
x <- seq(-4,4,0.05);
curve(dnorm(x, mean=0, sd=1), add=T);

mean(normal.data)

## [1] -0.003350634

  # 標準正規分布 N(0,1) から10,000個の乱数を抽出し
  # ヒストグラムと N(0,1) を描画して標本平均を求める
  # 反復ごとに標本平均はばらつく．

multisample1 <- numeric(100);
for (i in 1:100) multisample1[i] <- 
  + mean(rnorm(100, mean=0, sd=1));
hist(multisample1, freq=F); var(multisample1); sd(multisample1)

## [1] 0.01011974

## [1] 0.1005969

  # 正規母集団から毎回100個の正規乱数を抽出し，
  # 標本平均を求めるという試行を100回反復し，
  # 標本平均のヒストグラムを描画し，分散と標準偏差を計算する．

multisample2 <- numeric(100);
for (i in 1:100) multisample2[i] <- 
  + mean(rnorm(1000, mean=0, sd=1));
hist(multisample2, freq=F); var(multisample2); sd(multisample2)

## [1] 0.001052704

## [1] 0.03244541

  # 正規母集団から毎回1,000個の正規乱数を抽出し，
  # 標本平均を求めるという試行を100回反復し，
  # 標本平均のヒストグラムを描画し，分散と標準偏差を計算する．

multisample3 <- numeric(100);
for (i in 1:100) multisample3[i] <- 
  + mean(rnorm(10000, mean=0, sd=1));
hist(multisample3, freq=F); var(multisample3); sd(multisample3)

## [1] 8.167801e-05

## [1] 0.009037589

  # 正規母集団から毎回10,000個の正規乱数を抽出し，
  # 標本平均を求めるという試行を100回反復し，
  # 標本平均のヒストグラムを描画し，分散と標準偏差を計算する．

multisample4 <- numeric(100);
for (i in 1:100) multisample4[i] <- 
  + mean(rnorm(100000, mean=0, sd=1));
hist(multisample4, freq=F); var(multisample4); sd(multisample4)

## [1] 9.354947e-06

## [1] 0.003058586

  # 正規母集団から毎回100,000個の正規乱数を抽出し，
  # 標本平均を求めるという試行を100回反復し，
  # 標本平均のヒストグラムを描画し，分散と標準偏差を計算する．

# −−−−−−−−−−−−−−−−−−−
# 【正規乱数データからのリサンプリング】
# −−−−−−−−−−−−−−−−−−−

normal.data<- rnorm(10000, mean=0, sd=1);
hist(normal.data, freq=F);

mean1 <- mean(normal.data); mean1;

## [1] -0.006097712

sd1 <- sd(normal.data); sd1

## [1] 0.9988278

  # 標準正規分布 N(0,1) から10,000個の乱数を抽出し
  # そのヒストグラムと平均ならびに標準偏差を求める

# 【ノンパラメトリック・ブーツストラップ】−−−

nonparaboot <- numeric(100);
for (i in 1:100) nonparaboot[i] <- 
     mean(sample(normal.data, 10000, replace=T));
hist(nonparaboot, freq=F);

var(nonparaboot)

## [1] 0.000107567

  # データから重複を許して10,000個をリサンプリングし
  # その試行を100回反復する．毎回，標本平均を求め，
  # そのヒストグラムを描画し，標本平均の分散を計算する．

nonparaboot2 <- numeric(100);
for (i in 1:100) nonparaboot2[i] <- 
     mean(sin(sample(normal.data, 10000, replace=T)));
hist(nonparaboot2, freq=F);

var(nonparaboot2)

## [1] 4.201535e-05

  # たとえ“異常”な統計量であっても，リサンプリングに
  # よって，その誤差を求めることができる．たとえば，デ
  # ータを正弦変換（sin）した変量の標本平均を求めた上で，
  # そのヒストグラムを描画し，標本平均の分散を計算する．

# 【ジャックナイフ（delete-half jackknife）】−−−

jackknife <- numeric(100);
for (i in 1:100) jackknife[i] <- 
     mean(sample(normal.data, 5000, replace=F));
hist(jackknife, freq=F);

var(jackknife)

## [1] 9.937438e-05

  # データから重複を許さず5,000個をリサンプリングし
  # その試行を100回反復する．毎回，標本平均を求め，
  # そのヒストグラムを描画し，標本平均の分散を計算する．

# 【パラメトリック・ブーツストラップ】−−−

paraboot <- numeric(100);
for (i in 1:100) paraboot[i] <- 
     mean(rnorm(10000, mean=mean1, sd=sd1));
hist(paraboot, freq=F);

var(paraboot)

## [1] 0.0001088622

  # データから推定された平均と標準偏差を正規分布パラメーターとする．
  # その正規分布から10,000個の正規乱数を100回抽出し，
  # 標本平均を求め，ヒストグラムを描画し，標本平均の分散を計算する．

resampling2.R

minaka

Mon Jul 9 13:42:32 2018