本当に正規分布の正規四分位範囲が標準偏差と一致するのか SymPy になったので確かめてみた

第 3 四分位数と第 1 四分位数の差によって導かれる四分位範囲 (IQR: Interquartile Range) は、標準偏差と同様に統計的ばらつきの指標としてよく用いられ、標準偏差と比較して外れ値などの外乱に対して頑強 (ロバスト) であるために有用であるといわれる。 統計の知識がない一般の人でも、標準偏差より計算の理屈がわかりやすい¹ため、そういう意味でも有用といえる。

さらに、四分位範囲を「標準正規分布における四分位範囲」で除した値を正規四分位範囲 (NIQR: Normalized Interquartile Range) といい、正規分布の正規四分位範囲は標準偏差に一致するので、これまで標準偏差を利用した経験が多い人には、それに近いスケール感で統計的ばらつきを把握できる正規四分位範囲は輪をかけて有用そうだ。

しかし、この「正規分布の正規四分位範囲は標準偏差に一致する」という話は、単に Wikipedia から仕入れた情報に過ぎない。 不特定多数によって編集される Wikipedia が意外と信用ならないことはご存じの方も多いだろう。

そんなわけで、この「正規分布の正規四分位範囲は標準偏差に一致する」という命題が真か偽か、Python の記号計算ライブラリである SymPy を使いつつ確認してみよう。 Python が計算機代数システムになるとはおっちゃんびっくりだよ²。

確率分布におけるクォンタイル

一般的に四分位数を気にするとき、その対象は具体的な数値があることが多い (たとえば小学生の身長のデータとか、テスト成績のデータとか)。 こういった変量統計の分野ではともかく、正規分布という連続的な確率分布において四分位数 (および任意のクォンタイル) を求めるということはできるのだろうか?

この疑問に答えるにはそもそもクォンタイルとはなんだったのかを思いだす必要がある。 第 1 四分位数 (すなわち 0.25 クォンタイル) $Q_1$ は下位 25 パーセントのラインとなるわけだが、これはつまりすべてのデータを $Q_1$ 以下と $Q_1$ 以上に³分けたとき、それが $\frac{1}{4} : \frac{3}{4}$ という比率になるということになる。 そして、これこそがクォンタイルの定義である。

定義

すべての分布を $Q$ 以下と $Q$ 以上に分割した比率 $q : 1 - q$ をもって、値 $Q$ を $q$ クォンタイルとよぶ。

このように考えれば、確率分布におけるクォンタイルをどのように求めるかも予想がつくのではないだろうか。 なぜなら、確率密度函数とは以下のような性質を満たす函数だからである。

確率変数 $X$ についての確率密度函数 $f_X(x)$ は、以下の性質を持つ。

• 確率変数 $X$ が $X_1$ から $X_2$ までの値をとる確率は $\int_{X_1}^{X_2} f_X(x) dx$ (ただし、$X_1 < X_2$)
• $\int_{-\infty}^{\infty} f_X(x) dx = 1$

つまり、

\begin{eqnarray}
    \int_{-\infty}^Q f_X(x) dx &\geq& q \\
    \int_Q^{\infty} f_X(x) dx &\geq& 1 - q \\
\end{eqnarray}

を満たすような値 $Q$ を確率分布の $q$ クォンタイルとよべる。 以下では簡単のため
$$
\int_{-\infty}^Q f_X(x) dx = q
$$
を満たすような値 $Q$ のみを考えることにする。

正規分布における四分位数と四分位範囲が満たすべき性質

ここで、第 1 四分位数 $Q_1$ と第 3 四分位数 $Q_3$ を考えると、

\begin{eqnarray}
    \int_{-\infty}^{Q_1} f(x) dx &=& \frac{1}{4} \\
    \int_{-\infty}^{Q_3} f(x) dx &=& \frac{3}{4} \\
\end{eqnarray}

