« 消費者パネルデータの栄枯盛衰 | メイン | ちょっとした覚書：Stanのデータ型について »
2019年10月23日 (水)
たとえば、いまサイズ4の標本があり、ある変数の値が小さい順に100, 200, 300, 300であったとする。標本中央値は250ですわね。では40パーセンタイルは？
これは意外にややこしい話で、以前気が付いたんだけど、ソフトによって出力が異なるのである。
こんな細かいことに気を揉んでも一銭の得にもならん、くわばらくわばら...と思っていたが、このたび勤務先のお仕事の関係で関連するお問い合わせを受け、まあこれも機会であろうと思い、ちょっと調べてみた。
Hyndman, R.J. & Fan, Y. (1996) Sample Quantiles in Statistical Packages. American Statistician, 50, 361-365.
この話になるとよく引用される論文はこれだという印象だけがあったのだが、実際に手に入れてみたら、4pの短いコラムのような記事であった。
いわく。
分布関数を$F(x)$として、分位数とは
$Q(p) = F^{-1}(p) = inf \{x: F(x) \geq p \}, \ \ 0 < p < 1$
だが、これを標本においてどう定義するか。
以下、$i$番目の定義を$\hat{Q}_i(p)$とする。標本を$\{X_1, \ldots, X_n\}$とし、その順序統計量を$\{X_{(1)}, \ldots, X_{(n)}\}$とする[小さい順に並べて$3$番目の値を$X_{(3)}$と書きますってことでしょうね]。
$u$を超えない最大の整数を$L(u)$, $u$を下回らない最小の整数を$U(u)$と書く[←原文ではあまり見たことない記号を使っているので、このメモではこのように略記する。以下、$u$に負値が来ることはなさそうだから、要は$u$の端数を切り捨てたら$L(u)$, 切り上げたら$U(u)$ね]
標本分位数の定義はソフトによってさまざまである。一般化して書くと、
$\hat{Q}_i(p) = (1-\gamma) X_{(j)} + \gamma X_{(j+1)}$
ただし $\frac{j-m}{n} \leq p \leq \frac{j-m+1}{n}$、$m$は実数、$0 \leq \gamma \leq 1$。
$\gamma$は $j = L(pn + m)$と$g = pn + m - j$の関数である。
[著者らは頭が良すぎるせいか抽象化しすぎているので、具体例を使って、私にもわかるように書き換えよう。
標本$\{100, 200, 300, 300\}$(サイズ$n=4$と書く)の40パーセンタイル($p=0.4$と書く)とは、$j$番目の値と$j+1$番目の値の間の区間のどこかにある点である。その点からその区間の左端までの長さが区間の幅の$(1-\gamma)$倍であり、その点からその区間の右端までの長さが区間の幅の$\gamma$倍であるとしよう。ソフトによるちがいは、結局は$j$と$\gamma$の決め方の問題にすぎない。
$j$の決め方のひとつは、$np = 4 \times 0.4 = 1.6$ の端数を切り捨てた値とする、というものである。このデータ例だと$j=1$、つまり100と200の間ということになる。ソフトによっては端数を切り捨てる前に$m$を足すものもある。たとえば$m=0.5$なら$j=2$、つまり200と300の間ということになる。
$\gamma$の決め方にはいろいろあるんだけど、どの決め方にせよ、それは$j$と、切り捨てた端数$g$のふたつによって決まる。]
さて、標本分位数$\hat{Q}_i(p)$に期待される特性として次の6つを挙げることができる。
- P1: 連続量であること。
- P2: $Freq(X_k \leq \hat{Q}_i(p)) \geq L(pn)$であること。[Hyndmanさんのブログ記事の誤植リストに従って修正した。上の具体例でいうと、40パーセンタイルの左側にある値の個数が、$0.4 \times 4 = 1.6$の端数を切り捨てた値、すなわち1以上であってほしい。そりゃそうだ]
- P3: $Freq(X_k \leq \hat{Q}_i(p)) = Freq(X_k \geq \hat{Q}_i(1-p))$ [上の具体例でいうと、40パーセンタイルより左の値の個数は, 60パーセンタイルより右の値の個数と同じであってほしい。うん、そりゃそうあってほしいわね]
- P4: $\hat{Q}^{-1}_i (X)$が一意に定義されるとき、$k=1,\ldots,n$について $\hat{Q}^{-1}_i (X_{(k)}) + \hat{Q}^{-1}_i (X_{(n-k+1)}) = 1$ [えーと、たとえば左からk番目の値が20パーセンタイルだったら、右からk番目の値は80パーセンタイルであってほしいということだ。P3と同じく、データの扱いが左右対称であってほしいということであろう]
- P5: $\hat{Q}^{-1}_i (X)$が一意に定義されるとき、$\hat{Q}^{-1}_i (X_{(1)}) \geq 0$かつ$\hat{Q}^{-1}_i (X_{(n)}) \leq 0$ [最小値が0パーセンタイルというのは変だし最大値が100パーセンタイルというのも変だということであろう]
- P6: $\hat{Q}_i(0.5)$が通常定義される標本中央値(奇数個だったら真ん中の値、偶数個だったら真ん中にある二つの値の中点)であること。
大変ながらくお待たせしました。それでは！選手の入場です！
定義1.
$m=0$とする。$g > 0$のとき$\gamma = 1$, $g = 0$ のとき$\gamma = 0$とする。
この定義は以下の性質を持つ。
$Freq(X_k \leq \hat{Q}_1(p)) = U(pn)$
$Freq(X_k \geq \hat{Q}_1(1-p)) = L(pn+1)$
P2を満たす。また$n$が奇数の時P6を満たす。他の特性は満たさない。
[上の具体例でいうと、$j$は$0.4 \times 4 = 1.6$の端数を切り捨てた値, すなわち$1$。切り捨てが起きているので$\gamma = 1$、よって答えは$j+1$番目の値200となる。それ以下の値の個数は2、これは$1.6$の端数を切り上げた値に等しい。なお、60パーセンタイルは300、それ以上の値の個数は2となり、これは$1.6+1$の端数を切り捨てた値に等しい]
定義2.
$m=0$とする。$g > 0$のとき$\gamma = 1$, $g = 0$ のとき$\gamma = 0.5$とする。
この定義は以下の性質を持つ。
$Freq(X_k \leq \hat{Q}_2(p)) = Freq(X_k \geq \hat{Q}_2(1-p)) = U(pn)$
P2, P3のみを満たす。
[定義1とほぼ同じなのだが、$pn$が端数を持たなかったときの扱いが異なる。たとえば上の具体例での25パーセンタイルは、定義1では100となるが、定義2では100と200の中点すなわち150となる。なんだか奇妙な定義だが、たとえば20パーセンタイル以下の値の個数と80パーセンタイル以上の値の個数が常に一致し$U(0.2n)$となる、つまりP3の特性を持つというメリットがあるわけだ]
定義3.
$m=-0.5$とする。$g > 0$のとき$\gamma=1$とする。$g=0$の場合、jが偶数ならやはり$\gamma=0$, 奇数の場合は$\gamma=1$とする。
この定義はP2のみを満たす。
[上の具体例でいうと、40パーセンタイルを求めるとして、$j$は$np - m = 1.6-0.5= 1.1$の端数を切り捨てた値1となる。切り捨てが起きたので2番目の値200を採用する。切り捨てが起きなかったときがなんだかややこしいが、まあその値だかその次の値だかを使うわけだ]
以上の定義では、分位数は$j$番目の値か$j+1$番目の値かその中点であり、つまり非連続的である(P1を満たさない)。
ここからは区間ごとに線形関数を考えて内挿するタイプの定義で、すべてP1を満たす。
以下、$\hat{Q}_i (p_k) = X_{(k)}$と書く。
[つまりですね、横軸に$X_{(k)}$、縦軸に$\hat{Q}_i (p_k) = X_{(k)}$となる$p_k$をとり、点$k=1, \ldots, n$が布置する散布図を書いて、点と点を直線でつなぐ。で、たとえば水平線$Y=0.4$との交点を求めてこれを40パーセンタイルとしようという算段である。問題は$p_k$をどう定義するかだ]
定義4.
$p_k = k/n$として線形に内挿する。
この定義はP1, P2を満たす。
[上の具体例でいうと、40パーセンタイルを求めるなら、点$(100,1/4)$と点$(200,1/4)$を結ぶ線分と、水平戦$Y=0.4$との交点をもとめるわけね。40パーセンタイルは160となる。素直な発想だが、このやり方だと50パーセンタイルは200になり、通常の標本中央値と合致しなくなる、つまりP6は満たされない]
定義5.
$p_k = (k-0.5)/n$として線形に内挿する。
P1からP6のすべてを満たす。
[上の具体例でいうと、40パーセンタイルを求めるなら、点$(200,1.5/4)$と点$(300,2.5/4)$を結ぶ線分と、水平線$Y=0.4$との交点をもとめるわけだ。40パーセンタイルは 210。50パーセンタイルは250となり、通常の標本中央値と合致する]
定義4-5には分布の分位数を推定するという視点が欠けている。いっぽう定義6-8は、$F(X_{(k)})$のなんらかの特性値を$p_k$にするというもの。いま$F(X_k)$が一様分布とすれば、$F(X_{(k)})$はベータ分布$Beta(k, n-k+1)$になるので、そのなんらかの特性値を使う。[←この理屈がいまいちぴんとこないんだけど... $F(X_k)$は実際には一様分布じゃないっすよね...]
定義6.
$p_k = k / (n+1)$として線形に内挿する。
これは$F(X_{(k)})$の期待値である。P3以外の5つの特性を満たす。
[上の具体例でいうと、40パーセンタイルを求めるなら、点$(200,2/5)$と点$(300,3/5)$を結ぶ線分と$Y=0.4$との交点を求めるわけだ。40パーセンタイルは200, 45パーセンタイルだったら225]
定義7.
$p_k = (k-1) / (n-1)$として線形に内挿する。
これは$F(X_{(k)})$の最頻値である。P1,P2,P4,P6を満たす。
[上の具体例でいうと、40パーセンタイルを求めるなら、点$(200,1/3)$と点$(300,2/3)$を結ぶ線分と$Y=0.4$との交点を求めるわけだ。答えは220]
定義8.
$p_k = (k-1/3)/(n+1/3)$として線形に内挿する。
これは$F(X_{(k)})$の中央値の近似である。P3以外の5つの特性を満たす。
[上の具体例でいうと、40パーセンタイルを求めるなら、点$(200,1.666/4.333)$と点$( 300,2.666/4.333)$を結ぶ線分と、$Y=0.4$との交点を求めるわけってことか...答えは206.66]
定義9は視点が違って、$F$が正規分布だと仮定し、$F(E(X_{(k)}))$を求めるというもの。定義6では$F(X_{(k)})$の期待値を考えたが、ここでは$X_{(k)}$の期待値を考える。次の定義で近似できる由。
$p_k = (k-3/8)/(n+1/4)$として線形に内挿する。
P3以外の5つの特性を満たす。
[上の具体例でいうと...めんどくさいので省略するが、40パーセンタイルは207.5となる模様]
特性だけでいえば優れているのは定義5である。いっぽう推定という観点からは、定義8は$F$の形状を問わず$Q(p)$のmedian-unbiasedな推定値であり、これがもっとも優れている。
思えば、標本分散の定義もいろいろありうるが(直観的にはnで割るのが分かりやすいし、正規分布ならn+1で割ればMSE最小になる)、不偏性を重んじてn-1で割るというのが標準になっている。いっぽう標本分位数の定義は標準化されていない。由々しいことだ。ソフトメーカーのみなさん、定義8に統一しませんか。
... 細かいところはよくわかんないんですけど、勉強になりましたです。
ところで、第1著者Hyndman先生は2016年のブログ記事でこの論文について振り返っている。いわく、この論文はある時期からすごく引用されるようになった(そのきっかけのひとつはWikipediaに論文へのリンクが載ったからだが、実は先生がご自分で載せたのだそうだ。ここ笑っちゃった)。しかし思いの丈は果たされなかった。各ソフトウェアは標本分位数のデフォルトを定義8に統一するのではなく、むしろこれまで実装していなかった定義も実装するようになってしまった... とぼやいておられる。面白いなあ。
現時点での各種ソフトはどうなっているだろうか。
Rの場合。そもそもRのstats::quantile()を書いた人のひとりがHyndman先生その人であるからして、この論文の定義の番号に合わせ、たとえば quantile(x, type = 8) という風に書けばよい。ただしデフォルトは先生の意に反し、S言語での定義に合わせた定義7となっている。
私は使っていないけど、時々お問い合わせを受けるのがSPSSの挙動である。どうなっているのか。
手元にあるIBM SPSS 24のStatistic AlgorithmsというPDFによると(ここからダウンロードできる)、単変量の記述統計量を求めるEXAMINEコマンドは、分位数として5つの定義を実装している。FAQページにいわせるとデフォルトはHAVERAGEという定義だ。
PDFによればHAVERAGEの定義は次の通り。
標本サイズを$W$, $p$パーセンタイルを$x$, 重複を取り除いた$i$番目の値を$y_i$, その度数を$c_i$, 累積度数を$cc_i$と表記する。$tc_2 = (W+1)p$とし、$cc_i$で区切られる区間のうち$tc_2$が落ちる区間の下限を$cc_{k_2}$とする。$g^{*}_2 = tc_2 - cc_{k_2}$, $g_2 = g^{*}_2 / c_{k_2+1}$として、
- $g^{*}_2 \geq 1$のとき $x = y_{k_2+1}$
- $g^{*}_2 < 1$ かつ $c_{k_2+1} \geq 1$のとき$x = (1-g^{*}_2)y_{k_2} + g^{*}_2 y_{k_2+1}$
- $g^{*}_2 < 1$ かつ $c_{k_2+1} < 1$のとき$x = (1-g_2)y_{k_2} + g_2 y_{k_2+1}$
うむむ、ややこしい。
話がややこしくなっているのは、SPSSではローデータの個々の行に対して整数でない頻度ウェイトを与えることができるからである。
話を簡単にするために、頻度ウェイトは1以上の整数としよう(頻度ウェイトとは標本におけるその行の頻度であるから、本来は1以上の整数になるはずである)。このとき$y_i$の度数$c_i$は1以上の整数になるので、3本目の式は無視できる。
1本目、$tc_2$が区間下限$cc_{k_2}$より1以上離れていたら区間の上限を採用するというのは、たとえば$\{100, 200, 300, 300\}$の65パーセンタイルを求めるとき、$tc_2 = 0.65 * (4+1) = 3.25$が落ちる区間は累積度数2と累積度数4の間だけど、下から3人目以降はあきらかに300なんだから、$tc_2$が3以上だったら四の五の言わずに300だ、ということだと思う。SPSSでいう$y_i$はHyndmanさんのいう$X_{(i)}$とちがって、データの重複を取り除いたうえで小さい順に並べたときの$i$番目の値だから、こういう但し書きが必要になるわけだ。
話のポイントは2本目の式である。最初に挙げた具体例$\{100, 200, 300, 300\}$で考えると、その40パーセンタイルは、$tc_2 = 0.4 * (4+1) = 2$, これが落ちる区間の下限は2だから$g^{*}_2 = 0$で、2番目の式より200となる。45パーセンタイルなら$tc_2 = 0.45 * (4+1) = 2.25$, $g^{*}_2 = 0.25$となり、0.75*200 + 0.25 * 300 = 225。どちらの例でもHyndmanさんのいう定義6に一致する。
それもそのはずで、Hyndmanさんの説明は「点$(200,2/5)$と点$(300,3/5)$を結ぶ線分と水平線$Y=0.45$との交点を求める」であり、SPSSの説明は「点$(200,2)$と点$(300,3)$を結ぶ線分と水平線$Y=0.45 \times 5=2.25$との交点を求める」である。データ点の縦座標を5で割るか、欲しい$p$に5を掛けるかの違いで、結局同じ事をいっている。
このように、頻度ウェイトが正の整数であれば、SPSSの出力は定義6だと思う。
Excelはどうなっているんだろう... 上の具体例で求めてみたところ、セル関数PERCENTILE.INCは40パーセンタイルとして220を返し、PERCENTILE.EXCは200を返した。調べる気力が尽きたが、これ、なにやってんですかねえ。
論文：データ解析(2018-) - 読了：Hyndman & Fan (1996) 標本分位数の出力はソフトによってちがう