目的

ポリヤ・エッゲンベルガー分布の確率を計算する

使用法

(1) PolyaEggenberger(x, n, p, delta)
(2) PolyaEggenberger2(x, lambda, r)

引数

x       確率変数

n       標本サイズ
p       母比率
delta   追加する割合

lambda  n×p
r       n×delta

ソース

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

# ポリヤ・エッゲンベルガー分布の確率を計算する
# n, p, δ を与える場合
PolyaEggenberger <- function(        x,                      # 確率変数
                                n,                      # 標本サイズ
                                p,                      # 母比率
                                delta)                  # 追加する割合 δ
{
        exp(                                            # 対数で計算して最後に逆対数をとる
                sum(lchoose(n, x),
                sapply(0:x, function(i) ifelse(i == 0, 0, log(p+(i-1)*delta)-log(1+(i-1)*delta))),
                sapply((x+1):n, function(i) log(1-p+(i-x-1)*delta)-log(1+(i-1)*delta)))
        )
}
# λ, r を与える場合
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)
        )
}


使用例

> n <- 800	# n, p, δ を与える場合
> p <- 0.00373625
> delta <- 0.000322
> PolyaEggenberger(0, n, p, delta)
[1] 0.06964384
> PolyaEggenberger(1, n, p, delta)
[1] 0.1660618
> PolyaEggenberger(2, n, p, delta)
[1] 0.2148316
> PolyaEggenberger(3, n, p, delta)
[1] 0.1997851
> sum(sapply(0:15, PolyaEggenberger, n, p, delta))
[1] 0.9999915

> lambda <- 2.989	# λ, r を与える場合
> r <- 0.2576
> PolyaEggenberger2(0, lambda, r)
[1] 0.06998131
> PolyaEggenberger2(1, lambda, r)
[1] 0.1663280
> PolyaEggenberger2(2, lambda, r)
[1] 0.2146949
> sum(sapply(0:15, PolyaEggenberger2, lambda, r))
[1] 0.9999908

 解説ページ