裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

大きい数の積や商(その2)

2018年06月08日 | ブログラミング

超幾何分布


要因 a, b が独立な場合,以下のような分割表が得られる確率を求める

          B     not(B)   sum
A         x       o       m
not(A)  (k-x)     o       n
sum       k       o     (m+n)

超幾何分布により,求める確率は xCm * xCk-x / m+nCk

R の dhyper の引数はちょっと変で(?),x が生じる確率は dhyper(x, m, n, k)

          B     not(B)   sum
A         2       o       10
not(A)   (1)      o       5
sum       3       o      (15)

> x = 2; m = 10; n = 5; k = 3
> dhyper(x, m, n, k)
[1] 0.494505494505494
> choose(m, x) * choose(n, k-x) / choose(m+n, k)
[1] 0.494505494505495
> exp(lchoose(m, x) + lchoose(n, k-x) - lchoose(m+n, k))
[1] 0.494505494505495

しかし,

> x = 200; m = 2000; n = 5000; k = 600
> dhyper(x, m, n, k)
[1] 0.00104601828356522
> choose(m, x) * choose(n, k-x) / choose(m+n, k)
[1] NaN
> exp(lchoose(m, x) + lchoose(n, k-x) - lchoose(m+n, k))
[1] 0.00104601828356506

となり,普通のやり方では精度が足りない。

Rmpfr パッケージを用いて,単純に計算すると以下のようになる。

library(Rmpfr)
precBits=13
dhyper2 = function(x, m, n, k) {
  choose2 = function(a, b) {
    prod(mpfr(1:a, precBits)) / prod(mpfr(1:b, precBits)) / prod(mpfr(1:(a-b), precBits))
  }
  as.numeric(choose2(m, x) * choose2(n, k-x) / choose2(m+n, k))
}

> x = 200; m = 2000; n = 5000; k = 600
> dhyper(x, m, n, k)
[1] 0.00104601828356522
> dhyper2(x, m, n, k)
[1] 0.00104601828356522

dhyper2 の実行時間は 0.15 秒



コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村