コペルニクス原理と余命の確率分布

コペルニクス原理に基づく余命の予測

J. Richard Gott III によって論じられたように、コペルニクス原理に基づくと、ある事象が発生してから時間 T が経過している場合、その事象は今後最小で T/39 、最大で 39T 続くということが 95% の確率で予測できます。

これは「事象の寿命を予測したい」と思った瞬間が事象の始まったばかりの瞬間、あるいは逆にその事象が終わる瞬間である可能性は殆どなく、事象の始まりから終わりまでのどこかに均等な確率で位置しているに違いない、という理屈に基づいています。

時間 $T$ だけ続いている事象の余命は確率 $p$ で

$t_{min}(p) = T / (\frac{2}{1-p} - 1)$

から

$t_{max}(p) = T (\frac{2}{1-p} - 1)$

までの範囲に収まる、と一般化することができます。

このような理屈が通用するのは余命を推測する方法が他に一切見当たらない場合だけで、もっと良い推測方法がある場合には不適切であることに注意が必要です（そうでないと「老人は子供より長生きする」ということになってしまう）。

ここでは、SageMath の練習を兼ねて思いつきを計算していきたいと思います。手で計算した内容を SageMath の Jupyter notebook 環境に記述する、という形で検算にもなるし Markdown でメモも記録できる、という両得を狙っています。

予測余命とその正答率の関係

予測余命と、その予測があたっている確率の関係をグラフ化することを考えてみます。

$t_{min}$ の定義域は $0 < t_{min} < T$ で、$ t_{min} = T / (\frac{2}{1-p} - 1) $ より $p(t) = \frac{T - t_{min}}{T + t_{min}}$

SageMath で $t_{min}(p)$ の逆関数 $p(t_{min})$ 、 $t_{max}(p)$ の逆関数 $p(t_{max})$ を出してみると

T = var('T')
assume(T > 0)

t_min = var('t_min', latex_name='t_{min}')
assume(t_min > 0)
assume(t_min < T)

t_max = var('t_max', latex_name='t_{max}')
assume(t_max > T)

p = var('p')

solve(t_min == T / (2 / (1-p) - 1), p)

$ \left[p = \frac{T - {t_{min}}}{T + {t_{min}}}\right] $

solve(t_max == T * (2 / (1-p) - 1), p)

$ \left[p = -\frac{T - {t_{max}}}{T + {t_{max}}}\right] $

$T=1$ と置いたときの $t_{min}$, $t_{max}$ それぞれの値を 横軸を時間、縦軸を $p$ としたグラフにプロットしたのが下です。与えられた確率 $p$ に対応する横線を書くと余命の予測範囲をグラフとの交点で囲まれた区間で示すことができて、例えば

• $p = 1$ の場合は「必ず当たる予測をしろ、と言われたら、今この瞬間余命は尽きるかもしれないし、もしかしたら永遠に続くかもしれない、くらいしかわからない」という答え
• $p = 0$ の場合は「まあ、今まで $T$ だけ続いたから今後も $T$ だけ続くんでは？全く当たる保証はないけど」という答え

を得ることができます：

plot((1-t_min)/(1+t_min), t_min, 0, 1) + plot(-(1-t_max)/(1+t_max), t_max, 1, 10)

「今までと同じくらいの時間続くんじゃ？」となにげなく言ってみたところ「そんないい加減なことを言うな！もっとちゃんと考えろ！」と迫られてしまった、そんな人が「じゃあ、もうすこし確度を上げて、このくらいの範囲で...」と答えを変えていくやりかたに、科学的な裏付けを与えられるのが面白いところです。

余命の確率分布

ここからさらにこんなことを考えてみます。

「余命を司る確率密度関数 $\mathscr{P}(t)$ が存在して

$ \int_{t_{min}}^{t_{max}} \mathscr{P}(t) dt = p$

と書けるのだという仮定を考えてみます。ここで $t_{min}$, $t_{max}$ はそれぞれ正答確率が $p$ となるように寿命を予測したときの最小時間と最大時間で

$ t_{min} = T / (\frac{2}{1-p} - 1) $

$ t_{max} = T (\frac{2}{1-p} - 1) $

$t_{min}$, $t_{max}$ はあくまでも予測結果にすぎず、実際の余命は $t_{min}$ より短いかもしれないし、 $t_{max}$ より長いかもしれない。しかしいつか必ず寿命が尽きる、という事実から $\mathscr{P}(t)$ は

$ \int_{0}^{\infty} \mathscr{P}(t) dt = 1$

となるはず。この $\mathscr{P}(t)$ はどのような形をしているだろうか。」

