唯物是真 @Scaled_Wurm

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

毎日が天皇誕生日になるには何回天皇が交代する必要があるか(シミュレーション版)

今日は天皇誕生日ですが、以前「あと何回天皇が交代すれば毎日が天皇誕生日になるか(不謹慎)」の期待値を求める記事を書きました

毎日が天皇誕生日になるには何回天皇が交代する必要があるか - 唯物是真 @Scaled_Wurm
祝日と祝日の間に挟まれた日が、国民の休日で休みになるのを考慮していないという指摘を受けたので、今回はその場合の平均回数を求めます


厳密解をどうやって求めればよいか悩んでいたら「厳密解は諦めてシミュレーションでそれっぽい値を求めればよいのでは?」というアイディアをいただきました。ありがとうございます

国民の休日を考慮しないモデル

毎回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)

図が欲しいと言われたので図を追記


1200回交代するまでの途中時点の休日数を数えるのを、10000回繰り返して試した時の平均(中心の青線)と四分位数(全体でソートした時に25%と75%の位置に当たる数)の間の領域(黄色の部分)をプロット
f:id:sucrose:20141223134654p:plain
プロットに使ったソースコード

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()

参考

グラフ描く参考にした