桁落ち
桁落ちという言葉をご存知だろうか. 絶対値が非常に近い数同士の加減算で有効桁数が少なくなってしまうことを桁落ちいう. 例えば, 円周率とネイピア数の差について, 有効数字を 3 桁として計算すると次のようになる.
| 3 | . | 1 | 4 | ? | ? | |
| – | 2 | . | 7 | 2 | ? | ? |
| ― | ― | ― | ― | ― | ― | ― |
| 0 | . | 4 | 2 | ? | ? |
この例では, 有効桁数が 3 桁から 2 桁に減っている. このようなことは, プログラム上の計算でも起こりうる話だ. ここまでは, 以前に仮平均を用いた平均と分散の計算でも触れた話だ.
さて, 桁落ちを回避する方法に一般論は存在しない. むしろ, 回避できるケースが珍しいというのだ. 桁落ちを回避できるケースとして, 二次方程式の解を求める方法を考える.
二次方程式 ax2+bx+c=0 (a≠0) における解は, (-b+√(b2-4ac))/(2a) と (-b-√(b2-4ac))/(2a) の 2 つである
(平方根の底が 0 に等しいならば重解となる). 問題は, 4ac の絶対値が b2 と比べて非常に小さい値になるときである. この場合, b>0 のときは-b+√(b2-4ac) の精度がかなり落ちてしまうので, 結果として, 得られる解の精度も落ちてしまう. これを回避する方法として, 次の 2 つの考え方が知られている (計算の内容はほぼ同じ).
- もう一方の解 (-b-√(b2-4ac))/(2a) の精度はそれほど悪くならないので, 方程式の解と係数の関係を用いた次の計算を行う. α=(-b-√(b2-4ac))/(2a) とし, もう一方の解を β とすると c/a=αβ となることを利用し, β=c/(αa) を計算する.
- 分子を有理化して -2c/(b+√(b2-4ac)) を計算する.
このように, 非常に近い値同士の差を計算しないように式変形することが, 桁落ちを回避する根本的な手段となる.
さて, ここからが本題. 同一の二項分布 (母比率 π) に従う母集団から独立に n 個のサンプルを取り出したとき, そのすべてが陰性であったとする. このときの母比率 π の100(1-α)%信頼区間は [0, 1-α1/n] となる (大多数の文献では, サンプルに陽性のものと陰性のものが少なくとも 1 つずつある場合と処理を分岐しなくてもすむように, [0, 1-(α/2)1/n] としている. 本記事の主題ではないが, 指摘されそうな点であるので紹介しておく).
ここまでの前振りから想像できるかもしれないが, 非常に大きな n においてはα1/n の値は 1 に限りなく近い値となる. この値と 1 の差をとると, 桁落ちが発生してしまうため, コンピュータでは 0 に等しいという結果になってしまう. この計算において桁落ちを回避したいというのが本題である. 累乗が絡む計算において, 底はネイピアの数にしておいたほうが都合がいいことが多いと思うので次のように変形しておこう.
1-α1/n = 1-exp(ln(α))1/n = 1-exp(ln(α)/n)
t=-ln(α)/n とすれば, t がごく微小な正の数のときを考えればよいことになる.
以上を整理すると, 非常に小さな正の数 t に対して, 1-exp(-t) をなるべく高い精度で計算する方法があればぜひ教えてほしい.
