目的

ポリヤ・エッゲンベルガー分布への適合度の検定を行う

使用法

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

引数

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

ソース

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

# ポリヤ・エッゲンベルガー分布への適合度の検定
PolyaEggenbergerdist <- function(d)                     # 度数分布ベクトル
{
        PolyaEggenberger2 <- function(  x,              # 確率変数
                                        lambda,         # λ = n*p
                                        r)              # r = n*δ
        {
                exp(                                    # 対数で計算して最後に逆対数をとる
                        sum(sapply(0:(x-1), function(i) ifelse(i < 0, 0, log(1+i*r/lambda))))
                        +x*log(lambda)
                        -lgamma(x+1)
                        +(-lambda/r-x)*log(1+r)
                )
        }
#===
        data.name <- deparse(substitute(d))
        method <- "ポリヤ・エッゲンベルガー分布への適合度の検定"
        k <- length(d)                                  # 階級数
        n <- sum(d)                                     # データ数
        k1 <- k-1
        x <- 0:k1
        lambda <- sum(d*x)/n                            # 平均値（=λ）
        r <- sum(d*(x-lambda)^2)/n/lambda-1             # パラメータ r

        p <- sapply(x, PolyaEggenberger2, lambda, r)    # 確率
        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-3                                       # 自由度
        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, r=r), method=method,
                data.name=data.name, table=table), class=c("htest", "PolyaEggenbergerdist")))
}
# summary メソッド
summary.PolyaEggenbergerdist <- function(obj,           # PolyaEggenbergerdist が返すオブジェクト
                                digits=5)
{
        table <- obj$table
        colnames(table) <- c("階級", "度数", "確率", "期待値")
        cat("\n適合度\n\n")
        print(table, digits=digits, row.names=FALSE)
}
# plot メソッド
plot.PolyaEggenbergerdist <- function(  obj,            # PolyaEggenbergerdist が返すオブジェクト
                                ...)                    # 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 <- PolyaEggenbergerdist(d))

	ポリヤ・エッゲンベルガー分布への適合度の検定

data:  d 
X-squared = 0.6248, df = 8, p-value = 0.9997
sample estimates:
          n      lambda           r 
365.0000000   2.9917808   0.2682924 

> summary(a)

適合度

 階級 度数      確率  期待値
    0   27 0.0706286 25.7794
    1   61 0.1666061 60.8112
    2   77 0.2141257 78.1559
    3   71 0.1985646 72.4761
    4   54 0.1486017 54.2396
    5   35 0.0952554 34.7682
    6   20 0.0542416 19.7982
    7   11 0.0281137 10.2615
    8    6 0.0134934  4.9251
    9    2 0.0060738  2.2170
   10    1 0.0042954  1.5678

> plot(a)



 解説ページ