【Python】実装しながら理解するクォンタイル(パーセンタイル)

## クォンタイルってなんだよ

**クォンタイル**とは、**統計的代表値**のひとつで、日本語では**分位数**ともよびます。 $q$ クォンタイルは**分布を $q : 1 - q$ に分割する**値のことです。

:thinking::thinking::thinking::thinking::thinking::thinking::thinking::thinking::thinking::thinking::thinking::thinking:???????

えー、**パーセンタイル**という単語を知っている方はそっちで置き換えていただいて構いません。 **100 パーセンタイルが 1 クォンタイル**なだけです。

**中央値**、という単語を知っている方は、中央値が **0.5 クォンタイル**および **50 パーセンタイル**であるといえば理解しやすいかもしれません。

中央値という言葉に聞き覚えがない方もいるかとは思いますが、簡単に言えば **100 人のランキングで 50 位が中央値**です(雑)。 成績上位 10 パーセントとは 0.9 クォンタイル以上のスコアを取ったということですね。 ね、簡単でしょう?

## 実装すんだよ

クォンタイルについて理解できた(?)ところで早速実装していきましょう。

ま、とりあえずは当然のように利用する`math`パッケージをインポートしておきます。

```py:In[]
import math
```

今回はサンプルとして以下のような数値のクォンタイルを求めていきましょう。 テスト成績のイメージです。

```py:In[]
data = [78., 42., 100., 72., 89., 51., 75., 63., 82.,]
```

いきなりクォンタイルすべてを考えても:thinking:になるので、代表的なクォンタイルである中央値をヒントにしていきましょう。

まず、中央値を求める際には**順序どおり並べられている**必要があります。 そののち、それらの中で中央の値を取ります。 今回の要素数は 9 で、中央は (プログラムらしく 0 始まりで) 4 番目です。 これは $\frac{9 - 1}{2}$ と表せます。 他の要素数でも同様です。

先程の式は $\frac{1}{2}(9 - 1)$ と書けます。  中央値は前述の通り 0.5 クォンタイルなわけですが、**0 クォンタイルが最小値すなわち 0 番目**、**1 クォンタイルが最大値すなわち 8 番目**であることと合わせると、$n$ 個の要素の $q$ クォンタイルは、**昇順ソートを行ったうえで $q(n - 1)$ 番目の要素を取得することで求められる**、と一般化できそうです。

```py:In[]
def quantile(data, q):
    sorted_data = sorted(data)
    return sorted_data[q * (len(data) - 1)]
```

さぁ、早速サンプルデータの中央値を求めてみましょう!

```py:In[]
quantile(data, 0.5)
```

```:Out[]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-4-cef44d6c287e> in <module>
----> 1 quantile(data, 0.5)

<ipython-input-3-f523a2f1f0d2> in quantile(data, q)
      1 def quantile(data, q):
      2     sorted_data = sorted(data)
----> 3     return sorted_data[q * (len(data) - 1)]

TypeError: list indices must be integers or slices, not float
```

……。

**float との積をとった時点で float 型になってしまいます**ので、添え字として利用できなくなってしまいました。

では、int 型に変換すればいいのかというと、そう簡単な話でもありません。 **要素の数が偶数だった場合**、たとえば **3.5 番目**の要素を要求されます。 これに対して 3 番目の要素を渡してしまっては不誠実というものです。

中央値の求め方を知っている方はどのようにするか思い出してみましょう。 知らない方は想像してみてください。

そうです。 **3 番目の要素と 4 番目の要素の平均**をとるのです。

3.5 に対して 3 は「**それ以下の最大の整数**」であり、**床函数**で求められます。 3.5 に対して 4 は「**それ以上の最小の整数**」なので**天井函数**で求められます。 これらは``math``パッケージにあります。 これらは **int 型の値を返却する**のでそのまま添え字として利用可能です。

**「添え字」を非負な実数全体に拡張した関数**を実装していきましょう。 この関数は整数かどうかを判定して内部で場合分けします。 クォンタイルを求める関数では添え字を直接呼び出す代わりに、この「**拡張添え字関数**」に添え字を渡します。

```py:In[]
def complemented(data):
    def f(q):
        if isinstance(q, int) or q.is_integer():
            qi = int(q)
            return data[qi]
        else:
            qc = math.ceil(q)
            qf = math.floor(q)
            return (data[qf] + data[qc]) / 2
    return f

def quantile(data, q):
    sorted_data = complemented(sorted(data))
    return sorted_data(q * (len(data) - 1))
```

