これはLAPRAS Advent Calendar の7日目です。

はじめに

LAPRASでは、ITエンジニアのアウトプット活動などからその人のスキルを数値化した「技術力スコア」を提供しています。
今回はこのスコアに関連して、数学の小ネタをお話ししようと思います。

LAPRASのスコアは偏差値のように相対的に計算しており、平均3.0、分散0.5となるようにスケーリングしています。
（偏差値の場合は平均50、分散10となるようにスケーリング）

スコアが正規分布に従う場合、スコアの値と「上位何%」という数字は一対一に対応するはずです。
例えば、偏差値80（LAPRASスコア4.5）は上位0.13%、偏差値70（LAPRASスコア4.0）は上位2.28%、偏差値60（LAPRASスコア3.5）は上位15.87%、そして偏差値50（LAPRASスコア3.0）はちょうど上位50%になるはずです。

画像の出典: Wikipedia

しかしながら、現在のLAPRASスコアは厳密にこの値となっているわけではなく、かなりのズレが生じてしまっています。
逆に、これらの値をピッタリと合わせるようにスコアを計算する方法はあるでしょうか？

言い換えると、「全体の上位何％」という順位の情報から、正規分布に従うスコアを計算するにはどうしたらいいでしょうか？
今回の記事では、このことについて少し深ぼって考えたいと思います。

ナイーブな考え

上記の考えをそのままナイーブに適用することを考える。

ある人のスコアが、N人の中で下からn番目であるとき、その人のパーセンタイル値（下から何%に位置するか）は、
$$
\frac{n-1}{N-1} \times 100 [\%]
$$
で計算される。
これは、ビリの時には0%、トップの時には100%となる値である。

しったがって、正規分布の累積分布関数を $\Phi(x)$ とすると、正規分布に従うスコアが $x$ のとき、$N$が十分大きいならば

$$
\Phi(x) = \frac{n-1}{N-1}
$$

が成り立つはずである。そこで、$\Phi(x)$の逆関数を用いることで、順位からスコアを以下のように求めることができる。

$$
x = \Phi^{-1}\left( \frac{n-1}{N-1} \right)
$$

これは非常にナイーブな考え方で一見合理的だが、少々問題がある。
この式は$N\to\infty$の場合には成り立つが、Nが有限の場合には必ずしも成り立たない。特に、この式はトップの場合（n=N）には $x\to+\infty$ に、ビリの場合（n=1）には $x \to -\infty$ に発散してしまうため、これらの場合にスコアを考えることができない。
また、同じ1位だったとしても、10人の中での1位と100万人の中での1位では後者の方がよりスコアが高くなるようになるべきである。
このように、最高位や最下位も正しく扱えるような枠組みを考える必要がある。ここでは、確率モデルによる定式化を行う。

確率モデル

以下では、平均0、分散0.5の標準正規分布を考える。また、標準正規分布に従う能力値の指標を指して単に「スコア」とよぶ。
（偏差値やLAPRASスコアのスケールに変換するには、定数倍して並行移動すれば良い）

N人のユーザーのスコアが同じ標準正規分布に従っているとする。それぞれの値はわからないが、順位だけが観測されたものとする。このとき、下からn番目のユーザーのスコアの値はどれくらいになるか？
これを確率モデルとして考える。

言い換えると、標準正規分布に従うN個確率変数を考え、それらを小さい順に $x_1, \ldots, x_N$ としたとき、$n$番目の値 $x_n$ の分布を考える。

N人のうち1人に着目したとき、そのスコアの値が $x $ で、順位が下から$n$ 位としたとき、求める分布は、条件付き分布 $p(x|n)$ である。

これはベイズの定理から以下のように計算される:

$$
p(x|n) = \frac{p(n|x) \, p(x)}{p(n)}
$$

ここで、

$p(x) = \mathcal N(x)$
- 順位を指定しない場合、スコアは標準正規分布に従う
$p(n) = \frac{1}{N}$
- 対称性から、スコアが不明ならばどの順位になる確率も等しい
$p(n|x) = \binom{N-1}{n-1} \Phi(x)^{n-1} (1-\Phi(x))^{N-n} $
- 自分のスコアが $x$ のとき、自分が下からn番目になる確率
  = 標準正規分布からN-1回サンプリングした時、$x$ より低い値がn-1回得られる確率
  = 確率 $\Phi(x)$ の事象をN-1回観測した際の二項分布

ただし、$\mathcal N (x)$ は標準正規分布の密度関数、$\Phi(x)$ は標準正規分布の累積分布関数である。

以上から、求める分布は、

$$
\begin{aligned}
p(x|n)
&= \frac{(l+m+1) !}{ l! \, m!} \cdot \mathcal {N}(x) \cdot \Phi(x)^l \cdot \Phi(-x)^m && \cdots (1)
\end{aligned}
$$

