比例配分による中央値を出すために

6月にこの記事を投稿したところ、「比例配分の中央値を出さないことのもったいなさに気付かないとは」というコメントをいただきました。なるほど、と思いながら、そのまま1ヶ月ほどが経ってしまいました。

せっかくやるのであれば、自分が中央値に対する理解を深めないと、いけませんよね。

quantile()関数で比例配分について考える

quantile.r

quantile(c(160, 164, 166, 169, 175, 182))

例えばこの6つのデータの四分位数を求めるとき、$Q1 = 164, Q2 = 167.5, Q3 = 175$と、人間はあまり考えることなく解答することができます。ところが、実際のRの結果は以下のとおり異なります。これはエクセルでも同様の結果になります。

   0%   25%   50%   75%  100% 
160.0 164.5 167.5 173.5 182.0

これを手で計算してみました。

データは6つありますが、1番小さいデータを0番目としてカウントすると、1番大きなデータは5番目($n-1 = 6-1 = 5$) となります。この値$5$を$4$で割った値$1.25$が第一四分位数のある位置と考えます。これは1番目の値164と2番目の値166の間にあり、しかも164と166の間の0.25、つまり$\frac{1}{4}$にあたる位置にあります。これを計算すると、$164 + (166 - 164) \times \frac{1}{4}=164.5$ということになります。これが比例配分のイメージと私は捉えました。

中央値は第一四分位数や第三四分位数よりイージーに考えられる

ただ、四分位数を考えるときのように最小値の位置を$0$番目の値と考えるのは文系の私にとっては本当につらい。ですから、結局最後は$1$を足して最小値を$1$番目の値としてカウントし直して考えたいです。データを昇順に並べたとき、$x_1, x_2, \cdots, x_n$の中央値の位置が、最小値を1番目、最大値を$n$番目としたときにどこなのかを考えると以下のように計算できます。

\frac{n-1}{2} + 1 = \frac{n-1}{2} + \frac{2}{2} = \frac{n+1}{2}

結果的に、データの数に1を足したものを2で割って算出された値が中央値の位置になるという、文系人間の直感にも当てはまる結果となるわけです。

中央値をプログラミングする

以下は上記のことを反映させてプログラミングしました。if_else()関数でいくつか分岐を作っているのは「比例配分にする必要のない場合は、比例配分計算は行わない」というコンセプトに基づいています。

propotional_median.r

library(conflicted)
library(tidyverse)

## 日本の年齢別人口(2020年国勢調査)データ
df_med <- dplyr::tibble(
  value = seq(from = 0, to = 100),
  n = c(831824, 866525, 910005, 934063, 973665, 998664, 996576, 1020657, 1024325,
        1048871, 1057736, 1063185, 1083113, 1077379, 1069104, 1070370, 1113159,
        1123237, 1151389, 1159285, 1177049, 1174456, 1193935, 1193983, 1191883,
        1210380, 1212167, 1192761, 1209983, 1206673, 1234225, 1258775, 1299569,
        1335672, 1356353, 1407719, 1456774, 1476522, 1478920, 1491632, 1558309,
        1596338, 1656969, 1699648, 1779813, 1852394, 1956602, 1991137, 1954205,
        1895955, 1837431, 1807894, 1764215, 1758955, 1371356, 1689290, 1582910,
        1542017, 1491947, 1461318, 1469672, 1494495, 1450282, 1407740, 1475001,
        1516441, 1512118, 1598839, 1679855, 1768015, 1886160, 2052521, 2014143,
        1893386, 1165585, 1234832, 1488292, 1416571, 1430462, 1360771, 1206058,
        1017635, 1048747, 1035057, 989231, 893754, 794162, 739149, 662252, 580506,
        493824, 425042, 358386, 277613, 224151, 169988, 123057, 90075, 64897,
        44707, 79523)) %>%
  dplyr::mutate(cumsum = cumsum(n))
head(df_med)
tail(df_med)

# A tibble: 6 × 3
  value      n  cumsum
  <int>  <dbl>   <dbl>
1     0 831824  831824
2     1 866525 1698349
3     2 910005 2608354
4     3 934063 3542417
5     4 973665 4516082
6     5 998664 5514746
# A tibble: 6 × 3
  value      n    cumsum
  <int>  <dbl>     <dbl>
1    95 169988 122812002
2    96 123057 122935059
3    97  90075 123025134
4    98  64897 123090031
5    99  44707 123134738
6   100  79523 123214261

propotional_median.r

# Nの数
N <- tail(df_med$cumsum, 1)
# 中央値に位置する累計値half_nを求める。
half_n <- (N + 1) / 2
# half_nが小数を含む場合(N+1が奇数の場合)の切り捨てと切り上げ
lower_n <- floor(half_n)
upper_n <- ceiling(half_n)
# half_nの位置（何行目に位置するか）の確認
lower_row <- min(which(lower_n <= df_med$cumsum))
upper_row <- min(which(upper_n <= df_med$cumsum))

# 中央値が含まれる階級の下限値Lを求める
# (中央値が含まれる階級が48歳の場合は48-49の階級と考えられるため、
# 下限値Lは48となる。)
lower_class <- slice(.data = df_med, lower_row)
L <- lower_class$value
# 階級の上限値(=次の階級の値）U(48歳の場合は48-49なのでU=49)を求める。
# もし、中央値が最大の階級にある場合は、比例配分できないため
# U <- L + 1として、h = 1になるようにする。
U <- if_else(lower_row == nrow(df_med),
             true = L + 1,
             false = slice(.data = df_med, lower_row + 1)$value)
# 階級幅hを求める。今回の場合は1歳刻みなのでU=1。
h <- U - L
# Lの階級の度数fを求める。
f <- lower_class$n
# その階級の直前までの累計CFを求める。中央値が最初の階級に含まれる場合は、
# CF=0にする。
CF <- if_else(lower_row == 1,
              true = 0,
              false = slice(.data = df_med, lower_row - 1)$cumsum)
# 中央値medの計算
# lower_row == upper_rowの時は比例配分で中央値を求める。
# lower_row < upper_rowの時は([lower_row] + [upper_row]) / 2で求める。
med <- if_else(lower_row == upper_row,
               true = if_else(f == 1,
                              true = L,
                              false =  L + ((half_n - CF) / f) * h),
               false = (L + df_med$value[upper_row]) / 2)
print(med)

[1] 48.60554