あと何年生きられるのか

このページでは、余命分布についてあれこれ考える（ここは書きかけ）。

死亡年齢の分布

ゼロ歳児が $x$から$x+dx$歳の間に死亡する確率が $\phi(x) dx$ と表されるとしよう。 確率密度関数 $\phi(x)$の実際の形は、乳幼児期にひとつピークがあって、その後老年で再びピークが生じるような、双峰形になるようだが、ひとまずその詳細にはこだわらないことにする。 そうすると、平均寿命 $E$ は、死亡した年齢の平均であるから、 $$ E = \int_0^\infty x \phi(x) dx $$ と書ける。 現在のように医療が発達する以前は、乳幼児の死亡率が非常に高かったため、その影響で平均寿命はずいぶんと短かかったようである。

同じことであるが、$x$歳まで生存している確率分布（生存関数）は、 $$ \ell(x) = \int_x^\infty \phi(x) dx = 1 - \int_0^x \phi(x) dx $$ と表すことができ、これを使うと、平均寿命は $$ E = \int_0^\infty \ell(x) dx $$ とも書ける。

icon-pc　練習：平均寿命

生存関数の積分が平均寿命を与えることを示しなさい。

余命の分布

次に $y$歳まですでに生きた人が、その先、$t$年後に死亡する確率を考えてみよう。

$y$歳まで生きるという事象をAとすれば、その確率$P(A)$は、上述の生存関数 $\ell(y)$ そのものである。

$y+t$歳で死亡するという事象を$D$とすれば、$y$歳まで生きて、かつ、$y+t$歳で死亡する確率 $P(D,A)$というのは、 0歳児が$y+t$歳で死亡する確率に他ならないので、 $$ P(D,A) = \phi(y+t) $$ となる。

ここでベイズの定理 $$ P(D|A) = \frac{P(D,A)}{P(A)} $$ を使えば、 $$ p(y+t|y) = \frac{\phi(y+t)}{\ell(y)} $$ となる。

以上から、$y$歳まで生きた人の平均余命$E_y$は $$ E_y = \int_0^\infty t \, p(y+t|y) dt = \frac{1}{\ell(y)} \int_0^\infty t \, \phi(y+t) dt $$ と表せる。 ここで、$\phi(x) = -\ell'(x)$であることを思い出して変形すると、平均余命は生存関数のみをつかって $$ E_y = \frac{1}{\ell(y)} \int_0^\infty \ell(y+t) dt = \frac{1}{\ell(y)} \int_y^\infty \ell(t) dt $$ で与えられる。

icon-teacher　解説: 指数分布と無記憶性

死亡年齢の分布が指数分布 $\phi(x) = \lambda e^{-\lambda x}$ に従うと過程してみよう。 生物にこれを当てはめると、幼い時期ほど死亡する可能性が高く、平均寿命が $\lambda$ 年であるような（ちょっと悲しい）場合であるが、 原子核の崩壊のような非生物な現象ではよくある例と言える。

生存関数は、これを積分して、$\ell(x) = \int_x^{\infty} \phi(t) dt = e^{-\lambda x}$ となる。

そうすると、$y$歳時の平均余命を、上記に従って求めると $$ E_y = \frac{1}{\ell(y)} \int_y^\infty \ell(t) dt = \lambda $$ であり、$y$への依存性が完全にキャンセルされ、平均寿命と同じになる。 つまり、指数分布で死亡（崩壊）するようなプロセスというのは、どの瞬間を出発点にしても全く同様に見える（無記憶性を持つ）ことがわかる。

生命表

人口動態の調査によって知ることができる基本的なデータは、ある時点での年齢毎の人口、および、それぞれの年齢での（1年間の）死亡人数である。 原理的に、ゼロ歳児の死亡年齢分布や生存関数は調査時点で知りようがない。 そこで、各年齢での死亡率が一定である（あった）という定常性を仮定して、生存関数を近似して得られた結果が、「生命表」として政府機関から公表されている。

生命表では、年齢 $x$ の定常人口を $L_x$ とし、$x$歳 の年間死亡数 $d_x$ から、「中央死亡率」（$(x+0.5)$歳での死亡率）$M_x$を $$ M_x \approx \frac{d_x}{L_x} $$ で見積もっている。 $x$歳の人口を$\ell_x$, $x+1$歳の人口を$\ell_{x+1}$とすれば、$x+0.5$歳の人口は線形補間を行って $$ L_x = \frac{\ell_x + \ell_{x+1}}{2} $$ と近似できるが、$\ell_{x+1} = \ell_{x} - d_x$ であるから $$ L_x = \frac{ 2 \ell_x - d_x}{2} $$ となる。