いわば、「事象の生気（プネウマ）」的なものが存在して、それが時間とともに減衰していった結果、事象が継続する活力がなくなっていくとともに余命が尽きる、というイメージです。

余命について一切の仮定をおくことができない場合の「生気」はどのように推移するのかを計算で出そう、というオカルト的発想です。

冗談はさておき、 $\mathscr{P}(t)$ を計算してみます。

$ \int_{t_{min}}^{t_{max}} \mathscr{P}(t) dt = p$

を $t=T$ を基準として

$0<t<T$ の領域

$ \int_{t_{min}}^{T} \mathscr{P}(t) dt = \frac{1}{2} p$

$T<t<\infty$ の領域

$ \int_{T}^{t_{max}} \mathscr{P}(t) dt = \frac{1}{2} p$

の２つに分割してそれぞれ計算してみましょう。

まずは $0<t<T$ の領域について。

$ t_{min} = T / (\frac{2}{1-p} - 1) $ より $p(t) = \frac{T - t}{T + t}$ なので

$ \int_{t}^{T} \mathscr{P}(\tau) d\tau = \frac{1}{2}\frac{T - t}{T + t}$

$\tau=0$ からの定積分で書き直すと

$ \int_{0}^{T} \mathscr{P}(\tau) d\tau - \int_{0}^{t} \mathscr{P}(\tau) d\tau = \frac{1}{2}\frac{T-t}{T+t}$

両辺を $t$ で微分すると

$ \frac{d}{d t} \int_{0}^{T} \mathscr{P}(\tau) d\tau - \frac{d}{d t} \int_{0}^{t} \mathscr{P}(\tau) d\tau = \frac{1}{2}\frac{d}{d t} \frac{T - t}{T + t}$

$ \mathscr{P}(t) = -\frac{d}{d t} \frac{T-t}{T+t}$

T = var('T')
t = var('t')
(-1/2 * diff((T-t)/(T+t), t)).full_simplify()

$ \frac{T}{T^{2} + 2 T t + t^{2}} $

よって $0<t<T$ について

$ \mathscr{P}(t) = \frac{T}{(T+t)^2}$

そして、$T<t<\infty$ の領域について。

$ t_2 = T (\frac{2}{1-p} - 1) $ より $p(t) = \frac{t - T}{T + t}$ なので

$ \int_{T}^{t} \mathscr{P}(\tau) d\tau = \frac{1}{2}\frac{t - T}{t + T}$

$ \int_{0}^{t} \mathscr{P}(\tau) d\tau - \int_{0}^{T} \mathscr{P}(\tau) d\tau = \frac{1}{2}\frac{t - T}{t + T}$

両辺を $t$ で微分すると

$ \frac{d}{dt} \int_{0}^{t} \mathscr{P}(\tau) d\tau - \frac{d}{dt} \int_{0}^{T} \mathscr{P}(\tau) d\tau = \frac{1}{2}\frac{d}{dt} \frac{t - T}{t + T}$

$ \mathscr{P}(t) = \frac{d}{dt} \frac{t-T}{t+T} = \frac{T}{(T+t)^2}$

T = var('T')
t = var('t')
(1/2 * diff((t-T)/(t+T), t)).full_simplify()

$ \frac{T}{T^{2} + 2 T t + t^{2}} $

よって、$T<t<\infty$ について

$ \mathscr{P}(t) = \frac{T}{(T+t)^2}$

どちらも同じ形なので、$t > 0$ 全域について

$ \mathscr{P}(t) = \frac{T}{(T+t)^2}$

であることがわかります。
$T=1$ のときの $ \mathscr{P}(t)$ をプロットすると：

t = var('t')
plot(1/(1+t)**2, t, 0, 10)

検算

$ \int_{0}^{\infty} \mathscr{P}(t) dt = 1$ の確認：

T = var('T')
t = var('t')
assume(T > 0)
integral(T / (T + t)**2, (t, 0, infinity))

→ $1$

$ \int_{t_1}^{t_2} \mathscr{P}(t) dt = p $ の確認：

T = var('T')
assume(T > 0)
p = var('p')
assume(p > 0)
assume(p < 1)
t1 = T / (2 / (1-p) - 1)
t2 = T * (2 / (1-p) - 1)
t = var('t')
integral(T / (T + t)**2, (t, t1, t2)).full_simplify()

→ $p$

まとめ

コペルニクス原理に基づく余命の予測方法に基づいて、余命の確率分布を以下のように計算しました：

$ \mathscr{P}(t) = \frac{T}{(T+t)^2}$

さらに $T=1$ を時間の単位にして無次元化すると

$ \mathscr{P}(t) = \frac{1}{(1+t)^2}$

となり、事象の余命の確率分布は概ね時間の２乗に反比例して減衰していくことがわかりました。