苗字の増減を考えるモデルを通して物理の考え方を学んでみよう


# イントロ
少し前に、こんなツイートを見かけました。

https://x.com/mathlava/status/1755467188826452092

調べてみたら結構面白そうだったので、今回はこれを扱ってみます。

理論をざっくり学んだら、シミュレーションもやってみます。卒研生は、この記事を通してモデルとシミュレーションの感覚を掴んでみてください。対象が実験で再現可能なものでないので、物理ではないですが、我々のよく用いる、現象論的アプローチの感覚を身につけやすいかなと思います。

# モデル
「ある家系において、その姓を持つ子孫の数が、世代を経るごとにどのように増減していくか」という問題を確率過程論的に扱う、Galton-Watson過程というモデルを学びます。中でも今回は、最もシンプルなものを扱います。

### 仮定と定義
- 結婚すると、夫婦のうち一方が姓を継続して名乗り、もう一方が相手の姓に変更
    - 変わらない方を継承者と呼ぶ
- ある継承者が一生のうちにもうける子のうち、特に継承者になる子の数が、ある確率分布に従うと仮定
- 全員が共通の確率分布に従い、同時に幾つかのメンバーに分裂
- 分裂は過去の履歴や、他のメンバーの振る舞いとは無関係（Malkov過程）

## 予想

1組の親から生まれる継承者数の平均値が1を下回ったら、その苗字はいつか消えてしまいそうだと予想できます。これをぼんやり意識しながら、以下を読んでみてください。

## 導出
### 第j世代までに絶滅する確率
ある世代の継承者が、次世代の継承者を$k$人もうける確率を$p_k$を
```math
\begin{align}
p_k\ge0\\
\sum_{k=0}^{\infty} p_k = 1
\end{align}
```
と定義する（$p_0$は次世代の継承者が生まれない確率）。このとき、任意の1人の継承者が、生涯にもうける次世代継承者数の期待値$m$は
```math
\begin{align}
m = \sum_{k=0}^{\infty}k\cdot p_k
\end{align}
```
である。

また、j世代までに姓が絶滅する確率を$d_j$、いつか絶滅する確率を$d$とすると
```math
\begin{align}
d = \lim_{j\rightarrow\infty} d_j \tag{1.1}
\end{align}
```
と書ける。

