この記事について

この記事は、プログラミング初心者が覚えたてのモンテカルロ法を使ってみる話です。
登場するコードは超単純なのですが、長年の疑問を自分の手で解決できたことが非常に嬉しかったので記事にしました。
最終的に得られた確率(かなり綺麗な値)の数値的な意味については考察中なので、確率が得意という方は是非教えて下さい。

「時計」とは

一人用トランプゲームの一種で、ルールはとても単純です。
なお、プレイヤーの意志を挟む余地がない完全ルールベースの運ゲーなので、楽しさは保証できません。

[準備]
ジョーカーを抜きよく切った52枚のトランプを、画像のような形で13の山に分けます。
そしてこれを時計の文字盤になぞらえ、それぞれの山に数字を割り振ります。
余ってしまう中心の山には、ご都合主義で13を割り振ってください。
　　

[ルール]
1.　13の山の1番上のカードをめくります。
2.　めくったカードの番号を確認し、今度はその番号の山の1番上のカードをめくります。
3.　続けられる限り、2の作業を続けます。
4.　4枚目の13が出た時点で、それ以上めくるカードがなくなりゲーム終了となります。
5.　すべてのカードをめくることができたら成功、それ以外は一律失敗とします。

何が楽しいの？と思う方もいらっしゃると思いますが、13が出るたび寿命が縮んでいくようなこのスリル、意外と中毒性があります。
私自身、テスト前夜の現実逃避でトランプに手を伸ばし、「1回上がったら再開しよう」などと始めてしまったせいで今この記事を書くに至っています。

紙とペンと脳味噌で成功確率を考えてみた

直感的に、成功確率は1/13かなと思っていました。
イメージとしては、長さ52のリストの最後の要素が13なら成功、それ以外は失敗、という感じです。
しかし考えれば考えるほど、そう単純な話では収まらないような気がしてきました。
先ほどの例えでいえば時計は、リスト内を遷移する順番がリストの要素によって決定されるという仕組みです。
この複雑さ、紙とペンと脳味噌では少々荷が重いと判断し、習いたてのモンテカルロ法を使ってみることにしました。

モンテカルロ法とは

モンテカルロ法を一言で説明するなら、
実際に何回もシミュレーションしてみて、成功する確率を調べる手法です。
今回は、Pythonに時計を何回も何回もやらせて、成功回数/試行回数を計測します。
円周率の計算などで有名な手法なので、ご存じの方が多いかもしれません。

退屈なことはPythonにやらせよう

まず、カードをシャッフルし、52枚を並べる関数を作ります。

def shuffle(): #52枚のカードを並べる
  shuffle =[1,2,3,4,5,6,7,8,9,10,11,12,13]
  shuffle += shuffle + shuffle + shuffle
  random.shuffle(shuffle)
  return shuffle

ここからはゲームが始まるわけですが、
プレイヤーの意志を挟む余地がない完全ルールベースの運ゲーなので実装も簡単でした。
先に必要な情報をまとめておきます。

先ほど作った52枚分のリストとカードの対応は、
インデックス：山
0~3：1の山
4~7：2の山
...
(N-1)*4 + i (0≦i≦3)：Nの山
とします。

また、「開く」という行為はリストの要素を0に書き換えることで表します。

def tokei():
  cards = shuffle() #52枚のランダムな配列を用意
  flag =0
  number = 13 #最初に開く１３
  while flag == 0:
    #print(number) 
    if cards[(number-1)*4] != 0: #その番号が初めて開かれる場合
      next = cards[(number-1)*4]
      cards[(number-1)*4] = 0
      number = next
    elif cards[(number-1)*4+1] != 0: #その番号が２回目に開かれる場合
      next = cards[(number-1)*4+1]
      cards[(number-1)*4+1] = 0
      number = next
    elif cards[(number-1)*4+2] != 0: #その番号が３回目に開かれる場合
      next = cards[(number-1)*4+2]
      cards[(number-1)*4+2] = 0
      number = next
    elif cards[(number-1)*4+3] != 0: #その番号が４回目に開かれる場合
      next = cards[(number-1)*4+3]
      cards[(number-1)*4+3] = 0
      number = next
    else:                            #13が4枚出ちゃった場合or終わった場合
      flag =1
  return sum(cards)

いよいよモンテカルロ法の登場です。

def Monte_Carlo(n):
  N = 0
  for i in range(n):
    if tokei() == 0:
      N += 1
  return N/n

さしあたりnを10000で実行してみます。

Monte_Carlo(10000)

0.0787

それっぽい値が出てきました。
統計は苦手ですが、100万回くらい試行すれば十分でしょうか。

print("100回の試行では "+ str(Monte_Carlo(100)))
print("1000回の試行では "+ str(Monte_Carlo(1000)))
print("10000回の試行では "+ str(Monte_Carlo(10000)))
print("100000回の試行では "+ str(Monte_Carlo(100000)))
print("1000000回の試行では "+ str(Monte_Carlo(1000000)))
print("1/13 = " + str(1/13))

実行結果は以下の通りです。

100回の試行では 0.07
1000回の試行では 0.083
10000回の試行では 0.0736
100000回の試行では 0.07664
1000000回の試行では 0.076582
1/13 = 0.07692307692307693

統計は苦手ですが、これは大体1/13に収束しそうだと言ってもいいですかね。
そういうことにしましょう。

総括

自分の手でコードを書いて結論が出せて、滅茶苦茶楽しかったです。
whileが止まらなくなったり1枚目めくってゲームが終わっちゃったり、2枚目めくってゲームが終わっちゃったりしながらも、ググりながら頑張れました。
テスト前夜の現実逃避は怖いものです。
残る疑問は何故成功確率が1/13(もしくはそれに近い値)になるのか、ということですが、次のテスト前あたりで考えていきたいと思います。
これを読んでくださった方で、答えにたどり着いた方がいらっしゃいましたら、是非教えて下さい。