カイ二乗分布をモンテカルロシミュレーションプログラムで遊ぶ |
chi <- function(k, n){#kは標本数(自由度)、nは試行回数です
#chi()を読み込むモンテカルロシミュレーションプログラムです。
chi_sim <- function(k, n, r){#kは標本数、nは試行回数、rはシミュレーション回数です。
#実行例です。自由度1のカイ二乗分布を、試行回数1000、シミュレーション1000回で出します。
b <- chi_sim(1, 1000, 1000)
b <- data.frame(b)
colnames(b) <- c("value")
p <- ggplot(b, aes(x = value))
p <- p + geom_histogram(breaks = seq(0, 25, by = 0.1))
#上位5%の数値で線を引きます。order()を使ってデータを並べます。
sortlist <- order(b)
bsort <- b[sortlist, ] #sortlistの順にデータを並べ直します
bsort[1000*1000*0.95] #95%の値、すなわちp = 0.05となるカイ値を出します。
p <- p + geom_vline(xintercept = bsort[1000*1000*0.95], colour = "red") #geom_vline()で95%値を垂線で示します。
plot(p)
#自由度4のカイ二乗分布を、試行回数1000、シミュレーション1000回で出します。
bb <- chi_sim(4, 1000, 1000)
bb <- data.frame(bb)
colnames(bb) <- c("value")
pp <- ggplot(bb, aes(x = value))
pp <- pp + geom_histogram(breaks = seq(0, 25, by = 0.1))
#上位5%の数値で線を引くプログラム
sortlist <- order(bb)
bsort <- bb[sortlist, ]
bsort[1000*1000*0.95] #p = 0.05となるカイ値を出します。
pp <- pp + geom_vline(xintercept = bsort[1000*1000*0.95], colour = "red")
plot(pp)
自由度4のp=0.05を与えるカイ値は9.487、よくあった結果が得られています。
近いうちに自由度をいろいろ変えてグラフを動画で出力するプログラムをつくることを試してみようと思っています。何の役に立つかわかりませんが、原理や応用に関する理解は深まるように思います。