そのため、生命表では、上記の中央死亡率を使って、年齢$x$の生存者$\ell_x$が$x+1$歳を迎えるまでに死亡する確率（死亡率）$q_x$ を $$ q_x = \frac{d_x}{\ell_x} \approx \frac{M_x}{1 + \frac{1}{2} M_x} $$ と見積もっている。

また、生命表では各人口を10万人あたりの人数として掲載しているので、上記の$\ell_x$を10万で除した数が、生存関数 $\ell(x)$ に対応する。

残り年数のプロット

政府の統計データを公開しているサイト(e-Stat)から生命表のエクセルファイルをダウンロードすることができる。 e-Statの検索ページ(https://www.e-stat.go.jp/stat-search) にアクセスして、 絞り込み対象から「データセット」を選択して、キーワードを「生命表」にして検索すると、複数の年度の生命表がヒットするはずである。 ここでは、2024年の男性の簡易生命表のエクセルデータをダウンロードすることにする。 1life_table2024_male.xlsx という名前のファイルが出来ているはずである。

このエクセルファイルを読み込んで、年齢毎の余命をプロットするPythonコードの例を以下に示す。 生命表には各年齢の平均余命がすでに含まれているが、検算も兼ねて、ここでは敢えて生存数（生存関数）から計算するようにしてみた。

このコード中では、PandasとMatplotlibに加え、Pandasでエクセルファイルを読み込むためにOpenPyXLを使っているので、もしまだなら、 pipコマンドで追加しておく必要がある。 


import pandas as pd
import matplotlib.pyplot as plt

file_path = '1life_table2024_male.xlsx'

# Excelファイルの読み込み
# header=2: 3行目（インデックス2）をヘッダーとして読み込みます
# engine='openpyxl': .xlsxファイルを読み込むためのエンジンを指定
try:
    df = pd.read_excel(file_path, header=2, engine='openpyxl')
except FileNotFoundError:
    print(f"エラー: ファイル '{file_path}' が見つかりません。パスを確認してください。")
    exit()

# --- データの前処理 ---

# 1. 「年齢」列の欠損を除外
df = df.dropna(subset=['年齢'])
df['年齢'] = df['年齢'].astype(str)

# 2. 「年齢」に「年」が含まれる行のみを抽出（「0週」「2月」などの細かい内訳を除外）
df_plot = df[df['年齢'].str.contains('年')].copy()

# 3. 「年齢」から「年」の文字を削除して数値に変換する関数
def clean_age(val):
    try:
        # “105年” -> 105
        return int(val.replace('年', ''))
    except ValueError:
        return None

df_plot['age_num'] = df_plot['年齢'].apply(clean_age)
df_plot = df_plot.dropna(subset=['age_num'])
df_plot['age_num'] = df_plot['age_num'].astype(int)

# 4. 「生存数」を数値に変換した上で、10万人あたりからひとりあたりに変換。
df_plot['survival_rate'] = pd.to_numeric(df_plot['生存数'], errors='coerce') / 100000

# 5. 生存率を微分して、死亡年齢分布に変換
df_plot['age_at_death'] = - df_plot['survival_rate'].diff().fillna(0)

# 6. 余命を計算 
df_plot['life_expectancy'] = df_plot['survival_rate'].iloc[::-1].cumsum().iloc[::-1] / df_plot['survival_rate']

# --- プロット ---

plt.figure(figsize=(10, 6))

# 横軸：年齢(age_num)
# plt.plot(df_plot['age_num'], df_plot['survival_rate'], label='Survival Rate')
# plt.plot(df_plot['age_num'], df_plot['age_at_death'], label='Age distribution of deaths')
plt.plot(df_plot['age_num'], df_plot['life_expectancy'], label='Life Expectancy')

# グラフの装飾
plt.title("2024 Life Table for Males")
plt.xlabel('Age')
# plt.ylabel('Survival')
# plt.ylabel('Age at death')
plt.ylabel('Life Expectancy (Years)')
# plt.yscale('log')
plt.grid(True)
plt.legend()

# 表示
plt.show()

プロットを眺めると、0歳の時点での平均余命は81歳強であるが、年齢を経るに従って「残り時間」は単調に減少する様子が分かる。 60歳を越えたあたりから、グラフの勾配が少しずつなだらかになるのは、そこまで死なずに生きながらえた効果（しぶとさ？）と解釈することもできるだろう。 このページを制作している時点で、筆者は65歳であるが、そうすると、残り時間の平均は約20年ほどということになる。 この時間を大切にして生きたいものである。

同じデータとコードを使って、ゼロ歳児がどの年齢で死亡するかの分布$\phi(x)$をプロットしてみたのが以下である。 生まれてすぐの時期を過ぎると一旦減少するものの、成人を迎える頃に再び増加し、その後はしばらく横ばい状態が続く。 その後、大厄の頃を起点として老年期にはいると死亡率は急増し（80歳くらいまではほぼ指数関数的）、米寿の頃にピークを迎えている。

生命表の死亡率（$x$歳の人口のうちどれだけが1年間で死亡するか）を片対数プロットすると、以下のように、中年期以降、驚くくらい綺麗に指数関数的に増加していることがわかる。 試みに、2024年の男性の生命表について40歳から100歳までの範囲で死亡率を指数関数$d(x) = A \exp(b x)$ でフィッティングしてみると、$A = 1.32\times 10^{-5}$, $b = 0.103$であった。

icon-pc　練習：フィッティング

生命表のデータの死亡率の変化を指数関数でフィッティングするコードを作成し、パラメータを求めてみなさい。


老化のプロセス

事実として、壮年期から老年期までの一定の期間、死亡確率は年齢と共に指数関数的に増大しているように見える。 そこで、年齢$x$の生存者が、その後$x+\Delta x$までの間に死亡する確率密度を $$ q(x) = A e^{bx} $$ と仮定してみる。これを使うと、$T$歳まで生存する確率は、「死なずに済む確率」の積 $$ \ell(T) = \prod_{i}^{x_i = T - \Delta x} \left( 1 - A e^{b x_i} \Delta x \right) $$ となる。

その対数を取ると $$ \log \ell(T) = \sum_i \log\left( 1 - A e^{b x_i} \Delta x \right) \approx - \sum_i^{x_i=T} A e^{b x_i} \Delta x $$ であるので、和を積分で書き直すと（0歳から老化が始まるとして） $$ \log \ell(T) = - \int_0^T A e^{bx} dx =\frac{A}{b} \left( 1 - e^{bT} \right) $$ すなわち、生存関数は $$ \ell(T) = \exp\left( \frac{A}{b} \left( 1 - e^{bT} \right) \right) $$ と書ける。

わざわざ微分を取るまでもなく、その定義から $$ \phi(x) = \ell(x) \, q(x) $$ である。

これを用いると、死亡年齢の分布は $$ \phi(x) = - \frac{d \ell(x)}{d x} = A \exp\left[ b x - \frac{A}{b} (e^{bx} -1) \right] $$ となり、ゴンぺルツ分布として知られている。

いくつかのパラメータで、ゴンペルツ分布をプロットしてみた例（コードとグラフ）を以下に示す：


import numpy as np
import matplotlib.pyplot as plt

def gompertz_pdf(x, A, b):
    """
    ゴムペルツ分布の確率密度関数 phi(x)
    x: 年齢
    A: 初期死亡率パラメータ (Scale parameter)
    b: 老化率パラメータ (Shape parameter)
    """
    
    # ハザード関数部分: q(x) = A * e^(bx)
    hazard = A * np.exp(b * x)
    
    # 生存関数部分: ell(x) = exp( - (A/b) * (e^(bx) - 1) )
    survival = np.exp(-(A / b) * (np.exp(b * x) - 1))
    
    return hazard * survival

# 1. パラメータの設定
# 日本人男性に似せたパラメータセットと、比較用のセットを用意
params = [
    {"A": 1.318e-05, "b": 0.103, "label": "Japanese male-like(b=0.103)"},
    {"A": 1.318e-05, "b": 0.15,  "label": "Fast aging (b=0.15)"},
    {"A": 0.002,  "b": 0.05,  "label": "Slow aging with high infant death"},
]

# 2. x軸（年齢）のデータ作成
x = np.linspace(0, 120, 1000)  # 0歳から120歳まで

# 3. プロットの作成
plt.figure(figsize=(10, 6))

for p in params:
    y = gompertz_pdf(x, p["A"], p["b"])
    
    # 最頻値（モード）の計算（参考用）: x_mode = (1/b) * ln(b/A)
    # ただし定義域内にある場合のみ
    mode_val = (1/p["b"]) * np.log(p["b"]/p["A"])
    
    plt.plot(x, y, linewidth=2.5, label=p["label"])
    
    # 最頻値にマーカーを表示（グラフ内にある場合）
    if 0 <= mode_val <= 120:
        y_mode = gompertz_pdf(mode_val, p["A"], p["b"])
        plt.plot(mode_val, y_mode, 'o', color='black', markersize=5)

# グラフの装飾
plt.title(r"Gompertz distribution $\phi(x) = A \exp\left[ bx - \frac{A}{b} (e^{bx} -1) \right]$", fontsize=16)
plt.xlabel("Age $x$ (Years)", fontsize=14)
plt.ylabel("Death rate $\phi(x)$", fontsize=14)
plt.grid(True, which="both", linestyle="--", alpha=0.7)
plt.legend(fontsize=12)
plt.xlim(0, 120)
plt.ylim(bottom=0)

# 表示
plt.tight_layout()
plt.show()

死亡する確率が年齢と共に指数関数的に増大するというのは、（特に私のような高齢者には）少なからずショッキングな事実である。 加齢と共に何らかの不具合が自己触媒的に不具合の連鎖を生じさせているには違いないが、 健康診断の結果として顕在するかなり前からそうしたカスケードが進行しているとすると、もうどうしようもあるまい。

アルツハイマー病の診断のため、原因物質と目されているアミロイドβやタウタンパクを検出する手法が開発されているようだが、発症前の患者には検査を行わないようである。 これらのタンパクが検出されたからといって、発症までには一般にかなりの（猶予）期間があるらしく、「あなたはアルツハイマー病になる可能性が大ですよ」 と告げることが、その人にとってプラスになるかどうか、判断が難しいからだそうである。 同様に、遠からず、かなり早い時期に人生の残り時間が高い精度で予測できるようになるかもしれないが、「あなたの残り時間は21年です」と告げられても、黒澤明の映画『生きる』のようにはならない（できない）だろう。

モノの寿命

年齢 $x$ の人がさらに $dx$ 年のうちに死亡する確率（死亡率）を $q(x) dx$ としたとき、老年期にかけて$q(x)$は指数関数でよく表されることをみてきた。 一方、ヒトではなく、モノの故障や破壊については $$ q(x) = A x^{m-1} $$ のようなべき乗関数がよく当てはまると言われている。ここで、べき指数を1だけ引いているのは慣習による。 パラメータ$m$の値は、事象毎にかなりの幅があるようである。 $q(x)$は一般に危険率(hazard rate)と呼ばれている。 危険率が指数関数ではなくべき乗関数となるのは、故障や破壊が、カスケード的というよりは、散発的な事象によって突然生じることによる（例えば１個の部品の不良で全体がお釈迦になる）、という違いを反映している。

前節までと同様の計算を行うと、生存関数の対数が $$ \log \ell(T) = - \int_0^T A x^{m-1} dx = - \frac{A}{m} T^{m} $$ となるから、 $$ \ell(x) = \exp\left(- \frac{A}{m} x^{m} \right) $$ 製造から故障までの時間の分布は $$ \phi(x) = - \frac{d \ell(x)}{dx} = A x^{m-1} \exp\left( - \frac{A}{m} x^{m} \right) $$ となり、これはワイブル分布と呼ばれている。 故障率が経過年数に依らない（すなわち$m=1$の）場合、分布は単純な指数分布に帰着される。

もともとW. Weibullがこの分布を案出したのは、 故障時間ではなく、脆性物質の破壊強度の実験を説明するためだったようである。 多数の（同様の）試料を用意して、応力を加えて破壊実験を繰り返し、破壊時の応力の分布を調べると、それがワイブル分布で良く表現できることを彼は見出したのであった。

物質の強度は、それを構成する原子や分子の結合力から見積もられる値よりも遥かに小さい（弱い）ことが分かっていて、その原因は、 物質内に無数に存在する欠陥（ミクロな亀裂）にあると考えられている。 物質全体に徐々に応力を加えると、これらの欠陥のうちどれかが「活性化」することで全体が壊れる、というわけである。

物質の強度を議論する場合には、$x$は単純な時間経過ではなくて、試験で加える応力の系列であると解釈すればよい。 すなわち、0から徐々に応力を加えた場合、$x$では破断しなかった試料が、$dx$だけさらに余分の応力を加えたら壊れる比率が危険率$q(x)$に対応する。

ヒトとモノにとっての時間

死亡率（ハザード関数）の形は、ヒトとモノとでは異なるように見えるわけだが、そもそも、「時間」という変数が同じモノサシで良いかどうか、疑問が残る。 どちらも「物理的時間」の経過に伴って劣化が進行するには違いないが、その進み具合こそがその系にとって固有の「時計」である、という解釈も可能かもしれない。 そうすると、ヒトにとっての時計は（老年期は）物理時間を指数関数で変換したようなスケールで進み、 モノの時計は物理時間をべき乗関数で変換したように進む、ということになる。 歳を取ると時間が速く過ぎるように感じるのは、実際に、「時計」の進みが速いのかもしれない。