国立感染症研究所ののSARS-CoV-2の変異株B.1.1.529系統（オミクロン株）の潜伏期間の推定：暫定報告では，度数分布から中央値を求めるために，ガンマ分布でフィットして，中央値2.9日（95%信頼区間2.6-3.2）を出している。以下では特別なモデルを仮定せず各ビンの中では一様分布（分布関数が折れ線）として度数分布データから中央値（一般に分位数）およびその信頼区間を求める。
上述の「暫定報告」の「データ1」はオミクロン株曝露-発症間隔（0，1，2，3，4，5日）の度数分布を表す次のような $N = 35$ のデータである。
hist = [0, 2, 8, 15, 9, 1]
このようなデータの中央値を求めるのが問題である。
$N = 35$ であるから，18番目の値が中央値である。18番目は3のビンに入るので，3のビンの15個がすべて3に等しければ，中央値は3である。
そうではなく，ここでは3のビンのデータは2.5から3.5までの範囲に一様分布すると仮定しよう。ここで考えているデータでは，0日目の正午に暴露したと仮定するのと同じことである。もし3のビンのデータが3から4までの範囲とするのであれば，以下の結果に0.5を加えればよい。
このデータでは，3未満のビンには10個あるのだから，3のビンの15個中で8番目が中央値である。15個の8番目はちょうど中央であるので，こう考えても中央値は 3.0 になるとしてよさそうである。国立感染症研究所の結果（2.9日）ともほぼ一致する。
これは，次の累積分布で 35/2 のところを求めることに相当する。
より一般的に，度数分布で与えられる分布で，各ビン内では一様分布と仮定する場合，中央値は次の式で計算できる：
この場合にあてはめると 2.5 + (35/2 - 10) / 15 * 1 = 3.0 となり，さきほどの結果と一致する。
より一般に，階級値 0, 1, 2, 3, … についての度数のリスト hist が与えられたとき，このような中央値を求めるには，次のようにする：
def hmedian(hist):
n = sum(hist)
i = 0
s = hist[i]
while s < n / 2:
i += 1
s += hist[i]
return i + (n / 2 - s) / hist[i] + 0.5
これでいいと思っていたら，[1, 1, 0, 1, 1] のような真ん中が 0 の場合に，対称性を破る値を出してしまう。
次のようにすればよさそうである：
def hmedian(hist):
n = sum(hist)
i = 0
s = hist[i]
while s < n / 2:
i += 1
s += hist[i]
if s == n / 2 and hist[i + 1] == 0:
j = i + 2
while hist[j] == 0:
j += 1
return (i + j) / 2
return i + (n / 2 - s) / hist[i] + 0.5
ここで n / 2 を別の値に変えれば任意の分位数が求められる。このことを使って一般化したのが次の関数である。q に0より大きく1より小さい値を与える。例えば0.25を与えれば四分位数を出力する。
def hquantile(hist, q=0.5):
n = sum(hist)
t = n * q
i = 0
s = hist[i]
while s < t:
i += 1
s += hist[i]
if s == t and hist[i + 1] == 0:
j = i + 2
while hist[j] == 0:
j += 1
return (i + j) / 2
return i + (t - s) / hist[i] + 0.5
データがヒストグラムではなく自然数のリストで与えられているときは，例えば次のようにしてヒストグラムに直す：
import numpy as np a = [0, 1, 1, 3, 5] # リストで与えられたデータ hist, bin_edges = np.histogram(a, range(max(a) + 2))
bin_edges は無視すればよい。同様に，plt.hist() や ax.hist() でヒストグラムを描く際にも，通常は無視する戻り値を使えばよい：
h, bins, _ = ax.hist(a, range(int(np.nanmin(a)), int(np.nanmax(a)) + 2), align="left")
print("Median =", hmedian(h) + bins[0])
逆に，ヒストグラムをデータに直すには次のようにする：
np.repeat(range(len(hist)), hist)
中央値については中澤港先生の作られたRのパッケージ fmsb の truemedian という関数にすでに実装されていることを中澤先生にご指摘いただいた。truemedian は次のような実装になっている：
truemedian <- function (X, h = 1)
{
YY <- rep(0, length(X))
XX <- table(X)
q <- length(XX)
k <- 0
for (i in 1:q) {
L <- as.numeric(names(XX)[i]) - h/2
for (j in 1:XX[[i]]) {
k <- k + 1
YY[k] <- L + h * (2 * j - 1)/(2 * XX[[i]])
}
}
return(median(YY))
}
これを以下と同じ仕様にして Python で書けばこんな感じ：
def hmedian(hist):
a = np.array([])
for i, h in enumerate(hist):
if h != 0:
a = np.append(a, np.arange(1, 2 * h, 2) / (2 * h) + (i - 0.5))
return np.median(a)
ブートストラップを使えば分位数の信頼区間が簡単に求められる。
rng = np.random.default_rng(20220129)
hist = [0, 2, 8, 15, 9, 1]
m = len(hist)
n = sum(hist)
x = np.repeat(range(m), hist)
def func():
y = rng.choice(x, n)
h, _ = np.histogram(y, range(m + 1))
return hmedian(h)
r = [func() for _ in range(10000)]
np.quantile(r, [0.025, 0.975])
array([2.66666667, 3.33333333])
中央値の95%信頼区間はだいたい [2.67, 3.33] である。
この方法では区間の上下限は ..., 2.650, 2.656, 2.667, 2.675, 2.679, ..., 3.325, 3.333, 3.344, ... のような離散的な値しかとりえない。また，このようなブートストラップでは信頼区間が狭めに出ることが知られている（ファクタ $\sqrt{(n - 1) / n}$ くらい）。
遠藤先生がもっと簡単かつ正確な方法を考えてくださった：
rng = np.random.default_rng(20220201)
hist = [0, 2, 8, 15, 9, 1]
n = sum(hist)
x = np.repeat(range(len(hist)), hist)
def func():
return np.median(rng.choice(x, n) + rng.random(n)) - 0.5
r = [func() for _ in range(10000)]
np.quantile(r, [0.025, 0.975])
array([2.62720809, 3.37747261])
r = [func() for _ in range(100000)] np.quantile(r, [0.025, 0.975])
array([2.62591919, 3.37426923])
中央値の95%信頼区間はだいたい [2.63, 3.37] である。さきほどの結果よりほんの少し広くなった。
この方法なら信頼区間の上下限は稠密な値をとりうる。
ある点 $L$ が中央値の95%信頼区間の下端であれば，$L$ 以下に半数以上のデータが入る確率は 0.025 である。分布がヒストグラム（各ビン内では一様分布）で与えられていれば，$L$ 以下の面積の割合 $p$ がヒストグラムから求められる。2項分布を仮定すれば
\[ \sum_{r=\lceil n/2\rceil}^n {}_nC_r p^r (1-p)^{n-r} = 1 - \sum_{r=0}^{\lceil n/2\rceil - 1} {}_nC_r p^r (1-p)^{n-r} = 0.025 \]これを満たす $p$ を求めればよい。ここで第2辺の和は scipy.stats.binom.cdf()（R の pbinom()）で求められる。ただ，$n$ が偶数のときはこれでは正確な中央値が求められないので，とりあえず $n$ を 1 増やして奇数にしている（もっといい手はないか）。
import scipy.stats
import scipy.optimize
hist = [0, 2, 8, 15, 9, 1]
n = sum(hist)
def foo(p):
if n % 2 == 0:
n += 1
return 1 - scipy.stats.binom.cdf((n - 1) // 2, n, p) - 0.025
p = scipy.optimize.root_scalar(foo, bracket=(0, 0.5))
hquantile(hist, q=p.root), hquantile(hist, q=1-p.root)
(2.6264133025787975, 3.3735866974212025)
ブートストラップ（遠藤先生バージョン）とほぼ同じ結果になった。
こういうやりかたは既知なんだろうか。Rの DescTools パッケージの MedianCI という関数は，本質的な部分は次のようになっている：
k <- qbinom(p = (1 - conf.level)/2, size = n, prob = 0.5, lower.tail = TRUE) ci <- sort(x)[c(k, n - k + 1)] attr(ci, "conf.level") <- 1 - 2 * pbinom(k - 1, size = n, prob = 0.5)
念のため，簡単な場合にブートストラップと上述の2項分布の方法による $p$ の比較：
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats
import scipy.optimize
rng = np.random.default_rng(20220207)
def func(n):
r = [np.median(rng.random(n)) for _ in range(100000)]
return np.quantile(r, 0.025)
a1 = [func(i) for i in range(1, 41)]
def func(n):
if n % 2 == 0:
n += 1
def foo(p):
return 1 - scipy.stats.binom.cdf((n - 1) // 2, n, p) - 0.025
p = scipy.optimize.root_scalar(foo, bracket=(0, 0.5))
return p.root
a2 = [func(i) for i in range(1, 41)]
plt.plot(range(1, 41), a1, range(1, 41), a2, "o-")
plt.legend(["bootstrap", "binomial"])
ヒストグラムから中央値を求める際に、中央値の入る階級内では一様分布と仮定してきた。しかし、その両隣の階級の度数から、一様分布からの外れが推測できるかもしれない。具体的に、次のグラフのように両隣まで各階級の面積が等しい放物線で中央の階級内の密度を推定してはどうだろうか。下のスライダーを動かして実験されたい。
count[i-1] |
||
count[i] |
||
count[i+1] |
Last modified: