乱数の種
乱数は結果が予測できないことが望ましく毎回異なる結果が出るが、 seed
を元に値を生成するので set.seed()
を用いると同じ結果を再現することもできる。
set.seed(0)
runif(3)
runif(3)
1回目と2回目は違う結果
set.seed(0)
runif(6)
確率分布のシミュレーション
0.2
の確率で起こるベルヌイ試行を10回行った時の各回数の起こる確率を計算する。
```r
b01 <- tibble(x=0:10,y=dgeom(0:10, 0.2) )
ggplot(b01) +
geom_bar( aes(x= x, y=y ), stat=\identity\)
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
# 正規乱数のシミュレーション
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuc2V0LnNlZWQoMClcbmdncGxvdCh0aWJibGUoeCA9IHJub3JtKDEwMDAwKSApICkgK1xuICAgICAgICAgZ2VvbV9oaXN0b2dyYW0oIGFlcyh4PXgsIHk9YWZ0ZXJfc3RhdChkZW5zaXR5KSApICwgYmlucz01MCApICtcbiAgICAgICAgIGdlb21fZnVuY3Rpb24oZnVuID0gZG5vcm0sICAgY29sb3VyID0gXCJyZWRcIilcbmBgYCJ9 -->
```r
set.seed(0)
ggplot(tibble(x = rnorm(10000) ) ) +
geom_histogram( aes(x=x, y=after_stat(density) ) , bins=50 ) +
geom_function(fun = dnorm, colour = "red")

標本分布のシミュレーション
#シミュレーションのためのパラメータの設定
trial <- 20000
u1 <- 2; u2 <- 3
sigma1 <- 4; sigma2 <- 5
n1 <- 60; n2 <- 70
#
a <- vector("numeric", length = trial )
for(i in 1:trial ) {
x <- rnorm(n1, u1, sigma1); y <- rnorm(n2, u2, sigma2)
W1 <- sum( (x - mean(x) )^2 ) /sigma1^2
W2 <- sum( (y - mean(y) )^2 ) /sigma2^2
a[i] <- ( W1 / (n1-1) ) / ( W2 / (n2-1) )
}
df1 <-tibble( x = a )
ggplot(data=df1,aes(x=x)) +
geom_histogram(aes(y=after_stat(density) ), bins=50 ) +
geom_function(fun =df, args = c(df1 = n1-1, df2 = n2 -1), colour= "red" )

LS0tDQp0aXRsZTogIueiuueOh+WIhuW4gyINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNCg0KYGBge3IsZWNobz1GQUxTRSxtZXNzYWdlPUZBTFNFfQ0KbGlicmFyeSh0aWR5dmVyc2UpDQpgYGANCg0KIyDkubHmlbDjga7nqK4gDQoNCuS5seaVsOOBr+e1kOaenOOBjOS6iOa4rOOBp+OBjeOBquOBhOOBk+OBqOOBjOacm+OBvuOBl+OBj+avjuWbnueVsOOBquOCi+e1kOaenOOBjOWHuuOCi+OBjOOAgQ0Kc2VlZCDjgpLlhYPjgavlgKTjgpLnlJ/miJDjgZnjgovjga7jgacgYHNldC5zZWVkKClgIOOCkueUqOOBhOOCi+OBqOWQjOOBmOe1kOaenOOCkuWGjeePvuOBmeOCi+OBk+OBqOOCguOBp+OBjeOCi+OAgg0KDQpgYGB7cixzZXRzZWVkMX0NCnNldC5zZWVkKDApDQpydW5pZigzKQ0KcnVuaWYoMykNCmBgYA0KMeWbnuebruOBqDLlm57nm67jga/pgZXjgYbntZDmnpwNCmBgYHtyLHNldHNlZXQyfQ0Kc2V0LnNlZWQoMCkNCnJ1bmlmKDYpDQpgYGANCg0KIyDnorrnjofliIbluIPjga7jgrfjg5/jg6Xjg6zjg7zjgrfjg6fjg7MNCg0KYDAuMmDjgIDjga7norrnjofjgafotbfjgZPjgovjg5njg6vjg4zjgqToqabooYzjgpIxMOWbnuihjOOBo+OBn+aZguOBruWQhOWbnuaVsOOBrui1t+OBk+OCi+eiuueOh+OCkuioiOeul+OBmeOCi+OAgg0KYGBge3IgYmlub219DQpiMDEgPC0gIHRpYmJsZSh4PTA6MTAseT1kZ2VvbSgwOjEwLCAwLjIpICApDQpnZ3Bsb3QoYjAxKSArDQoJIGdlb21fYmFyKCBhZXMoeD0geCwgeT15ICksICBzdGF0PSJpZGVudGl0eSIpDQoNCmBgYA0KIyDmraPopo/kubHmlbDjga7jgrfjg5/jg6Xjg6zjg7zjgrfjg6fjg7MNCmBgYHtyIHJub3JtfQ0Kc2V0LnNlZWQoMCkNCmdncGxvdCh0aWJibGUoeCA9IHJub3JtKDEwMDAwKSApICkgKw0KICAgICAgICAgZ2VvbV9oaXN0b2dyYW0oIGFlcyh4PXgsIHk9YWZ0ZXJfc3RhdChkZW5zaXR5KSApICwgYmlucz01MCApICsNCiAgICAgICAgIGdlb21fZnVuY3Rpb24oZnVuID0gZG5vcm0sICAgY29sb3VyID0gInJlZCIpDQpgYGANCg0KIyDmqJnmnKzliIbluIPjga7jgrfjg5/jg6Xjg6zjg7zjgrfjg6fjg7MNCmBgYHtyIEZkaXN0fQ0KI+OCt+ODn+ODpeODrOODvOOCt+ODp+ODs+OBruOBn+OCgeOBruODkeODqeODoeODvOOCv+OBruioreWumg0KdHJpYWwgPC0gMjAwMDANCnUxIDwtIDI7ICAgICB1MiA8LSAzDQpzaWdtYTEgPC0gNDsgc2lnbWEyIDwtIDUNCm4xIDwtIDYwOyAgICBuMiA8LSA3MA0KIw0KYSA8LSB2ZWN0b3IoIm51bWVyaWMiLOOAgGxlbmd0aOOAgD3jgIB0cmlhbOOAgCkNCmZvcihpICBpbiAgMTp0cmlhbCApIHsNCiAgeCA8LSBybm9ybShuMSzjgIB1MSzjgIBzaWdtYTEpOyB5ICA8LSBybm9ybShuMizjgIB1MizjgIBzaWdtYTIpDQogIFcxIDwtIHN1bSggKHggLSBtZWFuKHgpICleMiAgKSAvc2lnbWExXjINCiAgVzIgPC0gc3VtKCAoeSAtIG1lYW4oeSkgKV4yICApIC9zaWdtYTJeMg0KICBhW2ldIDwtICAoIFcxIC8gKG4xLTEpICApICAvICggVzIgIC8gKG4yLTEpICkgDQp9ICANCmRmMSA8LXRpYmJsZSggeCA9IGEgKQ0KZ2dwbG90KGRhdGE9ZGYxLGFlcyh4PXgpKSArIA0KICBnZW9tX2hpc3RvZ3JhbShhZXMoeT1hZnRlcl9zdGF0KGRlbnNpdHkpICksIGJpbnM9NTAgICkgKw0KICBnZW9tX2Z1bmN0aW9uKGZ1biA9ZGYsICBhcmdzID0gYyhkZjEgPSBuMS0xLCBkZjIgPSBuMiAtMSksIGNvbG91cj0gInJlZCIgKSANCmBgYA0KDQoNCg0KDQo=