中野智文のブログ

データ・マエショリストのメモ

正規分布の分位点関数(パーセント点関数)の近似

背景

正規分布の分位点(パーセント点)は数値計算ライブラリの分位点関数を使えば簡単に求められる。 例えば、python だと、scipy.stats.norm の、ppt だ。 ところが、SQL には標準でないようである。 自分で関数を作ってもいいが、ここはSQLでかける近似式をつかってもいいかもしれない。

近似式を調べる

いろいろ調べたが、Hastings 1955 は有名だが、その後多くの近似式があるようである。 全部を調べるは大変で、サーベイっぽい論文などもあるのだが、最終的には、多く場合分けをして、多く次数を増やしていけば、近似精度がよくなるという話。 これはドメイン依存で、すなわち、どのような状況でどの程度近似精度が重要なのかと、SQLなどで記述する場合に複雑になっていってもいいのかというトレードオフの問題がある。 あと、近似式だから、正解はなく、多くの大量に考え出されているので、それら全てを評価するのもコストがかかる。

というわけで、私が探した中で、一番シンプルな式の山内 1965から検証する。

正規分布の分位点関数の近似式 山内 1965

戸田 英雄, 小野令美 「正規分布のパーセント点を求めるミニマックス近似」(1993) によると式は


\begin{align}
 \left[ y \left( 2.0611786 - \frac{5.7262204}{y + 11.640595} \right) \right] ^{ \frac{1}{2} }, \\
y = - \log 4P(1-P)
\end{align}

らしい。  y の部分がまとめられてしまっているが係数も3つで乗数もなく、入力の間違いの確認も容易である。 二乗根と対数ならSQLもサポートしている。

python での近似性能の検証

多くの場合分位点そのものに興味はない。 だから分位の近似の評価ではなく、近似された分位点から逆に計算された分位*1が元の分位と比べてどうだったかを調べる。

import numpy as np

def yamauchi1965(q):
    y = -np.log(4*q*(1-q))
    return np.sqrt(y*(2.0611786-5.7262204/(y+11.640595)))

def yamauchi1965_sql(q):
    return np.sqrt(-np.log(4*q*(1-q))*(2.0611786-5.7262204/(-np.log(4*q*(1-q))+11.640595)))
    # sqrt(-log(4*Q*(1-Q))*(2.0611786-5.7262204/(log(4*Q*(1-Q))+11.640595)))

# 検証用
from scipy.stats import norm
for i in [1, 10, 20, 50, 100, 1000, 10000, 100000]:
    q = 1.0-1.0/2/i
    z = norm.ppf(q)
    z1 = yamauchi1965(q) # これが近似の分位点z1
    q1 = norm.cdf(z1)    # 近似の分位点z1から、元の分位を逆計算q1
    r = 1-(1-q1)/(1-q)   # 分位(1-q)と近似から逆計算した分位(1-q1)の割合の違い(割合)からさらに1を引いて相対的な差を求める
    print("q: %f, z: %f, z1: %f, q1: %f, r: %f" % (q, z, z1, q1, r))

実行結果

q: 0.500000, z: 0.000000, z1: -0.000000, q1: 0.500000, r: 0.000000
q: 0.950000, z: 1.644854, z1: 1.645636, q1: 0.950081, r: 0.001613
q: 0.975000, z: 1.959964, z1: 1.960594, q1: 0.975037, r: 0.001471
q: 0.990000, z: 2.326348, z1: 2.326358, q1: 0.990000, r: 0.000027
q: 0.995000, z: 2.575829, z1: 2.575264, q1: 0.994992, r: -0.001635
q: 0.999500, z: 3.290527, z1: 3.288966, q1: 0.999497, r: -0.005561
q: 0.999950, z: 3.890592, z1: 3.890511, q1: 0.999950, r: -0.000335
q: 0.999995, z: 4.417173, z1: 4.420750, q1: 0.999995, r: 0.016409

0.5 がエラーになると思っていたが、特に問題はなさそう。 元の分位に戻したときの相対的な差は多くの場合1%以下に抑えられている。 個人的はこれで満足である。

ただし、標準正規分布表と比べたときに、一番下の段の方では、違いが出てくる。

まとめ

山内 1965

sqrt(-log(4*Q*(1-Q))*(2.0611786-5.7262204/(log(4*Q*(1-Q))+11.640595)))

ただし、Q は分位(割合・確率)。

ちなみに、mysql では、これらの整数を4 の代わりに 4.000000 としてやらないと、計算してくれないことがあるので注意

*1:分位という言葉は正しいかわからないが、ここでは分位点を調べるための割合、確率をいう