であるから、
$$
\int_{-\infty}^{Q_3} f(x) dx - \int_{-\infty}^{Q_1} f(x) dx = \int_{Q_1}^{Q_3} f(x) dx = \frac{1}{2}
$$
となる。 ここで確率分布が正規分布であった場合、第 1 四分位数と第 3 四分位数は平均 $\mu$ に対して対称になるため、$Q_1 = \mu - c$ および $Q_3 = \mu + c$  とおいて、
$$
\int_{\mu - c}^{\mu + c} \varphi_{\mu, \sigma^2}(x) dx = \frac{1}{2}
$$
である。 さらに、四分位範囲 $\mathrm{IQR}$ は $\mathrm{IQR} = Q_3 - Q_1 = (\mu + c) - (\mu - c) = 2c$ であるから、これを $c$ で解いた $c = \frac{\mathrm{IQR}}{2}$ を代入して、
$$
\int_{\mu - \frac{\mathrm{IQR}}{2}}^{\mu + \frac{\mathrm{IQR}}{2}} \varphi_{\mu, \sigma^2}(x) dx = \frac{1}{2}
$$
となる。 この方程式を解くことによって正規分布の四分位範囲を求めることができる。

SymPy による記号計算

ここまで来れたので、さっそく SymPy を使用して計算していく。

なお、以下では Jupyter 上での作業を想定している。 Jupyter にはセルの最後の行の式を評価した結果を出力する機能が存在するが、Jupyter + SymPy ではその出力を MathJax を利用した美麗な数式で表示することができる。

from sympy import *
init_printing()

このようにinit_printing関数を呼び出すことで、Jupyter 出力 (内部的にはIPython.core.display.display関数にあたるようだ) で MathJax を利用可能になる。 美麗な数式は精神に正の影響を与えることが知られている^{[独自研究]}ので、積極的に利用していこう。

ちなみに、import *は名前空間を気にすることなく無造作にインポートしてしまうので、「そんなことしたら pep8 に怒られちゃうだろ!」と教わった方が多いかとは思うが、SymPy は単体で完結する計算機代数システムを志向しているため、このような「独裁的」な手法が慣例化していると思われる。 この記事でもそれに従う。 長いものに巻かれるのが私の生き様だ。

まず、正規分布の確率密度函数を記述していく。 正規分布の確率密度函数は以下である。
$$
\varphi_{\mu,\sigma^2}(x) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left({-\frac{(x - \mu)^2}{2\sigma^2}}\right)
$$
これを、必要な記号を定義しながらそのまま記述していく。 ここで、$\sigma$ は非負なのでpositive=Trueとしておく。

x = Symbol('x')
mu = Symbol('\\mu')
sigma = Symbol('\\sigma', positive=True)

phi = (1 / sqrt(2 * pi * sigma**2)) * exp(-((x - mu)**2 / (2 * sigma**2)))
phi

$$
\frac{\sqrt{2} e^{- \frac{\left(- \mu + x\right)^{2}}{2 \sigma^{2}}}}{2 \sqrt{\pi} \sigma}
$$

SymPy 側の流儀に従って整理した状態で出力されるため入力した式とは一見形が異なるが、検算してみると同一の式になっているのがわかる。

続いてこの式を積分し、先程挙げた方程式を記述する。

IQR = Symbol('\\mathrm{IQR}', positive=True)

eq_iqr = integrate(phi, (x,mu - IQR / 2,mu + IQR / 2,)) - Rational(1, 2)
eq_iqr

$$
\operatorname{erf}{\left(\frac{\sqrt{2} \mathrm{IQR}}{4 \sigma} \right)} - \frac{1}{2}
$$

積分した結果登場した $\mathrm{erf}$ は、誤差函数とよばれるものらしい。 この辺、筆者は数学があまり得意ではないので詳しくはわからないが。

