正規分布の分位点関数(パーセント点関数)の近似
背景
正規分布の分位点(パーセント点)は数値計算ライブラリの分位点関数を使えば簡単に求められる。 例えば、python だと、scipy.stats.norm の、ppt だ。 ところが、SQL には標準でないようである。 自分で関数を作ってもいいが、ここはSQLでかける近似式をつかってもいいかもしれない。
近似式を調べる
いろいろ調べたが、Hastings 1955 は有名だが、その後多くの近似式があるようである。 全部を調べるは大変で、サーベイっぽい論文などもあるのだが、最終的には、多く場合分けをして、多く次数を増やしていけば、近似精度がよくなるという話。 これはドメイン依存で、すなわち、どのような状況でどの程度近似精度が重要なのかと、SQLなどで記述する場合に複雑になっていってもいいのかというトレードオフの問題がある。 あと、近似式だから、正解はなく、多くの大量に考え出されているので、それら全てを評価するのもコストがかかる。
というわけで、私が探した中で、一番シンプルな式の山内 1965から検証する。
正規分布の分位点関数の近似式 山内 1965
戸田 英雄, 小野令美 「正規分布のパーセント点を求めるミニマックス近似」(1993) によると式は
らしい。 の部分がまとめられてしまっているが係数も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:分位という言葉は正しいかわからないが、ここでは分位点を調べるための割合、確率をいう