```py:In[]
quantile(data, 0.5)
```

```:Out[]
75.0
```

良い感じに出ていますね。 要素数が偶数でも、

```py:In[]
quantile(data[:-1], 0.5)
```

```:Out[]
73.5
```

良い感じです。

**良いわけありません**。

我々が求めたいのは**任意のクォンタイル**であって、中央値だけではありません。 クォンタイルのとり方によっては 2.5 番目とか 3.5 番目とか 7.5 番目とかだけではなく、**3.1 番目**とか **6.7 番目**とか **1.14514 番目**とかをとる必要が出てきます。 3.1 番目に対しても 3 番目と 4 番目の平均を渡すべきでしょうか? ほぼ 3 なのに?

というか我々は平均という言葉に囚われすぎです。 あんなもん所詮足して割っているだけです。
$$
\frac{x_3 + x_4}{2}
$$
これは以下のように書き直せますね。
$$
\frac{1}{2}(x_3 + x_4)
$$
なにが $\frac{1}{2}$ だ小数にするぞ。
$$
0.5(x_3 + x_4)
$$
展開できます。
$$
0.5x_3 + 0.5x_4
$$
ここまでくるとなんとなく見えてきますね。 たとえば 3.1 番目は 3 に近いので以下のようにすればいいのではないでしょうか。
$$
0.9x_3 + 0.1x_4
$$
いいアイディアだとは思いますが、0.9 や 0.1 といった「重み」係数はどうやって求めましょう。 ここで重要なヒントは、$3.1 + 0.9 = 4$ であり、$3.1 - 0.1 = 3$ であることです。 逆から見れば、$4 - 3.1 = 0.9$ であり、$3.1 - 3 = 0.1$ であるということです。
$$
(4 - 3.1)x_3 + (3.1 - 3)x_4
$$
これを先程の床函数と天井函数を使って一般化すると、**整数でない $i$ 番目の値**を以下のようにして補完できます。
$$
(\lceil i \rceil - i)x_{\lfloor i \rfloor} + (i - \lfloor i \rfloor)x_{\lceil i \rceil}
$$
さあ、もう一度拡張添え字関数を書き直しましょう。

```py:In[]
def complemented(data):
    def f(q):
        if isinstance(q, int) or q.is_integer():
            qi = int(q)
            return data[qi]
        else:
            qc = math.ceil(q)
            qf = math.floor(q)
            return (qc - q) * data[qf] + (q - qf) * data[qc]
    return f
```

適当なクォンタイルを求めてみます。

```py:In[]
quantile(data, 0.333)
```

```:Out[]
68.976
```

よくわかりませんが、なかなかよさそうです。

ついでに**複数のクォンタイルを同時に求められる**ように拡張しておきましょう。

```py:In[]
def quantile(data, q):
    sorted_data = complemented(sorted(data))
    if hasattr(q, "__iter__"):
        return [sorted_data(i * (len(data) - 1)) for i in q]
    else:
        return sorted_data(q * (len(data) - 1))
```

```py:In[]
quantile(data, [0, 0.25, 0.5, 0.75, 1,])
```

```:Out[]
[42.0, 63.0, 75.0, 82.0, 100.0]
```

**パーセンタイル版**も用意しておきます。

```py:In[]
def percentile(data, p):
    if hasattr(p, "__iter__"):
        return quantile(data, [i / 100 for i in p])
    else:
        return quantile(data, p / 100)
```

```py:In[]
percentile(data, [0, 25, 50, 75, 100,])
```

```:Out[]
[42.0, 63.0, 75.0, 82.0, 100.0]
```

## 検算(?)

みんな大好き **NumPy** くんにはパーセンタイルを計算できる関数が用意されているので、比較してみましょう。

```py:In[]
import numpy as np

p = [0, 10, 25, 50, 75, 90, 100,]
print(F'NumPy:    {list(np.percentile(data, p))}')
print(F'Original: {percentile(data, p)}')
```

```:Out[]
NumPy:    [42.0, 49.2, 63.0, 75.0, 82.0, 91.2, 100.0]
Original: [42.0, 49.2, 63.0, 75.0, 82.0, 91.2, 100.0]
```

ピッタリ**一致**していますね。

## おまけ

### 統計的ばらつきの時間だオラァ!

せっかくクォンタイルが計算できるようになったので、**IQR** を計算できるようにしてみましょう。 IQR とは **Interquartile Range** の略で、**0.75 クォンタイルと 0.25 クォンタイルの差**です。 難しいことはありませんね。 よく見ると Qua**n**tile ではなく Qua**r**tile です。 日本語では**四分位範囲**といい、クォーターとかそっち系列なのがわかります。 おい強調構文のパーサおかしいぞ。 聞いているのか Qiita。

