目的

ポアソン分布への適合度の検定を行う

使用法

poissondist(d)
summary.poissondist(obj, digits=5)
plot.poissondist(obj, ...)

引数

d       度数分布ベクトル
obj     poissondist 関数が返すオブジェクト
digits  表示桁数
...     barplot に渡す引数

ソース

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

# ポアソン分布への適合度の検定
poissondist <- function(d)                              # 度数分布ベクトル
{
        data.name <- deparse(substitute(d))
        method <- "ポアソン分布への適合度の検定"
        k <- length(d)                                  # 階級数
        n <- sum(d)                                     # データ数
        k1 <- k-1
        x <- 0:k1
        lambda <- sum(d*x)/n                            # 平均値（=λ）

        p <- exp(-lambda)*lambda^x/factorial(x)         # 確率
        p[k] <- 1-sum(p[-k])                            # 最後の確率は合計が 1 になるように
        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, lambda=lambda), method=method,
                data.name=data.name, table=table), class=c("htest", "poissondist")))
}
# summary メソッド
summary.poissondist <- function(obj,                    # poissondist が返すオブジェクト
                                digits=5)
{
        table <- obj$table
        colnames(table) <- c("階級", "度数", "確率", "期待値")
        cat("\n適合度\n\n")
        print(table, digits=digits, row.names=FALSE)
}
# plot メソッド
plot.poissondist <- function(   obj,                    # poissondist が返すオブジェクト
                                ...)                    # 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)
}       


使用例

> d <- c(27, 61, 77, 71, 54, 35, 20, 11, 6, 2, 1)
> (a <- poissondist(d))

	ポアソン分布への適合度の検定

data:  d 
X-squared = 14.1437, df = 8, p-value = 0.0781
sample estimates:
         n     lambda 
365.000000   2.991781 

> summary(a)

適合度

 階級 度数      確率   期待値
    0   27 0.0501980 18.32226
    1   61 0.1501813 54.81618
    2   77 0.2246548 81.99899
    3   71 0.2240393 81.77434
    4   54 0.1675691 61.16272
    5   35 0.1002660 36.59709
    6   20 0.0499957 18.24841
    7   11 0.0213680  7.79932
    8    6 0.0079910  2.91673
    9    2 0.0026564  0.96958
   10    1 0.0010805  0.39437
 
> plot(a)



 解説ページ