では、第$j$世代までに姓が絶滅する確率を考えていく。
- 第0世代の1人から始まることを考えるので、$d_0=0$
- 第1世代までに絶滅する確率は、第0世代が1人も継承者を設けない確率で、$d_1 = p_0$
- 第$j+1$世代までに姓が絶滅するとき
これは、第1世代に$k$人の継承者が存在していたとする(これが起こる確率は$p_k$）と、第1世代の$k$人の継承者を起点とする各々の家系が全て第$j$世代目までに絶滅することに他ならない

![5CF1FE9F-9680-411C-ADD1-2C3140D53B90_4_5005_c.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3737640/7c033111-b066-8061-3d34-b73e3bc8e0cc.jpeg)


Malkov過程を考えているので、それぞれ独立な$k$個の家系が、第$j$世代目までに絶滅する確率は$d_j$であり、$k$個の家系が全て第$j$世代目までに絶滅する確率は、$(d_j)^k$である。
以上より、第$j+1$世代までに姓が絶滅する確率$d_{j+1}$は、
```math
\begin{align}
&d_1 = p_0\\
&d_{j+1} = \sum_{k=0}^{\infty}p_k (d_j)^k \quad (j\in \mathbb{N})
\end{align}
```
となる。

ここで、確率分布${p_k}$の積率母関数(moment generating function)$f(t)$は
```math
f(t) = \sum_{k=0}^{\infty}p_k \cdot t^k \tag{1.2}
```
で定義されていたことを思い出すと、
```math
\begin{align}
&d_{j+1} = f(d_j)\\
&d_0 = 0 \tag{2.1}\\ 
&d_1 = f(d_0) = p_0
\end{align}
```
と書ける。これを図示すると以下

![B6E79279-105F-4DFE-8557-23EFCC5EFB4F_1_201_a.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3737640/7e51b29c-6847-85e7-db31-cb21289c58b4.jpeg)


のようになっていて、振る舞いは以下のようである。
1. $d_0=0$
1. $d_1=f(d_0)$はy軸の切片
1. $y=t$の直線のy座標が$d_1$になる位置から真下におろして、$x=d_1$を見つける
1. $f(d_1)$を見つける
1. 以下繰り返し

これを繰り返すと、$d_i$は$y=f(t)$と$y=t$の最初の交点へ向かい、それを超えないことが分かる。


### 幾何分布に従う次世代継承者数
ここからは、$p_k$の分布が
```math
\begin{align}
p_k = (1-a)a^k \quad (k = 0, 1, 2 \cdots, 0<a<1)\\
a: \text{次世代継承者を1人もうけると選択し続ける確率}
\end{align}
```

で与えられる場合について考えてみましょう。この分布の"お気持ち"としては、

「ある継承者の次世代継承者が$k$人である」$\iff$ 「継承者をもうけると$k$回連続で選択[^2]→もうけるのをやめる」

なので、上式の右辺のように書けます[^3][^4]。

[^3]:$a$が小さい方が早く絶滅することになるのは、感覚に合う
[^4]:$k=0$では$p_0=1-a$だが、$p_0$は次世代継承者を一人ももうけなかった場合に対応する確率なので、うまく表現できている

[^2]:双子などは考慮しない

ここで、ようやく"物理"が入りました。我々に馴染みのある表現では、構成方程式みたいなものでしょうか。逆に卒研生は、このように「性質をうまく切り取って表現したものが構成方程式」というイメージを持つと良いのではないかと思います。

この確率分布を仮定した場合の、絶滅確率について考えてみましょう。$p_k = (1-a)a^k $より、
```math
\begin{align}
f(t) &= \sum_k p_k t^k\\
&= \sum_{k=0}^{\infty} (1-a)a^k t^k\\
&\downarrow a\le 1, t\le 1\\
&= (1-a)\frac{1}{1-at}
\end{align}
```
となり、
```math
\begin{align}
&d_0 = 0\\
&d_1 = p_0 = 1-a \tag{3.1}\\ 
&d_j = f(d_{j-1}) = \frac{1-a}{1-ad_{j-1}}
\end{align}
```
が今回の解くべき方程式として得られます。
この漸化式は有名な形ではなく、一見解析解は存在しなさそうですが、テキストには
```math
\begin{align}
&d_j = (1-a)\frac{(1-a)^j-a^j}{(1-a)^{j+1}-a^{j+1}} \quad \left(a \ne \frac{1}{2}\right)
 \tag{3.2}\\
&d_j = \frac{j}{j+1} \quad \left(a = \frac{1}{2}\right)
\end{align}
```
が解として示されています。実際代入してみると、式(3.1)を満たします。


# 結果
式(3.2)で得られた解をプロットして、以下を得ます。

![Galton_ana.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3737640/5017fb77-c63d-7440-e0e2-2c7e9b931815.png)



確かに、aが小さいと必ず絶滅しそうです。

# シミュレーション
解析解が見つからない場合は、数値的に（近似的に）解くというアプローチがよく用いられます。
今回の例だと、式(3.1)の漸化式へ実際に値を代入していけばある$j$に対する$d_j$を具体的に求めることができます。
卒研生は、プログラミングに慣れていないと思うので、ある$a, j$に対して$d_j$を出力するようなコードを書いてみましょう。数値的なアプローチが難しければ、解析解を出力してもOKです。

<details><summary>ヒント</summary>

数値解を求めるpythonのコードを置いておくので、C++で書き直してみてください。
```py
def d_j(j, a):
    if j==0: return 0
    elif j==1: return 1-a
    else: return (1-a)/(1-a*d_j(j-1, a))

# a=0.2の時、第10世代までに絶滅する確率を出力
print(d_j(10, 0.2))
```

</details>

<details><summary>実装例</summary>

```py
# 1種系モデルにおける解析解と数値解をプロット

import matplotlib.pyplot as plt

# 定数の設定
n = 50     # 追いかける世代数
a_list = [0.2, 0.5, 0.7]     # 試したいaの候補

# 数値的に解く
def d_j(j, a):
    if j==0: return 0
    elif j==1: return 1-a
    else: return (1-a)/(1-a*d_j(j-1, a))

# 解析解
def analytical(j, a):
    if a==0.5: return j/(j+1)
    else: return (1-a)*((1-a)**j - a**j)/((1-a)**(j+1) - a**(j+1))

# グラフの設定
plt.rcParams['figure.autolayout'] = True
plt.rcParams['figure.figsize'] = [5, 4]
plt.rcParams['font.size'] = 12
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['axes.linewidth'] = 1.0
plt.rcParams['lines.markersize'] = 6
fig, ax = plt.subplots()
ax.set_xlim(0, n)
ax.set_ylim(0, 1)
ax.set_xlabel("$j$", fontsize=18)
ax.set_ylabel("$d_j$", fontsize=18)
ax.set_title('changes in d_j')

# selected = "ana"
selected = "num"

# 実際にシミュレーション
for a in a_list:
    x = list(range(n))
    y = [d_j(i, a) for i in range(n)]
    if selected == "num":
        ax.plot(x, y, label='a={}'.format(a))
        
    
    ya = [analytical(i, a) for i in range(n)]
    if selected == "ana":
        ax.plot(x, ya, label="a={}".format(a))


ax.legend(loc='lower right')
plt.savefig(f"Galton_{selected}.png")
```
</details>

数値的に解くと、以下のような図が得られました。先ほどの図と、目では違いが分からない程度に一致しています。
このように、数値的に解くのも有効なアプローチであるということが実感してもらえたのではないでしょうか。

![Galton_num.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3737640/85374201-2b10-63b2-6047-c84ffe7c61e0.png)



# モデルの改良
今回は、非常にシンプルなモデルを扱いましたが、例えば継承者の性別を取り入れたりすることで、より現実的なモデルにすることも可能です。
また、苗字の推移以外にも、例えば中性子の核連鎖反応[^5]などもGalton-Watson型分岐過程が観察されるそうです[^6]。これはまさしく物理ですね。

[^5]: 原子核と中性子が衝突した結果、ランダムな個数の中性子が放出される。放出された中性子は、それぞれ他の原子核に衝突し、新たな中性子の放出を引き起こす。
[^6]: 参考文献に(Feller, 1957)としか書いてなかったのですが、「An Introduction to Probability Theory and Its Applications」かなと思います。

# まとめ
今回は、Galton-Watson型分岐過程を扱いました。
なんとなく現象論的アプローチのイメージが掴めたでしょうか。同時に、シミュレーションの面白さを少しでも感じてもらえたら嬉しいです。

---
## おまけ
素朴な予想通り、次世代継承者数の期待値$m$が$m\le1$なら$d=1$であり、$m>1$なら$d<1$であることを確かめます。

### 最終絶滅確率
$j \rightarrow \infty$では$d_j = d_{j+1}$なので、式(1.1)と式(2.1)より、
```math
\begin{align}
d = f(d) \tag{4.1}
\end{align}
```
となります。
また、上で見たように、$d$は$y=f(t)$と$y=t$の最初の交点です。ここで、$y=t$と$y=f(t)$は$0 \le t \le 1$において少なくとも１つ、高々２つの交点を持ちます(Appendix参照)。
それぞれ図示すると、以下のようになっています。

![DDF1AFFB-51C9-4B40-BD32-9788BE8CE557_1_201_a.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3737640/9b343651-5db2-eeb5-9db1-7746db3fda0a.jpeg)


(ⅰ)
最初の交点を見ると、$d=1$
また、$f'(1) \le 1 \iff m \le 1$

以上より、期待値$m \le1$ $\implies$ $d=1$

(ⅱ)
(ⅰ)と同様に、$m>1$ $\implies$ $d<1$

[まとめ]

1人の継承者がもうけると期待される次世代継承者の人数$m$について、
```math
\begin{align}
&m \le1 \implies d=1  \text{ （必ず消滅）}\\
&m>1 \implies d<1 \text{　（永続確率は$1-d$）}
\end{align}
```
であることがわかりました。これは、最初の素朴な予想に一致します。


ちなみに、$d<1$の時、任意の継承者から始まる家系の分布が永遠に続く確率は$1-d$で与えられます。これは、任意の継承者を祖先と考えれば理解できそうです。

### 幾何分布における最終絶滅確率
$f(t) = \frac{1-a}{1-at}$の場合、1階微分は$f'(t) = \frac{a(1-a)}{(1-at)^2}$なので、
```math
m = f'(1)=\frac{a}{1-a} \tag{4.2}
```
です。

また、式(4.1)より、
```math
\begin{align}
&y = f(t) = \frac{1-a}{1-at}\\
&y = t
\end{align}
```
の交点が$d$を与えるので、これを解くと、
```math
\begin{align}
t = 1, \frac{1-a}{a}
\end{align}
```
が交点と分かります。$t<1$より、
```math
\begin{align}
&d = 1 \quad (a \le \frac{1}{2})\\
&d = \frac{1-a}{a} \quad (a > 1)
\end{align}
```
を得ます。これを図示すると

![Galton_d-a.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3737640/10ab7e56-fdf7-adf8-c4e8-4197d6a0708a.png)

のようになっており、$a$が大きくなると絶滅確率が小さくなるという予想に合います。

また、式(4.2)より

```math
\begin{align}
&d = 1 \quad (m \le 1)\\
&d = \frac{1}{m} \quad (m > 1)
\end{align}
```
と書くこともできます。

$a=\frac{1}{2}$で$m=1$だったので、これが境界になるというのは、やはり最初の素朴な予想に一致します。

## Appendix
積率母関数は、以下のような特徴を持つ。
```math
\begin{align}
&f(0) = p_0 \ge 0\\
&f(1) = \sum_{k=0}^{\infty} p_k = 1\\
    &\left. \frac{1}{n!}\frac{d^n f}{dt^n}\right|_{t=1} = p_n  \\
&\left. \frac{df}{dt}\right|_{t=1} = \sum_{k=0}^\infty k \cdot p_k = m \quad \text{:確率変数の期待値}
\end{align}
```
母関数が定義される範囲において、$p_k\ge 0$である。よって$df/dt$、$d^2f/dt^2$は必然的に非負である。よって$y=f(t)$は下に凸の増加関数である。この時、下図を見ると、少なくとも一つ、高々２つの交点を持つことは明らか。
![773EC6FE-6209-4610-8DE6-1135773E9CB5_4_5005_c.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3737640/005ce947-6fbf-c89b-77fd-09a69e9e321b.jpeg)


# 参考文献
- 佐藤葉子, 瀬野裕美 姓の継承と絶滅の数理生態学 京都大学学術出版