この記事について、$\frac{N+1}{2}$を使用するのは誤りではないかとのご指摘がありましたので、再度掲載いたします。基本的な考え方は前記事をお読みください。
岩原信九郎:『教育と心理のための推計学,日本文化科学社』,2001.(新訂版第32刷)P.41 にこのような記事があるそうです。
分布表から中央値を求めるときは、同点のある場合と同じ方法を用いる。分布表は級間(単位)を大きくしてその間の値をすべて同点にしたものだからである。しかしこのときは次の式を用いることが出来る。
Mdn = l + h\left(\dfrac{\frac{N}{2} - F}{fm} \right)
$l \cdots$中央値のある級間の下限点
$F \cdots l$以下の総度数
$fm \cdots$ 中央値のある級間の度数
$h \cdots$ 級間の幅
$N \cdots$ 総度数。奇数でも偶数でもよい
とのことです。そこで、比例配分を求める部分だけ修正を行ったコードを下記に示します。
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))
# 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 + ((N / 2 - CF) / f) * h),
false = (L + df_med$value[upper_row]) / 2)
print(med)
[1] 48.60554