「ノーマル・チップスを用いたシミュレーションを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軸にする。
----
*問題の解答
※解答は例です。もっと効率的でスマートな方法があるやもしれません。
**問題 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)
少し長くなりますが、以下の様にすると読みやすくなり、条件を変えたときの観察もしやすくなります。
try <- 50 #試行回数
n1 <- 5 #分母の分散を計算するときのサンプルサイズ
n2 <- 5 #分子の分散を計算するときのサンプルサイズ
x.var.rate <- numeric(try) #try個分の長さを持った空の数値ベクトル。データ記録用。
for(i in 1:try){
var.1 <- var(sample(norm.tips, n1, replace=T)) #分母にする分散
var.2 <- var(sample(norm.tips, n2, replace=T)) #分子にする分散
F.val <- var.1/var.2 #分散比F
x.var.rate[i] <- F.val #分散比Fの記録
}
library(MASS) #truehist関数を使うためにMASSパッケージを読み込む
truehist(x.var.rate) #分散比Fの分布を(ちょっといい感じの)ヒストグラムとして描画
例えばここでtryというオブジェクトに試行回数を入れています。試行回数はデータ記録用のベクトルを作成するときと、for文を作るときの2回使用します。オブジェクトに試行回数を入れるようにしておけば、試行回数を変えたくなったときはそこだけの変更で事足ります。2回以上使用する数値はミスを防ぐためにオブジェクトにしておくべきでしょう。
また、1回しか使用しないn1、n2についてもオブジェクトにし、さらにコメントでその意味の解説を付け加えています。こうしておくと後でプログラムを見たときにもn1やn2がどういった意味を持った数値なのかが一目瞭然です。また、tryなどと一緒にプログラムの最初のほうにまとめておくと、後で数値をイロイロ変えて実験したくなったときに便利です。
というわけでプログラム中で使用する数値(特に変更の可能性、または余地があるもの)は使用回数にかかわらずオブジェクトにしておきましょう。
**問題 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))
こいつも可読性重視で長めに書き換えると次のようになります。
try <- 1000 #試行回数
color <- 2 #箱の中のボールの種類
lot <- 50 #一回の試行で引くボールの個数
box <- 1:color #箱の作成（ボールの割合は各色等しい)
x.diff <- numeric(try) #データを記録する入れ物の作成
for(i in 1:try){
result <- sample(box, lot, replace=T) #ボールをlot個抽出
chi <- 0 #ずれの測度(カイ二乗値)をリセット
for(j in 1:color){ #期待値とのずれの測度を累積してカイ二乗値を計算
chi <- chi + (length(result[result==j])-(lot/color))^2/(lot/color)
}
x.diff[i] <- chi #カイ二乗値を記録
}
library(MASS) #いーかんじのヒストグラムを描くためにMASS読み込み
truehist(x.diff, xlim=c(0,10), ylim=c(0,0.3)) #いーかんじのヒストグラムを描く
par(new=T) #一回だけ上書きを許可
curve(dchisq(x,color),0,10, ylim=c(0,0.3)) #自由度がボールの種類分のカイ二乗分布を描く
ついでに箱の中のボールの色を増やせるようにしました。ここでは、ボールの色を1、2、3、…といった数値に対応させています。ここから抽出するとこれらの値が混ざったベクトルが返ってきますから、
> sample(box, 10, replace=T)
[1] 1 1 1 1 2 1 2 1 2 1
このベクトルの中にある1の数、2の数、3の数…をそれぞれカウントしてやればいいわけです。これは[[添え字による指定>http://www39.atwiki.jp/stat_semi/pages/18.html#id_8d8776d4]]とlength関数を組み合わせることで行うことが可能です。length関数は引数として与えられたベクトルの要素の個数をカウントします。ですから、例えば次のように入力すると、サンプリングの結果の中で1が出た回数だけをカウントできます。
> result <- sample(box, 10, replace=T)
> length(result[result==1])
[1] 7
「resultというオブジェクトの中で値が1に等しいものだけを数えろ！」という命令をしているわけですね。
で、あとは教科書の定義にしたがって「食い違いの測度」すなわちカイ二乗値を算出しています。和を計算するのにfor文を使うことで項の増加に対応しています。数式で言うと$$\Sigma$$を使うのと同じで、たくさん書くのが大変なのでちょっとズルをしているわけです。
んでそのあとはいつもどおりヒストグラムを書くわけですが、ついでに自由度＝色数のカイ二乗分布の密度関数を上から重ね描きするようにしてみました。最初のcolorを色々変えてみると、実際にカイ二乗値の分布がカイ二乗分布に従うことが観察できると思います。
ここでは「2種類のボールの抽出」を「n種類のボールの抽出」へ一般化してプログラムを書き直しました。このように「特定の状況の解決」のためのプログラムを「一般的な状況の解決」のために書き直しておくと後で再利用がしやすくなります。要するになるべくたくさんのパラメータを簡単に再設定できるようにしておけということです。
**問題 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)
単に本文中の関数を100回繰り返しただけです。
ポイントは
-記録するのは平均値ではなく、「平均値のバラつき」
-「平均値のバラつき」は「サンプル数」に対応してひとつだけ記録する
-「平均値のバラつき」と「サンプル数」の対応を観察したいので、描くのは散布図
-散布図には「平均値のバラつき」ベクトルだけ与えれば、x軸は自動的に1, 2, 3…というベクトルが補完される
といったところでしょうか。この式はfor文の中でfor文をまわしているので繰り返し回数を少し増やしただけでもかなり処理時間を食います。まあ「少し」といっても0をいくつ増やすか、といったオーダの話なので無茶をしないかぎりは大丈夫です。
----
#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)
少し長くなりますが、以下の様にすると読みやすくなり、条件を変えたときの観察もしやすくなります。
try <- 50 #試行回数
n1 <- 5 #分母の分散を計算するときのサンプルサイズ
n2 <- 5 #分子の分散を計算するときのサンプルサイズ
x.var.rate <- numeric(try) #try個分の長さを持った空の数値ベクトル。データ記録用。
for(i in 1:try){
var.1 <- var(sample(norm.tips, n1, replace=T)) #分母にする分散
var.2 <- var(sample(norm.tips, n2, replace=T)) #分子にする分散
F.val <- var.1/var.2 #分散比F
x.var.rate[i] <- F.val #分散比Fの記録
}
library(MASS) #truehist関数を使うためにMASSパッケージを読み込む
truehist(x.var.rate) #分散比Fの分布を(ちょっといい感じの)ヒストグラムとして描画
例えばここでtryというオブジェクトに試行回数を入れています。試行回数はデータ記録用のベクトルを作成するときと、for文を作るときの2回使用します。オブジェクトに試行回数を入れるようにしておけば、試行回数を変えたくなったときはそこだけの変更で事足ります。2回以上使用する数値はミスを防ぐためにオブジェクトにしておくべきでしょう。
また、1回しか使用しないn1、n2についてもオブジェクトにし、さらにコメントでその意味の解説を付け加えています。こうしておくと後でプログラムを見たときにもn1やn2がどういった意味を持った数値なのかが一目瞭然です。また、tryなどと一緒にプログラムの最初のほうにまとめておくと、後で数値をイロイロ変えて実験したくなったときに便利です。
というわけでプログラム中で使用する数値(特に変更の可能性、または余地があるもの)は使用回数にかかわらずオブジェクトにしておきましょう。
**問題 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))
こいつも可読性重視で長めに書き換えると次のようになります。
try <- 1000 #試行回数
color <- 2 #箱の中のボールの種類
lot <- 50 #一回の試行で引くボールの個数
box <- 1:color #箱の作成（ボールの割合は各色等しい)
x.diff <- numeric(try) #データを記録する入れ物の作成
for(i in 1:try){
result <- sample(box, lot, replace=T) #ボールをlot個抽出
chi <- 0 #ずれの測度(カイ二乗値)をリセット
for(j in 1:color){ #期待値とのずれの測度を累積してカイ二乗値を計算
chi <- chi + (length(result[result==j])-(lot/color))^2/(lot/color)
}
x.diff[i] <- chi #カイ二乗値を記録
}
library(MASS) #いーかんじのヒストグラムを描くためにMASS読み込み
truehist(x.diff, xlim=c(0,10), ylim=c(0,0.3)) #いーかんじのヒストグラムを描く
par(new=T) #一回だけ上書きを許可
curve(dchisq(x,color-1),0,10, ylim=c(0,0.3)) #自由度がボールの種類分-1のカイ二乗分布を描く
ついでに箱の中のボールの色を増やせるようにしました。ここでは、ボールの色を1、2、3、…といった数値に対応させています。ここから抽出するとこれらの値が混ざったベクトルが返ってきますから、
> sample(box, 10, replace=T)
[1] 1 1 1 1 2 1 2 1 2 1
このベクトルの中にある1の数、2の数、3の数…をそれぞれカウントしてやればいいわけです。これは[[添え字による指定>http://www39.atwiki.jp/stat_semi/pages/18.html#id_8d8776d4]]とlength関数を組み合わせることで行うことが可能です。length関数は引数として与えられたベクトルの要素の個数をカウントします。ですから、例えば次のように入力すると、サンプリングの結果の中で1が出た回数だけをカウントできます。
> result <- sample(box, 10, replace=T)
> length(result[result==1])
[1] 7
「resultというオブジェクトの中で値が1に等しいものだけを数えろ！」という命令をしているわけですね。
で、あとは教科書の定義にしたがって「食い違いの測度」すなわちカイ二乗値を算出しています。和を計算するのにfor文を使うことで項の増加に対応しています。数式で言うと$$\Sigma$$を使うのと同じで、たくさん書くのが大変なのでちょっとズルをしているわけです。
んでそのあとはいつもどおりヒストグラムを書くわけですが、ついでに自由度＝色数-1のカイ二乗分布の密度関数を上から重ね描きするようにしてみました。最初のcolorを色々変えてみると、実際にカイ二乗値の分布がカイ二乗分布に従うことが観察できると思います。
ここでは「2種類のボールの抽出」を「n種類のボールの抽出」へ一般化してプログラムを書き直しました。このように「特定の状況の解決」のためのプログラムを「一般的な状況の解決」のために書き直しておくと後で再利用がしやすくなります。要するになるべくたくさんのパラメータを簡単に再設定できるようにしておけということです。
**問題 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)
単に本文中の関数を100回繰り返しただけです。
ポイントは
-記録するのは平均値ではなく、「平均値のバラつき」
-「平均値のバラつき」は「サンプル数」に対応してひとつだけ記録する
-「平均値のバラつき」と「サンプル数」の対応を観察したいので、描くのは散布図
-散布図には「平均値のバラつき」ベクトルだけ与えれば、x軸は自動的に1, 2, 3…というベクトルが補完される
といったところでしょうか。この式はfor文の中でfor文をまわしているので繰り返し回数を少し増やしただけでもかなり処理時間を食います。まあ「少し」といっても0をいくつ増やすか、といったオーダの話なので無茶をしないかぎりは大丈夫です。
----
