ダメ出し：例題をやってみる その2

今回は，R プログラミングではなく，統計学の方でダメ出し
例題をやってみる　その2　で，
> ここでは、抗生物質を使うとプラセボより感染率は小さくなるはずだ、という期待が込められているので片側検定を行なっていると思われる。
のに
> （ここではカイ二乗検定を行なっている。）
とは，何を考えているのだろうか?そのようなときには，prop.test を使う。
そして，alternarive="less" を指定する。そうすると，
>ここでp=0.02と返ってきているが、論文では0.01と言っている。
などと，寝言を言わなくても，ちゃんと p= 0.01043 になる。

---

 

> prop.test(c(32, 53), c(488, 484), alternative="less")

    2-sample test for equality of proportions with continuity
    correction

data:  c(32, 53) out of c(488, 484)
X-squared = 5.3389, df = 1, p-value = 0.01043
alternative hypothesis: less
95 percent confidence interval:
 -1.00000000 -0.01212703
sample estimates:
    prop 1     prop 2
0.06557377 0.10950413

ダメ出し：例題をやってみる その1の続き ブートストラップ法

相変わらず，「空白などという無駄なものは置くものか」という，読む気のしないグッチャリプログラム

元のプログラムは 98 秒かかる（作者の環境では 3 分ぐらいかかるそうだ）が，

        pvalue <- sum(apply(boot, 1, mean) > 37)/l

を

        pvalue <- sum(rowMeans(boot) > 37)/l

に変えてやるだけで，6 秒ほどで計算が終わる。まあ，約 15 倍速くなると言うわけだ。
二重の for も，どうにかしたいと思うが，ヘタに手を入れると，かえって遅くなるのでそのまま。

以下はそのプログラム。

system.time({
set.seed(888)
tm <- c(37, 36, 37.1, 37.1, 36.2, 37.3, 36.8, 37, 36.3, 36.9, 36.7, 36.8) # 本の値
L <- (1:10)*1000 # リサンプリングする回数
trials <- 100 # 1リサンプリングあたりのシミュレーション回数
data <- matrix(0, trials, length(L))
for(l in L) {
    for(i in 1:trials) {
        boot <- matrix(sample(tm, size=l*length(tm), replace=TRUE), nc=length(tm))
        pvalue <- sum(rowMeans(boot) > 37)/l
        data[i, which(L == l)] <- pvalue
    }
}
matplot(t(data), type="p", pch=1, cex=0.2, col=1, axes=FALSE)
axis(1, 1:length(L), L)
axis(2)
lines(1:length(L), apply(data, 2, mean), col=4, lwd=3)
})
   ユーザ   システム       経過 
    6.249      0.163      6.400

ついでに，11/28　MIKUセミナーも
同じく apply(..., 1, mean) を rowMeans にするだけで倍速になる

> system.time({
+ set.seed(888)
+ L <- 100000 # リサンプリングする回数
+ tm <- c(37, 36, 37.1, 37.1, 36.2, 37.3, 36.8, 37, 36.3, 36.9, 36.7, 36.8) # 本の値
+ data <- matrix(sample(tm, size=L*length(tm), replace=TRUE), nc=length(tm)) # リサンプリングします
+ datamean <- rowMeans(data) # 平均値
+ sum(datamean > 37)/L # が、帰無仮説37度より高いか
+ hist(datamean, nclass=1000) # ヒストグラム
+
+ res <- cumsum(rowMeans(data) > 37)/(1:L)
+ plot(res, cex=0.1) # リサンプリングごとの、mu>37の確率。
+ })
   ユーザ   システム       経過 
     2.401      0.045      2.448