「ノーマル・チップスを用いたシミュレーションをRで再現する」の編集履歴（バックアップ）一覧はこちら
追加された行は緑色になります。
削除された行は赤色になります。
#contents()
*ノーマル・チップスの作成
テキストp72、「シューハートのノーマル・チップス」をベクトルとして作成する。
norm.tips <- rep(0:60, c( 1, 1, 1, 1, 1, 2, 2, 3, 4, 4,
5, 7, 8, 9,11,13,15,17,19,22,
24,27,29,31,33,35,37,38,39,40,
40,40,39,38,37,35,33,31,29,27,
24,22,19,17,15,13,11, 9, 8, 7,
5, 4, 4, 3, 2, 2, 1, 1, 1, 1, 1))
これで、norm.tips というオブジェクトの中に998 枚のノーマル・チップスが入ったことになる。rep 関数は第一引数に与えたベクトルを第二引数に与えた回数だけつなげたベクトルを作成する関数で、第一引数に0～60まで1ずつ増えるベクトル(＝カードの番号) を指定し、第二引数でそれぞれの値を何回繰り返すか(＝カードの枚数) を指定している。
ノーマル・チップスの平均と標準偏差を確認する。平均の計算はmean 関数、標準偏差の計算はsd 関数で行う。
> mean(norm.tips)
[1] 30
> sd(norm.tips)
[1] 9.959396
つぎにヒストグラムを描く。ヒストグラムはhist 関数で描けるが、見た目がいいのでturehist 関数を使う。truehist 関数を使う前にはlibrary 関数でMASS パッケージを読み込んでおく必要がある。
library(MASS)
truehist(norm.tips)
テキストにあるとおり、「近似的に、平均値が30 で、標準偏差が10 の正規分布(p73)」となっていることがわかる。
*サンプリング
ノーマル・チップスからのサンプリングはsample 関数により行う。
sample(norm.tips, 5, replace=T)
ここで1 つ目の引数はサンプリングするオブジェクトを指定する。ここではnorm.tips から抽出するのでnorm.tips を指定している。2 つ目の引数はサンプリングする個数で、ここでは5 個(カード5 枚) の抽出を指定している。3 つめの引数は重複サンプリングを許すか否かの指定で、T(＝ TRUE) は許すという指定である。テキスト中の「取り出した一枚は再び箱の中に戻し(p74)」という指定に従った。
*チップ実験1、$$\bar{x}$$の分布
p74 からのチップ実験は次のようなものである。
+チップを1 枚引き、番号を記録してつぼに戻す
+1 の作業を5 回繰り返し、番号を5 つ記録する
+5 つの番号の平均値を計算し、記録する
+1～3 の作業を100 回繰り返し、平均値を100 組分記録する
+平均値¯x の平均、標準偏差はいくらか？また分布の形はどうなるか？
これをR 上のノーマル・チップスで行うには、sample 関数を100 回繰り返せばよい(あるいはもっと多く、例えば1 万回にしてもさほどの時間はかからない)。関数を繰り返し実行させ、各実行結果を記録するにはfor構文とを使う。
x.mean <- numeric(100)
for(i in 1:100){
x.mean[i] <- mean(sample(norm.tips, 5, replace=T))
}
for 構文前のnumeric(100) 関数は平均値を記録するためのオブジェクト(入れ物) を作成している。100 とすれば100 個の入れ物になるが、別にnumeric(0) などとしてもR は自動で入れ物を広げてくれる。ただしオブジェクトの拡張には少しだけ時間がかかるので、必要なサイズがわかっている場合はそれにあわせて指定すべきである。
for 構文は(i in M) という部分と式という2 つの部分からなる。(i in M) の部分でまずi は繰り替えし処理の間に変化させる引数である。i という記号である必要は無く、j でもk でも何でもいい。1:100 の部分に指定したベクトルが最初から順に引数に代入されていく。つまり、1:100 と指定した今回の場合なら、「i に1 を代入して式を実行」「i に2 を代入して式を実行」「i に3 を代入し
て式を実行」…「i に99 を代入して式を実行」「i に100 を代入して式を実行」といった具合に繰り返し処理が実行される。
今回、i という記号は添え字に利用している。x.mean[i] はx.mean オブジェクトのi 番目を参照する。つまり、i 番目に計算されたサンプルの平均値がx.mean のi 番目に代入される。あとは、得られた100 個の¯x について平均や標準偏差を調べ、ヒストグラムを描けば目的を達する。
**問題 1
テキストp89 からの分散比に関するチップ実験を行え。
-分散の計算はvar 関数によって行う。
-比を計算するので、サンプリングは1 回の繰り返し中に2 回必要である。
**問題 2
テキストp100 からの「玉の抽出実験」を行え。
-黒を0、白を1 などとして考える。サンプリング結果の和が白玉の個数となる。
-和はsum 関数により計算する。
**問題 3
p79 の32 図、「試料の数$$n$$を増したときの$$\bar{x}$$の制度」をチップ実験により再現し、「サンプル数が$$n$$の時の$$\bar{x}$$の標準偏差」をプロットせよ。また、同様に「サンプル数が$$n$$のときの$$\bar{x}$$ の分散」にサンプル数$$n$$を
乗じたものをプロットせよ。結果は何を意味するか？
-チップ実験1 で使ったfor 構文を100 回繰り返せばよい。
-100 回繰り返す間に、サンプル数を1～100 まで1 ずつ増やす。
-記録するのは「$$\bar{x}$$」ではなく、「サンプル数が$$n$$のときの$$\bar{x}$$の標準偏差」である。
-プロットはplot関数により行う。plot(x, y) とすれば、xをx軸に、yをy軸にして散布図が描かれる。ここではチップ実験で得られた標準偏差をy 軸に、その標準偏差が得られたときのサンプル数をx軸にする。
----
#contents()
*ノーマル・チップスの作成
テキストp72、「シューハートのノーマル・チップス」をベクトルとして作成する。
norm.tips <- rep(0:60, c( 1, 1, 1, 1, 1, 2, 2, 3, 4, 4,
5, 7, 8, 9,11,13,15,17,19,22,
24,27,29,31,33,35,37,38,39,40,
40,40,39,38,37,35,33,31,29,27,
24,22,19,17,15,13,11, 9, 8, 7,
5, 4, 4, 3, 2, 2, 1, 1, 1, 1, 1))
これで、norm.tips というオブジェクトの中に998 枚のノーマル・チップスが入ったことになる。rep 関数は第一引数に与えたベクトルを第二引数に与えた回数だけつなげたベクトルを作成する関数で、第一引数に0～60まで1ずつ増えるベクトル(＝カードの番号) を指定し、第二引数でそれぞれの値を何回繰り返すか(＝カードの枚数) を指定している。
ノーマル・チップスの平均と標準偏差を確認する。平均の計算はmean 関数、標準偏差の計算はsd 関数で行う。
> mean(norm.tips)
[1] 30
> sd(norm.tips)
[1] 9.959396
つぎにヒストグラムを描く。ヒストグラムはhist 関数で描けるが、見た目がいいのでturehist 関数を使う。truehist 関数を使う前にはlibrary 関数でMASS パッケージを読み込んでおく必要がある。
library(MASS)
truehist(norm.tips)
テキストにあるとおり、「近似的に、平均値が30 で、標準偏差が10 の正規分布(p73)」となっていることがわかる。
*サンプリング
ノーマル・チップスからのサンプリングはsample 関数により行う。
sample(norm.tips, 5, replace=T)
ここで1 つ目の引数はサンプリングするオブジェクトを指定する。ここではnorm.tips から抽出するのでnorm.tips を指定している。2 つ目の引数はサンプリングする個数で、ここでは5 個(カード5 枚) の抽出を指定している。3 つめの引数は重複サンプリングを許すか否かの指定で、T(＝ TRUE) は許すという指定である。テキスト中の「取り出した一枚は再び箱の中に戻し(p74)」という指定に従った。
*チップ実験1、$$\bar{x}$$の分布
p74 からのチップ実験は次のようなものである。
+チップを1 枚引き、番号を記録してつぼに戻す
+1 の作業を5 回繰り返し、番号を5 つ記録する
+5 つの番号の平均値を計算し、記録する
+1～3 の作業を100 回繰り返し、平均値を100 組分記録する
+平均値¯x の平均、標準偏差はいくらか？また分布の形はどうなるか？
これをR 上のノーマル・チップスで行うには、sample 関数を100 回繰り返せばよい(あるいはもっと多く、例えば1 万回にしてもさほどの時間はかからない)。関数を繰り返し実行させ、各実行結果を記録するにはfor構文とを使う。
x.mean <- numeric(100)
for(i in 1:100){
x.mean[i] <- mean(sample(norm.tips, 5, replace=T))
}
for 構文前のnumeric(100) 関数は平均値を記録するためのオブジェクト(入れ物) を作成している。100 とすれば100 個の入れ物になるが、別にnumeric(0) などとしてもR は自動で入れ物を広げてくれる。ただしオブジェクトの拡張には少しだけ時間がかかるので、必要なサイズがわかっている場合はそれにあわせて指定すべきである。
for 構文は(i in M) という部分と式という2 つの部分からなる。(i in M) の部分でまずi は繰り替えし処理の間に変化させる引数である。i という記号である必要は無く、j でもk でも何でもいい。1:100 の部分に指定したベクトルが最初から順に引数に代入されていく。つまり、1:100 と指定した今回の場合なら、「i に1 を代入して式を実行」「i に2 を代入して式を実行」「i に3 を代入し
て式を実行」…「i に99 を代入して式を実行」「i に100 を代入して式を実行」といった具合に繰り返し処理が実行される。
今回、i という記号は添え字に利用している。x.mean[i] はx.mean オブジェクトのi 番目を参照する。つまり、i 番目に計算されたサンプルの平均値がx.mean のi 番目に代入される。あとは、得られた100 個の¯x について平均や標準偏差を調べ、ヒストグラムを描けば目的を達する。
*練習問題
**問題 1
テキストp89 からの分散比に関するチップ実験を行え。
-分散の計算はvar 関数によって行う。
-比を計算するので、サンプリングは1 回の繰り返し中に2 回必要である。
**問題 2
テキストp100 からの「玉の抽出実験」を行え。
-黒を0、白を1 などとして考える。サンプリング結果の和が白玉の個数となる。
-和はsum 関数により計算する。
**問題 3
p79 の32 図、「試料の数$$n$$を増したときの$$\bar{x}$$の制度」をチップ実験により再現し、「サンプル数が$$n$$の時の$$\bar{x}$$の標準偏差」をプロットせよ。また、同様に「サンプル数が$$n$$のときの$$\bar{x}$$ の分散」にサンプル数$$n$$を
乗じたものをプロットせよ。結果は何を意味するか？
-チップ実験1 で使ったfor 構文を100 回繰り返せばよい。
-100 回繰り返す間に、サンプル数を1～100 まで1 ずつ増やす。
-記録するのは「$$\bar{x}$$」ではなく、「サンプル数が$$n$$のときの$$\bar{x}$$の標準偏差」である。
-プロットはplot関数により行う。plot(x, y) とすれば、xをx軸に、yをy軸にして散布図が描かれる。ここではチップ実験で得られた標準偏差をy 軸に、その標準偏差が得られたときのサンプル数をx軸にする。
----
*問題の解答
※解答は例です。もっと効率的でスマートな方法があるやもしれません。
**問題 1
※ノーマル・チップスはnorm.tipsというオブジェクトに格納済みとする。
x.var.rate <- numeric(50)
for(i in 1:50){
x.var.rate[i] <- var(sample(norm.tips,5,replace=T))/var(sample(norm.tips,5,replace=T))
}
hist(x.var.rate)
**問題 2
x.diff <- numeric(10000)
box <- c(0,1)
for(i in 1:10000){
black <- sum(sample(box, 50, replace=T))
white <- 50 - black
x.diff[i] <- (25-black)^2/25 + (25-white)^2/25
}
library(MASS)
truehist(x.diff,xlim=c(0,20),ylim=c(0,1.5))
**問題 3
x.mean.var <- numeric(100) #平均値のバラつきの大きさ(分散)の入れ物
for(j in 1:100){
x.mean <- numeric(100)
for(i in 1:100){
x.mean[i] <- mean(sample(norm.tips, j, replace=T))
}
x.mean.var[j] <- var(x.mean)
}
plot(x.mean.var)
plot(x.mean.var*1:100)
----
