目的

二項分布への適合度の検定を行う

使用法

binomdist(d, x, size)
summary.binomdist(obj, digits=5)
plot.binomdist(obj, ...)

引数

d       度数分布ベクトル
x       階級値ベクトル
size    試行回数
obj     binomdist 関数が返すオブジェクト
digits  表示桁数
...     barplot に渡す引数

ソース

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

# 二項分布への適合度の検定
binomdist <- function(       d,                                      # 度数分布ベクトル
                        x,                                      # 階級値ベクトル
                        size)                                   # 試行回数
{
        data.name <- paste( "\n度数分布ベクトル：", deparse(substitute(d)),
                        "\n階級値ベクトル：　", deparse(substitute(x)),
                        "\n試行回数：　　　　", deparse(substitute(size)), sep="")
        method <- "二項分布への適合度の検定"
        k <- length(d)                                               # 階級数
        if (k != length(x)) {
                stop("度数ベクトル d と階級値ベクトル x の長さが違います")
        }
        else if (any(floor(d) != d)) {
                stop("度数ベクトル中に小数値があります")
        }
        else if (any(d < 0)) {
                stop("度数ベクトル中に負の値があります")
        }
        else if (any(x > size)) {
                stop("階級値ベクトル中に試行数より大きい数値があります")
        }
        else if (any(x < 0)) {
                stop("階級値ベクトル中に負の数値があります")
        }
        else if (any(floor(x) != x)) {
                stop("階級値ベクトル中に小数値があります")
        }
        n <- sum(d)                                          # サンプルサイズ
        prob <- sum(d*x)/n/size                                      # 母比率
        p <- dbinom(x, size, prob)                           # 二項分布の確率
        e <- n*p                                             # 期待値
        table <- data.frame(x, d, p, e)
        rownames(table) <- paste("c-", x, sep="")
        while (e[1] < 1) {                                   # 1 未満のカテゴリーの併合
                d[2] <- d[2]+d[1]
                e[2] <- e[2]+e[1]
                d <- d[-1]
                e <- e[-1]
                k <- k-1
        }
        while (e[k] < 1) {                                   # 1 未満のカテゴリーの併合
                d[k-1] <- d[k-1]+d[k]
                e[k-1] <- e[k-1]+e[k]
                d <- d[-k]
                e <- e[-k]
                k <- k-1
        }
        chisq <- sum((d-e)^2/e)                                      # カイ二乗統計量
        df <- k-2                                            # 自由度
        p <- pchisq(chisq, df, lower.tail=FALSE)             # P 値
        names(chisq) <- "X-squared"
        names(df) <- "df"
        return(structure(list(statistic=chisq, parameter=df,
                p.value=p, estimate=c(n=n, probability=prob), method=method,
                data.name=data.name, table=table),
                class=c("htest", "binomdist")))
}
# summary メソッド（適合度に関する結果を表示する）
summary.binomdist <- function(       obj,                            # binomdist が返すオブジェクト
                                digits=5)                       # 表示桁数
{
        table <- obj$table
        colnames(table) <- c("級", "度数", "確率", "期待値")
        cat("\n適合度\n\n")
        print(table, digits=digits, row.names=FALSE)
}
# plot メソッド（観察度数と理論度数を図示する）
plot.binomdist <- function(  obj,                            # binomdist が返すオブジェクト
                                ...)                            # barplot 関数へ渡す引数
{
        table <- obj$table
        nr <- nrow(table)
        barplot(table$d, space=0, ...)                          # 観察度数を barplot で描く
        old <- par(xpd=TRUE)
        points(1:nr-0.5, table$e, pch=3)                        # は理論度数を，記号 + で示す
        text(1:nr-0.5, -strheight("H"), table$x)                # 階級表示
        par(old)
}


使用例

(a <- binomdist(c(2, 14, 20, 34, 22, 8), 0:5, 5))
summary(a)
plot(a)

出力結果例

> (a <- binomdist(c(2, 14, 20, 34, 22, 8), 0:5, 5))
	二項分布への適合度の検定

data:  
度数分布ベクトル：c(2, 14, 20, 34, 22, 8)
階級値ベクトル：　0:5
試行回数：　　　　5 
X-squared = 4.0076, df = 4, p-value = 0.405
sample estimates:
          n probability 
    100.000       0.568 

> summary(a)

適合度

 階級 度数     確率   期待値
    0    2 0.015046   1.5046
    1   14 0.098913   9.8913
    2   20 0.260105  26.0105
    3   34 0.341989  34.1989
    4   22 0.224826  22.4826
    5    8 0.059121   5.9121

> plot(a)



 解説ページ