ダメ出し：添え字はそのまま使う

「楽しい　統計学セミナー　- 証拠を科学する方法 -」の「第2講」 において

試行回数 N=10～100 に対する二項分布の確率関数を

x <- numeric(91)
for (i in 10:100) {
 x[i - 9] <- choose(i, 3) * (1/10)^(3) * (1- (1/10))^(i - 3)
}
plot(x)
abline(v = 20.5, col="blue")

のように，書いている。10 ~ 100 を添え字として 1 ～ 91 を使おうということで numeric(91) とか for (i in 10:100) とか x[i - 9] とか，様々な注意が必要。しかし，このようなことはちょっと油断するとバグの混入原因。縦軸（確率）を納める変数名が x なのも混乱を増幅するが。

実際，掲載されたプログラムが
for (i in 10:91) {
 x[i - 9] <- > for (i in 10:100) {
となっていたり，確率が最大になる index の値を図に示すと，index=20, 21 のときに最大になるというように見えるが，求める数値は 29 と 30 なのだ，のように，ぼろが出る。

添え字は，0 やマイナスの値にならない限りそのまま使う方がよい。

y <- numeric(100)
n <- 10:100
for (i in n) {
 y[i] <- choose(i, 3) * (1/10)^(3) * (1- (1/10))^(i - 3)
}
plot(n, y[n])
abline(v = 29.5, col="blue")

なお，二項分布の確率関数は dbinom で計算できる

> x <- 10:100
> y <- dbinom(3, x, 1/10)
> plot(x, y, type="l", ylim=c(0, 0.25))
> (index <- which.max(y))
[1] 20
> abline(v = x[index] + 0.5, col="blue")
> (max.value <- y[index])
[1] 0.2360879322323426