大きい数の積や商

二項分布 B(x, n, p) = nCx p^x (1-p)^(n-x) の計算
R では，dbinom(x, n, p) で速く，簡単に求まる
これを，自分で計算してみるとどうなるかをやってみる。

式をそのまま記述すると以下のようになる。

options(digits = 15)

f = function(x, n, p) {
  return(choose(n, x) * p ^ x * (1 - p) ^ (n - x))
}

n が小さいうちは，これでも十分である。

> x = 0
> n = 10
> p = 0.5
> dbinom(x, n, p)
[1] 0.0009765625
> f(x, n, p)
[1] 0.0009765625

しかし，n が大きくなると，これじゃダメだ。

> x = 4000
> n = 10000
> p = 0.4
> dbinom(x, n, p)
[1] 0.00814316030659453
> f(x, n, p)
[1] NaN

==============================================

一つのやりかたは，掛算・割り算を log の足し算・引き算で計算し，結果の exp をとる。

g = function(x, n, p) {
  g0 = function(m) {
    ifelse(m <= 1, 0, sum(log(2:m))) # R の仕様からくる問題を避けるために関数を定義
  }
  return(exp(g0(n) - g0(x) - g0(n - x) + x * log(p) + (n - x) * log(1 - p)))
}

> dbinom(x, n, p)
[1] 0.00814316030659453
> g(x, n, p)
[1] 0.00814316030660582

このやり方では，有効数字は 10 桁ほどである。

==============================================

もう一つのやり方は，n の階乗は n+1 に対するΓ関数なので，それを利用して対数Γ関数を使う。

h = function(x, n, p) {
  return(exp(lgamma(n + 1) - lgamma(x + 1) - lgamma(n - x + 1) +
         　　x * log(p) + (n - x) * log(1 - p)))
}

> dbinom(x, n, p)
[1] 0.00814316030659453
> h(x, n, p)
[1] 0.00814316030654657

このやり方では，有効数字は 11 桁ほどである。

==============================================

R では nCx の対数は lchoose 関数で計算できるので，少し簡単に書ける。

v = function(x, n, p) {
  return(exp(lchoose(n, x) + x * log(p) + (n - x) * log(1 - p)))
}

> dbinom(x, n, p)
[1] 0.00814316030659453
> v(x, n, p)
[1] 0.00814316030660582

このやり方では，有効数字は 10 桁ほどである。

==============================================

nCx は整数値（約分すると分母は 1 になる）ということで，約分後の分子の log の和を求めて lchoose と同じ答えを目指す。

gcd = function(m, n) {
  while ((temp <- n %% m) != 0) {
    n <- m
    m <- temp
  }
  return(m)
}

lchoose2 = function(n, k) {
  k = min(k, n - k)
  if (k == 0) {
    return(0)
  }
  numerator = (n - k + 1):n
  denominator = 1:k
  for (den in seq_along(denominator)) {
    for (num in seq_along(numerator)) {
      GCD = gcd(denominator[den], numerator[num])
      if (GCD != 1) {
        denominator[den] = denominator[den] / GCD
        numerator[num] = numerator[num] / GCD
        if (denominator[den] == 1) {
          break
        }
      }
    }
    numerator = numerator[numerator != 1]
  }
  return(sum(log(numerator)))
}

w = function(x, n, p) {
  return(exp(lchoose2(n, x) + x * log(p) + (n - x) * log(1 - p)))
}

> dbinom(x, n, p)
[1] 0.00814316030659453
> w(x, n, p)
[1] 0.00814316030659842

このやり方では，有効数字は 12 桁に増加するが，計算時間が 5 秒近くかかる。
dbinom の結果の違いは，計算途中での計算精度が不足しているからである。

==============================================

R には多倍長演算をサポートするパッケージがいくつかあるが，そのうちの Rmfpr を使ってみる。
計算精度を適当（適切！）に指定して以下のような結果になる。

積・商を対数の和・差にする場合。

w2 = function(x, n, p) {
  library(Rmpfr)
  precBits = 70
  P = mpfr(p, precBits)
  lchoose = sum(log(mpfr((n-x+1):n, precBits))) - sum(log(mpfr(1:x, precBits)))
  as.numeric(exp(lchoose + x * log(P) + (n - x) * log(1 - P)))
}

積・商をそのまま計算する場合。

w3 = function(x, n, p) {
  library(Rmpfr)
  precBits = 70
  P = mpfr(p, precBits)
  choose = prod(mpfr((n-x+1):n, precBits)) / prod(mpfr(1:x, precBits))
  as.numeric(choose * (P ^ x) * (1 - P) ^ (n - x))
}

> dbinom(x, n, p)
[1] 0.00814316030659453
> w2(x, n, p)
[1] 0.00814316030659453
> w3(x, n, p)
[1] 0.00814316030659453

いずれでも， bdinom と同じ精度の結果が得られる。
計算時間は，積・商をそのまま計算する w3 の方が速く 0.14 秒ほどである（dbinom は瞬時）。

==============================================

二項分布 B(x, n, p) = nCx p^x (1-p)^(n-x) の計算の場合は dbinom があるので，以上のような事をやる必要はさらさらないが，「ものすごく大きな数とものすごく小さな数の積がほどほどの値になる場合」や「ものすごく大きな 2 つの数の商がほどほどの値になる場合」には，「積と商を対数の和と対数の差をとって結果の逆対数をとる」というのが必要になることもある。
そういうときには Rmfpr を使うのがもっとも効果的のようである。
... と，そういうこと。