/var/log/messages

debugging with sixth sense

パラメトリックブートストラップ検定

何故か 5-4 から再開。つうか

良書ですね。Python で云々とか緑な本とか読んだりしたいです。

とりあえず、何となく 5-4-8 以降をもごもごしつつ控えを。

データを用意して統計モデルを作る

以下?

1
2
3
4
5
6
7
d2 <- data.frame(
   length = c(2, 4, 10, 12, 6, 8),
   medicine = c("A", "A", "B", "B", "C", "C")
)

modelANOVA <- lm(length ~ medicine, data=d2)
modelNull <- lm(length ~ 1, data=d2)

標本から F 比を計算

凄い、

1
anova(modelANOVA)

から F value だけ取り出せるのか。

1
2
> anova(modelANOVA)$"F value"[1]
[1] 16

R は配列参照が 1 origin なのですね。

シミュレーションデータ作成

こう、なのかな。

1
2
3
> simData <- cbind(
+ simulate(modelNull, 1),
+ d2$medicine)

で、以下なのかな。

1
2
3
4
5
6
7
8
9
> colnames(simData) <- c("length", "medicine")
> simData
     length medicine
1  7.073384        A
2  9.022781        A
3  1.486284        B
4  2.100310        B
5 10.555131        C
6  5.317685        C

ええと、値が変わらないのですが良いのかどうか。simulate すると値は変わりますね。

F 比を多数計算する

合わせると以下なのかな。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Nsim <- 10000
sim <- simulate(modelNull, Nsim)

simFValue <- numeric()

for(i in 1:Nsim){
  simData <- cbind(
    sim[i],
    c("A", "A", "B", "B", "C", "C")
  )
  colnames(simData) <- c("length", "medicine")
  model <- lm(length ~ medicine, data=simData)
  simFValue[i] <- anova(model)$"F value"[1]
}

で、値を見ても意味分からないので、以下で図示。

1
2
3
4
hist(
  subset(simFValue, simFValue < 20),
  main="F 比のヒストグラム"
)

スクリプトがバグッててアレ。繰返しにブラケット付けるの忘れてました。

p 値の算出

16 以上の F 比の数。

1
2
> sum(simFValue>=16)
[1] 262

Nsim で割れば p 値が出る模様。

1
2
> sum(simFValue>=16)/Nsim
0.0262

正規分布から派生した確率分布

正規分布偉い、らしい。とは言え、「母集団分布が推定できるわけではない」らしい。

つうかこの 5-5 章はちゃんと確認した方が良さげ。何度か目を通して別途必要があればエントリ入れます。

Comments