```py:In[]
def iqr(data):
    q = quantile(data, [0.25, 0.75,])
    return q[1] - q[0]
```

```py:In[]
iqr(data)
```

```:Out[]
19.0
```

この指標をなにに使えるかと言うと、標準偏差のように**値のばらつきを示す**のによく利用されます。 標準偏差が平均と似た仕組みで求められているのに対し、IQR は中央値と同様にクォンタイルから求められているため、**外れ値などに強い**、ナウい言い方をすると**ロバスト性**があります。

しかし、標準偏差となんとなく尺度が異なるのでなんか使いづらいイメージがないでもないです。 一方で **NIQR** という指標は**標準正規分布のそれがちょうど 1 になります**。 **Normalized Interquartile Range** です。 日本語では**正規四分位範囲**です。 ~~長過ぎる。~~

```py:In[]
def niqr(data):
    return iqr(data) / 1.3490
```

```py:In[]
niqr(data)
```

```:Out[]
14.084507042253522
```

`1.3490`とは**標準正規分布のおおよその IQR** です。 確率分布にクォンタイル?と思うかもしれませんが、こちらについては~~いずれ別の記事で触れてみたいなあと思います。~~ **2019-10-30 追記**: [書きました](https://qiita.com/BlueRayi/items/ce0f9b3cbb57bfed855f)。

**正規分布では NIQR と標準偏差が等しい**ので、標準偏差と同じ感覚で使えるんじゃないかなと思います。 統計全然やったことないのでよくわかりませんが。

### クォータイルのリーダー「お前たちの場合分けって、醜くないか?」

:sushi:「そうか? 俺には整数を場合分けせずにクォンタイルを求める力がある」

床函数 $\lfloor x \rfloor := \max\\{n \in \mathbb{Z} | n \leq x\\}$ は「$x$ **以下**の最大の整数」を返す函数ですから、整数を渡すと引数をそのまま返却します。 ですから、先ほど挙げた式を整数の添え字 (たとえば $i = 3$ など) に適用するとおかしなことになります。

```math
\begin{eqnarray}(\lceil i \rceil - i)x_{\lfloor i \rfloor} + (i - \lfloor i \rfloor)x_{\lceil i \rceil} &=& (3 - 3)x_3 + (3 - 3)x_3\\&=& 0x_3 + 0x_3 \\&=& 0\end{eqnarray}
```

これが拡張添え字関数で場合分けをしている理由なわけですが、なんというかあちこち凸凹で**美しくありません**。 ~~我が魔王はちょっと座っててください。~~

ですが、「$x$ **未満**の最大の整数」を返す函数 $\lfloor x \rfloor^{\ast} := \max\\{n \in \mathbb{Z} | n < x\\}$ を考え、これを床函数の代わりに置くと、以下に示すように**整数の添え字でも適用できる**ようになります。

```math
\begin{eqnarray}(\lceil i \rceil - i)x_{\lfloor i \rfloor^{\ast}} + (i - \lfloor i \rfloor^{\ast})x_{\lceil i \rceil} &=& (3 - 3)x_2 + (3 - 2)x_3\\&=& 0x_2 + 1x_3 \\&=& x_3\end{eqnarray}
```

この函数は $\lfloor x \rfloor^{\ast} < x \leq \lfloor x \rfloor^{\ast} + 1$ という関係を満たすわけですが、天井函数が $\lceil x \rceil - 1 < x \leq \lceil x \rceil$ という関係を満たすことを考えるとわかるように、**$\lfloor x \rfloor^{\ast} = \lceil x \rceil - 1$ となります**。

これを利用することで、場合分けをせずに拡張添え字関数を定義することができます。

```py:In[]
def complemented(data):
    def f(q):
        qc = math.ceil(q)
        qf = qc - 1
        return (qc - q) * data[qf] + (q - qf) * data[qc]
    return f
```

## あとがきポエム

やっぱりこうやって車輪の再発明をしていくのは、本を読んでわかったつもりになるよりも数十倍身につきやすい気がします。

今回の記事では数式をいきなり提示してそれをコードに落とし込むだけではなく、**その数式自体を直感で納得しやすいところから順に積み上げていく**ようにしたつもりです。 数式は天から降ってくるのではなく己のうちから湧き上がってくるものなのだということを感じていただければ幸いです(?)。