数値計算の定石，そして，R の定石

NHKスペシャルを題材に、Ｒコードの最適化を考える－その２
http://markovchainmontecarlo.hatenablog.com/entry/2013/12/27/000000

だけど。ネタ振りしているのだとは思うのだけど，そうでもないのかなとも思う。

試しに5,000個中250個丁度の貸し倒れが起きる確率を上記算式に当てはめてみる。
> n <- 5000
> k <- 250
> p <- 0.2
> prod(1:n) / prod(1:k) / prod(1:(n-k)) * p^k * (1-p)^(n-k)
[1] NaN

え？NaN？
非数 (Not a Number) が出てきたぞ？

全然ダメじゃん。

ということで、実際に5,000このサイコロを何度も振ることで確率論を求める方法にシフトしてみます。
その中で最適化を考える事にしてみましょう。

ではなく，この場合は，
> exp(lchoose(n, k)+k*log(p)+(n-k)*log(1-p))
[1] 2.619655e-206
とするべきですね。コンピュータによる数値演算の定石です。
そもそも，そんな小細工しなくても，R はちゃんと答えを出してくれます。

> dbinom(k, n, p) # たったこれだけ
[1] 2.619655e-206