今日は天皇誕生日ですが、以前「あと何回天皇が交代すれば毎日が天皇誕生日になるか(不謹慎)」の期待値を求める記事を書きました
毎日が天皇誕生日になるには何回天皇が交代する必要があるか - 唯物是真 @Scaled_Wurm
祝日と祝日の間に挟まれた日が、国民の休日で休みになるのを考慮していないという指摘を受けたので、今回はその場合の平均回数を求めます
さらに、挟まれた日が国民の休日になるというのを考えると、もっとずっと複雑になるな。(考える気はない)http://t.co/AuibRNF969
— Hiroshi Manabe (@takeda25) 2014, 4月 30
厳密解をどうやって求めればよいか悩んでいたら「厳密解は諦めてシミュレーションでそれっぽい値を求めればよいのでは?」というアイディアをいただきました。ありがとうございます
@Scaled_Wurm ああ、ここでシミュレータと言ったのは、ある日をランダムで祝日にする→国民の休日になってたらそれを反映する→ある日をランダムで祝日にする、を全部祝日で埋まるまでシュミレーションするだけです。何回かかったかは平均値とれば大体わかるかと
— りすお@社二病 (@risuosan) 2014, 12月 23
国民の休日を考慮しないモデル
毎回1年365日のいずれかの日が乱数で選ばれて天皇誕生日になるときの、すべてが天皇誕生日になって終了するまでの回数の平均値を求めます
下記のコードを実行したら10000回の平均はおよそ2366回になりました(実行するたびに変わります)
これは前回の記事で求めた理論値の2364回にかなり近い値です
このことからシミュレーションでも近い値がでると考えられるので、国民の休日を考慮したシミュレーションを試してみます
import random result = [] for z in xrange(10000): N = 365 count = 0 days = [0] * N t = 0 while count < N: i = random.randint(0, N - 1) if days[i] == 0: days[i] = 1 count += 1 t += 1 result.append(t) print 'mean:', float(sum(result)) / len(result)
国民の休日を考慮したシミュレーション
下記のコードを実行したら100000回の平均はおよそ1174回になりました(実行するたびに変わります)
これは国民の休日を考慮しない場合のモデルの半分弱の回数になります
このことから一つの年号が3,40年続くとすると、3,4万年経てば毎日が休日になることがわかります
import random result = [] for z in xrange(100000): N = 365 count = 0 days = [0] * N t = 0 while count < N: i = random.randint(0, N - 1) if days[i] == 0: days[i] = 1 count += 1 if days[i - 2] == 1 and days[i - 1] == 0: days[i - 1] = 1 count += 1 if days[(i + 2) % N] == 1 and days[(i + 1) % N] == 0: days[(i + 1) % N] = 1 count += 1 t += 1 result.append(t) print 'mean:', float(sum(result)) / len(result)
図が欲しいと言われたので図を追記
3万年かぁ…遠いなぁ。横軸が時間、縦軸が休日の数(幅付き)のグラフがほしいところ。
— berobero (@berobero11) 2014, 12月 23
1200回交代するまでの途中時点の休日数を数えるのを、10000回繰り返して試した時の平均(中心の青線)と四分位数(全体でソートした時に25%と75%の位置に当たる数)の間の領域(黄色の部分)をプロット
プロットに使ったソースコード
import random import collections import numpy as np import matplotlib.pyplot as plt result = [[] for t in xrange(1200)] for z in xrange(1000): N = 365 count = 0 days = [0] * N for t in xrange(1200): i = random.randint(0, N - 1) if days[i] == 0: days[i] = 1 count += 1 if days[i - 2] == 1 and days[i - 1] == 0: days[i - 1] = 1 count += 1 if days[(i + 2) % N] == 1 and days[(i + 1) % N] == 0: days[(i + 1) % N] = 1 count += 1 result[t].append(count) data = np.array(result) mean = data.mean(axis=1) quantile1 = np.percentile(data, 25, axis=1) quantile3 = np.percentile(data, 75, axis=1) plt.xlabel('#change') plt.ylabel('#holiday') plt.plot(xrange(1, 1201), mean, lw=2) plt.fill_between(xrange(1, 1201), quantile1, quantile3, facecolor='yellow', alpha=0.5) plt.grid() plt.show()