SciPyのF分布ppfが抱える数値精度問題

概要

scipy.stats.f.ppf（F分布のパーセント点関数）が、特定の条件下で不正確な結果を返す可能性があるという問題についての報告です。具体的には、確率が1に非常に近い値（裾の領域）と大きな自由度の組み合わせにおいて、ppf関数の結果と、累積分布関数（cdf）を最適化して得られる結果との間に大きな乖離が見られます。

この問題は2021年頃に発見されましたが、最近のバージョンでも同様の現象が確認されたため、改めて報告するものです。

なお、本件は、既に、https://github.com/scipy/scipy/issues/20835 で報告されており、恐らく、scipy 1,17で対策はマージされると思います。

ENH: special: boostify F distribution functions

問題の詳細

scipy.stats.f.ppf に欠陥がある可能性が疑われます。最適化関数と累積分布関数（cdf）の組み合わせで得られる結果が、ppf関数の戻り値と一致しません。

再現コード

以下のコードは、sps.f.ppf を直接呼び出した場合と、spo.brentq を使って sps.f.cdf から逆関数を求めた場合の結果を比較します。

import scipy.stats as sps
import scipy.optimize as spo

# ppf関数を直接呼び出す
print(sps.f.ppf(0.99999995, 1, 50000, loc=0, scale=1))

# cdfと最適化関数brentqを使って逆関数を求める
print(spo.brentq(lambda x: sps.f.cdf(x, 1, 50000) - 0.99999995, 1, 1e10))

コードの出力

3333.8385803475894
29.725915444860586

出力結果からわかるように、2つの計算結果は大きく異なっています。

エラーメッセージ

この問題では、エラーメッセージは出力されません。コードは正常に終了しますが、得られる結果の信頼性に疑問があります。

原因の考察

このPythonコードと出力は、F分布のパーセント点関数（ppf）と、それを最適化関数（brentq）を使って求めるアプローチとの間の数値的な不安定性を示唆しています。

sps.f.ppf(0.99999995, 1, 50000)
これはF分布のパーセント点関数（ppf）を直接呼び出し、累積確率が0.99999995となるF値を求めます。この計算では、F値が 3333.83… という非常に大きな値になっています。
spo.brentq(lambda x: sps.f.cdf(x, 1, 50000) - 0.99999995, 1, 1e10)
こちらはscipy.optimizeモジュールのbrentq関数を使い、同じ問題を数値的に解いています。brentqは、与えられた関数の値が0になる点（根）を見つけるための堅牢なアルゴリズムです。ここでは、「sps.f.cdf(x, 1, 50000) - 0.99999995」という式が0になるx、つまり累積確率が0.99999995になるF値を探索しています。
この計算では、F値が 29.72… という、ppfの直接呼び出しとは全く異なる、より妥当と思われる値が得られます。

この乖離は、scipy.stats.f.ppfが特定の条件下で正確な答えを返さない可能性があることを示しています。

F分布のような統計分布の逆関数（ppf）を、極端な確率値（この例では1に非常に近い0.99999995）で計算する場合、数値的な精度の問題が発生しやすくなります。特に自由度が大きい場合、分布の裾（テール）の部分が非常に平坦になるため、計算過程でのわずかな丸め誤差などが、結果として得られるF値の大きな違いにつながることがあります。

一方、brentqは指定された範囲内で根を探索する、より頑健な数値計算アルゴリズムです。関数の振る舞いをより正確に追跡できるため、このようなケースではより信頼性の高い結果を生成すると考えられます。

問題の可視化

この問題は、以下のコードでscipy.special.fdtri（F分布の逆関数）の挙動をプロットすることで、より明確に可視化できます。自由度d2が50000と大きい場合に、確率pが0に近づく（つまり累積確率が1に近づく）領域で、関数の値が急激に変動していることがわかります。

import numpy as np
from scipy import special
import matplotlib.pyplot as plt

d1, d2, p = 1, 50000, np.logspace(-16, -1)
# fdtriは1-pを引数に取るため、pが0に近づくと累積確率が1に近づく
fdtri = special.fdtri(d1, d2, 1-p)

plt.semilogx(p, fdtri)
plt.xlabel("p (1 - Cumulative Probability)")
plt.ylabel("F-value (fdtri)")
plt.title("fdtri behavior for extreme probabilities (d1=1, d2=50000)")
plt.grid(True)
plt.show()

F分布の逆関数における極端な値での変動

このグラフは、F分布の逆関数（ppf）の計算が、極端な確率値や特定の自由度の組み合わせにおいて、依然として数値的な不安定さを抱えていることを示唆しています。この例では、自由度d2=50000という条件下で、裾の領域で値が不安定になっている様子が確認できます。

環境情報

以下は、この問題が確認された環境の情報です。

Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /usr/local/include
    lib directory: /usr/local/lib
    name: openblas
    openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=64
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.27
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /usr/local/include
    lib directory: /usr/local/lib
    name: openblas
    openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=64
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.27
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.12.0
Compilers:
  c:
    commands: cc
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  c++:
    commands: c++
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.10
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  pythran:
    include directory: ../../tmp/pip-build-env-mnl4e8vy/overlay/lib/python3.10/site-packages/pythran
    version: 0.15.0
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
Python Information:
  path: /opt/python/cp310-cp310/bin/python
  version: '3.10'

概要#

問題の詳細#

再現コード#

コードの出力#

エラーメッセージ#

原因の考察#

問題の可視化#

環境情報#

概要