唯物是真 @Scaled_Wurm

プログラミング(主にPython2.7)とか機械学習とか

毎日が天皇誕生日になるには何回天皇が交代する必要があるか

下記の記事を読んで、子供の頃に「天皇が交代するたびに祝日が増えれば毎日が祝日になるなぁ(不謹慎」と思っていたのを思い出しました

せっかくなので回数の期待値(平均)を計算してみた

2月29日(うるう日)を考えない場合

\(E(x)\)を残り\(x\)個の日付が埋まるまでの回数の期待値とします。
このときもちろん\(E(0)=0\)です
残り\(x\)個のとき、残りの日付が\(1\)個減る確率は\(\frac{x}{365}\)となります
この場合期待値は確率の逆数として考えてよいので、\(E(x) = E(x - 1) + \frac{365}{x}=\sum_{k=1}^{x}\frac{365}{k}\)となる
計算すると\(E(365)\approx 2364\)回となる。
以下の表によると日本の男性の平均寿命は\(79\)年なので(この先も男性とは限らないですが)、あと\(19\)万年弱ぐらい経過すれば毎日が天皇誕生日になることがわかる

追記

年数について重大な指摘がされていました

\(E(x)\)のグラフ(横軸が\(x\))
\(x\)が大きい時にはあまり回数を必要としない(縦の減少幅が小さい)が、\(x\)が小さくなると多くの回数が必要になる
f:id:sucrose:20140430234217p:plain

ちなみにこういうタイプの問題はクーポンコレクター問題という名前で知られています。

ソースコード

from pylab import *
s = 0
result = [s]
for k in xrange(1, 365):
    s += 365.0 / k
    result.append(s)

xlabel('#day')
xticks(xrange(1, 366, 50), xrange(365, 0, -50))
plot(xrange(1, 366), result[::-1])
show()

2月29日(うるう日)を考える場合

うるう日を加えると計算が複雑になります(少し合っているか不安

\(E(x, y)\)を残り\(x\)個の日付(うるう日以外)と残り\(y\)個の日付(うるう日)が埋まるまでの回数の期待値とします。
このとき\(E(0, 0)=0\)です
まずうるう日以外とうるう日になる確率を考えます。
うるう日は\(400\)年に\(97\)回出現するので、うるう日になる確率は\(p_l=\frac{97}{365\times400+97}\)
うるう日以外の日付が残り\(x\)個あるとき、\(x\)個のうちのいずれかが選ばれる確率は\(p_n(x)=\frac{400x}{365\times400+97}\)

うるう日を考えない場合、\(x>0\)について次の式が成り立つ
$$E(x, 0) = p_n(x)E(x - 1, 0)+(1-p_n(x))E(x, 0)+1$$
これを式変形すると上でやった計算と同様の式になる
$$E(x, 0) = E(x - 1, 0) + \frac{1}{p_n(x)}$$

次にうるう日を考えた場合、まず\(x=0\)について
$$E(0, 1) = p_lE(0, 0)+(1-p_l)E(0, 1)+1$$
これを式変形すると\(E(0, 1) = \frac{1}{p_l}\)
次に\(x>0\)について
$$E(x, 1) = p_lE(x, 0)+p_n(x)E(x - 1, 1)+(1-p_l-p_n(x))E(x, 1)+1$$
式変形すると
$$E(x, 1)=\frac{p_lE(x, 0) + p_n(x)E(x - 1, 1)+ 1}{p_l+p_n(x)}$$

以上より動的計画法を使えば簡単に計算できて、期待値\(E(365, 1) \approx 2693\)回となる

ソースコード

pl = 97.0 / (365 * 400 + 97)
pn = 400.0 / (365 * 400 + 97)

noleap = 0
leap = 1 / pl

for x in xrange(1, 366):
    noleap += 1 / pn / x
    leap = (pl * noleap + pn * x * leap + 1) / (pl + pn * x)

print leap