と求まる。
ただし、

$l = n-1$ は「自分よりも順位が下の人の人数」
$m = N-n$ は「自分よりも順位が上の人の人数」

である。また、$\Phi(-x) = 1-\Phi(x)$ という関係を用いた。

(1)式を $x$ について全域で積分すると、

$$
\begin{aligned}
\int_{-\infty}^{\infty} p(x|n) \,dx
&= \frac{(l+m+1) !}{ l! \, m!} \int_{-\infty}^{\infty} \mathcal {N}(x) \cdot \Phi(x)^l \cdot \Phi(-x)^m \, dx \\
&= \frac{(l+m+1) !}{ l! \, m!}\int_{0}^{1} z^l (1-z)^m \, dz
& (\because z = \Phi(x)) \\
&= \frac{(l+m+1) !}{ l! \, m!} \, B(l+1, m+1) & (Bはベータ関数) \\
&= 1
\end{aligned}
$$

となり、確かに確率分布になっていることが確認できる¹。

最頻値を求める

スコアの推定値として、確立分布(1)の最頻値(最大値)を用いることを考える。
これは $\frac{d}{dx} \left( \log p(x|N) \right) = 0$ となる $x$ を計算することで求まる。

$$
\begin{align}
\frac{d}{dx} \left( \log p(x|N) \right) = -x + l \cdot \frac{\mathcal{N(x)}}{\Phi(x)} - m \cdot \frac{\mathcal{N(x)}}{\Phi(-x)} = 0
\;\;\;\; \cdots (2)
\end{align}
$$

$h(x) := \frac{x \Phi(x) \Phi(-x)}{\mathcal{N(x)}}$ としてこの式を整理すると、

$$
\Phi(x) = \frac{l - h(x)}{m+l} \;\;\;\; \cdots (3)
$$

となる。

ここで、関数 $h(x)$ は下図のようなシグモイド型の関数であり、とくに$-1 \leq h(x) \leq 1$ である。

$l, m$ の値が大きい（つまり、自分の上にも下にも十分な人数がいる）場合、$|h(x)| < 1 \ll l, m$ だから、(3) 式において $h(x)$ は無視してしまってよい。

このとき、最頻値は

$$
x = \Phi^{-1}\left(\frac{l}{l+m}\right) = \Phi^{-1}\left(\frac{n-1}{N-1}\right) \;\;\;\; \cdots (4)
$$

と求まる。
これは、前述の「ナイーブな考え」の結果と一致している。
ただし、上で考察したように、この式は最下位($l=0$)の場合はにスコアが$-\infty$に、最上位($m=0$)の場合は$\infty$ に発散してしまうため、最下位や最上位付近では $h(x)$ の効果を無視することはできない。

Newton法による数値解

最下位や最上位付近でも正しいスコアを計算するため、(2)式をNewton 法によって数値的に求める。

(2) 式を少し変形して、
$$
f(x) := (m+l) \, \Phi(x) + h(x) - l =0
$$

この式の微分

$$
f'(x) = (m+l) \, \mathcal{N}(x) + (1+x^2) \frac{\Phi(x)\Phi(-x)}{\mathcal{N}(x)} - x(2\Phi(x) -1)
$$

を用いて、更新式

