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

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

???????

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

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

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

実装すんだよ

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

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

In[]

import math

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

In[]

data = [78., 42., 100., 72., 89., 51., 75., 63., 82.,]

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

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

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

In[]

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

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

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 型の値を返却するのでそのまま添え字として利用可能です。

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

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))

In[]

quantile(data, 0.5)

Out[]

75.0

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

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}
$$
さあ、もう一度拡張添え字関数を書き直しましょう。

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

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

In[]

quantile(data, 0.333)

Out[]

68.976

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

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

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))

In[]

quantile(data, [0, 0.25, 0.5, 0.75, 1,])

Out[]

[42.0, 63.0, 75.0, 82.0, 100.0]

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

In[]

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

In[]

percentile(data, [0, 25, 50, 75, 100,])

Out[]

[42.0, 63.0, 75.0, 82.0, 100.0]

検算(?)

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

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 クォンタイルの差です。難しいことはありませんね。よく見ると Quantile ではなく Quartile です。日本語では四分位範囲といい、クォーターとかそっち系列なのがわかります。おい強調構文のパーサおかしいぞ。聞いているのか Qiita。

In[]

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

In[]

iqr(data)

Out[]

19.0

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

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

In[]

def niqr(data):
    return iqr(data) / 1.3490

In[]

niqr(data)

Out[]

14.084507042253522

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

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

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

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

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

\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\}$ を考え、これを床函数の代わりに置くと、以下に示すように整数の添え字でも適用できるようになります。

\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$ となります。

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

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

あとがきポエム

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

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