適合度の検定（exact test）　　　　　Last modified: Jan 14, 2006

目的
　カイ二乗分布を用いる適合度の検定の exact test である。
　データの大きさが固定された，可能性のある全ての度数分布表を発生させ，その度数分布表についてカイ二乗分布を用いる適合度検定の検定統計量を求め，観察された度数分布表における検定統計量より大きい場合に，その度数分布表の生起確率（多項分布に基づく）を累積したものを，正確な P 値とするものである。
　なお，度数分布表によっては計算量が多くなり実用的な時間内に計算が終わらない可能性があるので，注意されたい。そのような場合にはネット上で計算すればもう少し計算時間は短い。 
使用法

gft(x)
gft(x，p)

引数

x    度数
p    理論比。合計が 1 でなくても良い p=c(9, 3, 3, 1) のように指定できる
     省略されたときには，一様性の検定 p=rep(1/length(x), length(x)) を行う

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/gft.R", encoding="euc-jp")

# 適合度の検定（exact test）
gft <- function(o,                                                           # 度数
                p=rep(1/length(o), length(o)))                                  # 理論比
{

        printf <- function(fmt, ...)                                         # C 言語の printf をシミュレートする
        {
                cat(sprintf(fmt, ...))
        }

        gen_tab <- function(y)                                                       # 度数表の発生
        {
                if (y == 1) {
                        x2 <- sum((o-e)^2/e)
                        if (x2 >= x0 || abs(x2-x0) <= EPSILON) {
                                w2 <<- w2+exp(temp-sum(fact[o+1]))*prod(p^o)
                        }
                }
                else {
                        gen_tab(y-1)
                        while (o[1]) {
                                o[y] <<- o[y]+1
                                o[1] <<- o[1]-1
                                gen_tab(y-1)
                        }
                        o[1] <<- o[1]+o[y]
                        o[y] <<- 0
                }
        }

        EPSILON <- 1e-10
        total <- sum(o)                                                              # サンプルサイズ
        p <- p/sum(p)                                                                # 理論比
        e <- p*total                                                         # 期待値
        x0 <- sum((o-e)^2/e)                                                 # カイ二乗値
        n <- length(o)                                                               # カテゴリー数
        p_val <- pchisq(x0, n-1, lower.tail=FALSE)                           # P 値
        printf("カイ二乗値は %g，自由度は %i，P値は %g\n", x0, n-1, p_val)      # 結果出力
        fact <- lfactorial(0:(total+1))                                              # 対数を取った階乗表
        temp <- fact[total+1]
        w2 <- 0
        o <- numeric(n)
        o[1] <- total
        gen_tab(n)                                                              # 正確検定
        printf("正確なP値は %g\n", w2)
}


使用例

gft(c(10, 12, 9, 4, 13, 8))
gft(c(10, 12, 9, 4, 13, 8), p=rep(1, 6))
gft(c(29, 12, 8, 2), p=c(9, 3, 3, 1))

出力結果例

> gft(c(10, 12, 9, 4, 13, 8))
カイ二乗値は 5.5，自由度は 5，Ｐ値は 0.357946
正確なＰ値は 0.370005

> gft(c(10, 12, 9, 4, 13, 8), p=rep(1, 6))
カイ二乗値は 5.5，自由度は 5，Ｐ値は 0.357946
正確なＰ値は 0.370005

> gft(c(29, 12, 8, 2), p=c(9, 3, 3, 1))
カイ二乗値は 1.32244，自由度は 3，Ｐ値は 0.723811
正確なＰ値は 0.741471
・ 解説ページ

・ 直前のページへ戻る　　・ E-mail to Shigenobu AOKI

Made with Macintosh