$$
x_{n+1} = x_{n} - \frac{f(x_n)}{f'(x_n)}
$$

が収束するまで更新する。
Python で実装は以下のようになる。

from scipy import stats

def _f_df(x, l, m):
    """ Newton 法の更新幅 f(x) / f'(x) の値を求める
    """
    pdfx = stats.norm.pdf(x)
    cdfpx = stats.norm.cdf(x)
    cdfmx = stats.norm.cdf(-x)

    f = (m + l) * cdfpx + x * cdfpx * cdfmx / pdfx - l
    df = (m + l) * pdfx + x * (1 - 2 * cdfpx) + (1 + x * x) * cdfpx * cdfmx / pdfx
    return f / df

def newton(l, m, x0=0.0, iterations=40):
    """ Newton 法により、順位から標準正規分布のスコアを計算する
    """
    x = x0
    for i in range(iterations):
        x -= _f_df(x, l, m)
    return x

このアルゴリズムは、$N \leq 1,000,000 $ 程度であれば、反復20回以内で収束した。

以下に、Newton法による数値解(newton)と上記の「ナイーブな考え」の式(4)による結果(naive)との比較を示す。

確かに、中央付近では両者の結果がよく一致しているが、上位と下位の両端ではその差は大きくなり、naive では最上位と最下位で値が発散してしまうところを、Newton法の数値解では有限の値に落ち着いている。

また、実際に正規分布からN個のスコアをサンプルした実験結果と比較してみた。

図は「N=50人中のトップの人のスコア」の分布（1万回サンプリングしたヒストグラム）であり、破線はNewton法による最頻値の推定値である。
この図を見ると、確かにNewton法による結果は最頻値を正しく推定できていることがわかる。

近似値の計算

計算量と数値的な安定性を高めるため、Newton法などの反復法によらない陽な形式での近似値を計算したい。

導出は省略するが、$h(x)$ の $x→\infty$ における振る舞いを$\Phi(x)$を用いて近似することで、
$$
h(x) \approx 1 - \frac{1}{2\log(1-\Phi(x))}
$$
という関係が成り立っている。
これを利用して、
$$
h(x) \approx \frac{(l-m)\left( 1-\frac{1}{2 \log(l+m+2)} \right)}{l + m + 2\left( 1-\frac{1}{2 \log(l+m+2)} \right)}
$$
という近似を用いるものとする。

このとき、$\varepsilon_N := 1-\frac{1}{2 \log(l+m+2)} = 1-\frac{1}{2 \log(N+1)} $ とすれば、最頻値の近似値を

$$
x = \Phi^{-1} \left( \frac{l-h(x)}{l+m} \right)\approx \Phi^{-1}\left(\frac{l+\varepsilon_N}{ l + m + 2 \varepsilon_N}\right)
$$

によって求めることができる。
これは、「ナイーブな考え」の式(4) において $l, m$の値にそれぞれ$\varepsilon_N$ の値を加えたものになっており、つまりは「仮想的に最上位の上と最下位の下にそれぞれ $\varepsilon_N$ 人 ($0 < \varepsilon_N < 1$) がいる」と考えることに相当する。

Newton 法の結果と比べてみると、この近似解(approx)はNewton法の数値解(newton)とほとんど重なっており、良い近似となっていることがわかる。
この誤差は$N$ が大きくなるにつれて小さくなり、 $N \geq 100$では 0.05以下となるため、実用上は十分な精度であるといえる。

現実世界への適用

せっかくなので、現実世界の問題に当てはめてみる。

よく「人の身長の分布は正規分布になる」と言われる。これに今回の考えを当てはめて「身長が正規分布に従っているならば、1位の人の身長はどの程度と予測されるか？」を考えてみる。

政府の「国民健康・栄養調査(2019年)」によると、20歳以上の日本人男性の平均身長は167.7、標準偏差は6.9とのことである²。また、現在の20歳以上の日本人男性の人口は約5,050万人とのことである³。

これらから、Newton法により1位の人の身長を推定すると、

In [1]: newton(50500000, 0) * 6.9 + 167.7
Out[1]: 205.63639314207018

となり、約206cmが推定値と求まる。これは偏差値に換算すると偏差値105に相当する。

存命中の日本人として身長が最も高いのは元バスケットボール選手の岡山恭崇さんの230cmと言われており⁴、これは偏差値140に相当する値であある。

なぜこのように予測と実際に差が生じたのだろうか？
これは元々の「身長は正規分布に従う」というのが厳密には正しくないためであると考えられる。
身長が正規分布に従うのは中心極限定理による漸近的なもので、これは「身長は独立な無数の因子の足し合わせによるもの」という仮定に立っているものである。
分布の中心付近ではこの仮定がある程度成り立っているとしても、たとえば「ごくまれに起こる身長が非常に大きくなる因子」が存在する場合、分布の端では正規分布への収束は遅くなる。

このように、「現実世界の観測値が厳密に正規分布に従っている」という仮定は分布の端付近ではなかなか成り立たないものであるということがわかった。
しかしながら、観測値そのものではなく、仮想的に「理想的な正規分布に従う隠れ変数」を考え、その隠れ変数をノンパラメットリックに推定しているものであると考えれば、現実世界観測値がどのような分布に従っていたとしても適用できる変数の正規化手法として有用なのではないかと考えられる。

まとめと宣伝

今回は、LAPRASのスコアや偏差値に関連して、「スコアがxxxのとき、上位xxx%です」という説明から、逆に、順位からスコアを計算する方法について考察しました。
ナイーブな手法では、最高位や最下位のスコアを決定することができませんでしたが、「正規分布に従うN人中n位の人のスコア」を確率モデルとして扱うことで、分布の中央付近だけでなく端付近でも適用可能な理論的枠組みと最頻値の計算方法を提案しました。

LAPRASでは、この考え方を取り入れた、スコアを安定的でわかりやすいものにする改修を行う予定ですので、お楽しみにしてください。

$z=\Phi(x)$ の変数変換により、この分布はベータ分布 $\beta(z; l+1, m+1)$ とることがわかる。 ↩
https://www.e-stat.go.jp/dbview?sid=0003224177 ↩
人口推計（2022年（令和4年）10月1日現在） https://www.stat.go.jp/data/jinsui/2022np/index.html 別表1 ↩
https://ja.wikipedia.org/wiki/%E5%B2%A1%E5%B1%B1%E6%81%AD%E5%B4%87 ↩