ここで、eq_iqr = integrate(phi, (x,mu - IQR / 2,mu + IQR / 2,)) - 1 / 2としてしまうと、1 / 2の部分がintとintの商となるため Python 組み込み型の演算規則に従ってfloat型となってしまう。 floatが入り込むと有限精度が邪魔をしてきれいな記号計算を行うことができなくなってしまう。 int同士の商は/演算子を使用せず、SymPy の有理数オブジェクトRationalとして運用すること。

この方程式を $\mathrm{IQR}$ について解く。 方程式を解く SymPy の関数solveは、複数の解が求まることを考慮してリストで返却されるが、今回求める解は一意のため、先頭の要素のみを取り出す。

IQR_N = solve(eq_iqr, IQR)[0]
IQR_N

$$
2 \sqrt{2} \sigma \operatorname{erfinv}{\left(\frac{1}{2} \right)}
$$

これが正規分布の四分位範囲である。 $\mathrm{erfinv}$ は誤差函数の逆函数だが筆者は略。

この時点で $\sigma$ が全体の係数なので……などと考えていけばおのずと結論は出そうだが、ここは最後まで計算していこう。

標準正規分布における四分位範囲 $\mathrm{IQR}_{\mathcal{N}(0, 1)}$ を求めるために、これに $\mu = 0, \sigma = 1$ を代入する。 $\mu$ に関してはすでに消えているが。

IQR_N_0_1 = IQR_N.subs([(mu,0,), (sigma,1,),])
IQR_N_0_1

$$
2 \sqrt{2} \operatorname{erfinv}{\left(\frac{1}{2} \right)}
$$

ここで正規四分位範囲 $\mathrm{NIQR}$ について考える。 $\mathrm{NIQR} = \frac{\mathrm{IQR}}{\mathrm{IQR}{\mathcal{N}(0, 1)}}$ であるから、これを $\mathrm{IQR}$ について解いた $\mathrm{IQR} = \mathrm{NIQR} \cdot \mathrm{IQR}{\mathcal{N}(0, 1)}$ を先の方程式に代入する。 あーもうめちゃくちゃだよ。 Qiita くん、パーサはちゃんと作ろう!

$$\mathrm{NIQR} = \frac{\mathrm{IQR}}{\mathrm{IQR}_{\mathcal{N}(0, 1)}}$$

であるから、これを $\mathrm{IQR}$ について解いた $\mathrm{IQR} = \mathrm{NIQR} \cdot \mathrm{IQR}_{\mathcal{N}(0, 1)}$ を先の方程式に代入する。

NIQR = Symbol('\\mathrm{NIQR}', positive=True)

eq_niqr = eq_iqr.subs(IQR, NIQR * IQR_N_0_1)
eq_niqr

$$
\operatorname{erf}{\left(\frac{\mathrm{NIQR} \operatorname{erfinv}{\left(\frac{1}{2} \right)}}{\sigma} \right)} - \frac{1}{2}
$$

最後に、この方程式を $\mathrm{NIQR}$ について解く。

NIQR_N = solve(eq_niqr, NIQR)[0]
NIQR_N

$$
\sigma
$$

見事、正規分布の正規四分位範囲が標準偏差に等しいことが証明できた。

おまけ

SymPy は式を任意精度で計算することもできる。 前回の記事で Wikipedia から引っ張ってきた値で決め打ちしていた「標準正規分布における四分位範囲」を 500 桁まで計算してみよう。

IQR_N_0_1.evalf(500)

$$
1.3489795003921634864044540290826143707738088300997237913241875771896973585648861821908900894803355691470944268025862559823609241772506202548596796690440429997246378194772736295225048764221238103767277303955351652277519615363293040517301300404141772969519027307835818042849164133990052372200109437856077024718496431614331559376641877640706752297852554106153468551705697059662463626713575464038601510702286202536700805101035448972471611842629431350954384053091400476131362345086159068534983331089117546
$$

つよい (確信)。

SymPy はなかなか使い勝手がいいライブラリだと思う。