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

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

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

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

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

[^definition]: あくまでも主観だが、標準偏差が「各データの平均との差を 2 乗して、さらにそれらの平均をとり、その平方根をとったもの」という説明になるのに対し、四分位範囲は「全体の上位 25 パーセントラインと下位 25 パーセントラインの差」となり、後者のほうが理解しやすそうである。
[^sympy]: まだ私も SymPy を遊び尽くせていないのでアレなのだが、いずれ SymPy を紹介する記事も書いてみたい。

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

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

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

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

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

> 確率変数 $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$

つまり、

```math
\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$ のみを考えることにする。

[^than]: ここでいう「以下」とか「以上」は、かなり曖昧で不正確で語弊がありまずい単語チョイスではあるのだが、直感に訴える意味を持ってこのようにしている。

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

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

```math
\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 を利用した美麗な数式で表示する**ことができる。

```py:In[]
from sympy import *
init_printing()
```

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

ちなみに、`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`としておく。

```py:In[]
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 側の流儀に従って整理した状態で出力されるため**入力した式とは一見形が異なる**が、検算してみると**同一の式**になっているのがわかる。

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

```py:In[]
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`は、複数の解が求まることを考慮して**リストで返却される**が、今回求める解は一意のため、**先頭の要素のみを取り出す**。

```py:In[]
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$ に関してはすでに消えているが。

```py:In[]
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 くん、パーサはちゃんと作ろう!:anger:

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

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


```py:In[]
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}$ について解く。

```py:In[]
NIQR_N = solve(eq_niqr, NIQR)[0]
NIQR_N
```

$$
\sigma
$$

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

## おまけ

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

```py:In[]
IQR_N_0_1.evalf(500)
```

$$
1.3489795003921634864044540290826143707738088300997237913241875771896973585648861821908900894803355691470944268025862559823609241772506202548596796690440429997246378194772736295225048764221238103767277303955351652277519615363293040517301300404141772969519027307835818042849164133990052372200109437856077024718496431614331559376641877640706752297852554106153468551705697059662463626713575464038601510702286202536700805101035448972471611842629431350954384053091400476131362345086159068534983331089117546
$$

**つよい (確信)**。

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