乱数の生成 sample table hist replicate

乱数を発生させて、現実をシミュレーションしてみることは大変有益です。

> sample(c("A","B","C"),10,rep=T)
 [1] "B" "A" "A" "A" "A" "A" "C" "B" "B" "B"

?sampleで調べると、次のように使うと出てきます。

sample(x, size, replace = FALSE, prob = NULL)

xの成分が、でたらめにsize個並んだベクトルが返ります。
replaceはxの成分を繰り返し使うかどうかの指定で、省略するとFALSEすなわち、繰り返しなしとなりますから、sizeはlength(x)以下である必要があります。上の例では、rep=Tとしてあります。繰り返しを許すので、個数はいくつでもよくなります。replaceはrepと省略しました。誤解のない範囲で省略できますが、この辺のことはまた後日。　

> x=sample(1:6,100,rep=T)
> x
  [1] 6 2 2 4 6 1 3 4 2 1 4 1 2 2 5 4 2 3 2 3 6 6 5 4 3 3 1
 [28] 3 5 5 4 2 2 6 2 6 4 2 5 4 1 6 2 4 2 4 3 1 2 3 3 5 4 3
 [55] 6 6 1 3 1 3 5 1 3 6 5 5 5 4 1 5 3 3 5 5 3 5 6 1 3 4 6
 [82] 3 6 1 4 1 1 6 6 5 4 3 4 2 1 1 1 4 6 4

さいころを１００回振ることをシミュレーションしました。

> table(x)
x
 1  2  3  4  5  6 
17 15 19 18 15 16 

それぞれの目の出た回数です。

> hist(x)

ヒストグラムです。階級はRが自動調整しますが、今回はその自動調整が裏目に出ています。感覚的にとらえるにはグラフを描かせることが有効です。

> mean(x)
[1] 3.47

平均の理論値は３．５です。

さいころ２個の目の和の分布です。

> y=sample(1:6,100,rep=T)+sample(1:6,100,rep=T)
> table(y)
y
 2  3  4  5  6  7  8  9 10 11 12 
 3  9 10  7 11 23  9 13  7  4  4 
> hist(y)
> mean(y)
[1] 6.86

さいころを１００回振って平均を出すということを、何回も繰り返してみます。

> mean(sample(1:6,100,rep=T))
[1] 3.66
> replicate(10,mean(sample(1:6,100,rep=T)))
 [1] 3.46 3.19 3.17 3.58 3.71 3.44 3.39 3.41 3.48 3.61
> d=replicate(1000,mean(sample(1:6,100,rep=T)))
> hist(d)

最初のは、１００回サイコロを振って平均を出すということを１行で書いたものです。
次は、それを１０回繰り返しました。replicate(回数,命令)で命令を回数繰り返します。
その次は、繰り返しを１０００回にして、その結果を変数dに代入しました。 最後にヒストグラムを描きました。

おおかた平均は３．５近辺になりますが、３以下から４まで、分布していることがわかります。ちょうど３．５にはならないのです。これが現実というものです。１００回振ってもこれだけ幅があるのです。

問題１　コインを投げることをシミュレーションします。コインを１００回投げ、表は"H"、裏は"T"として、結果を文字列ベクトル s としなさい。表裏はHead Tail といいます。

問題２　表が裏の２倍出やすいコインをシミュレーションするにはどうしますか？

問題３　普通のコインを５枚投げて、表の出る回数を数えてください。
　　（文字列ベクトルに対して=="H"でHかどうかがわかります。sumで総和をとります。）

問題４　問題３を１０００回繰り返し、５枚投げたときの表の回数の分布を調べてください。