【無限の猿定理】ルービックキューブを猿に解かせてみた(Part1: 普通の猿)

はじめに

無限の猿定理をご存知だろうか。wikipediaによると、以下の通りである。

無限の猿定理（むげんのさるていり、英語: infinite monkey theorem）とは、十分長い時間をかけてランダムに文字列を作り続ければ、どんな文字列もほとんど確実にできあがるという定理である。
(中略)
ここではそのような厳密な議論には立ち入らず、この定理の言わんとすることを初等的に考察してみる。話を簡単にするため、タイプライターのキーがちょうど100個あるとする (この数は実際のキー配列での数に近い)。例えば「monkey」という1単語からなる文章は全部で6文字あるので、ランダムにキーを叩いて、「monkey」とタイプされる確率は、

{\displaystyle {\begin{aligned}{\frac {1}{100}}\times {\frac {1}{100}}\times {\frac {1}{100}}\times {\frac {1}{100}}\times {\frac {1}{100}}\times {\frac {1}{100}}&=\left({\frac {1}{100}}\right)^{6}&={\frac {1}{10000000000000000}}\end{aligned}}} (1兆分の1)

である (ここでは、タイプの一様性と独立性を仮定している)。これは非常に小さな確率であるが、0ではないため、猿が「ランダムに6個のキーを打つ」という操作を非常に多くの回数繰り返せば「monkey」という文字列がタイプされる確率は100%に非常に近くなる。
(wikipedia 無限の猿定理より引用)

では、この定理をルービックキューブに置き換えてみよう。ルービックキューブをただひたすら適当に回し続けて偶然揃う確率は0ではないので、十分に長い時間回し続ければ、ほとんど確実に揃えられる。
本記事では、猿にルービックキューブを解いてもらうが、やはり普通の3x3x3のルービックキューブでは難しすぎるので、2x2x2のルービックキューブを解いてもらう。

では、本当にルービックキューブを適当に回して揃えることは可能なのか... 観ていきましょう。

実装

本記事は、

を参考に作成されている。この記事を[1]とする。一部コードも引用しているため、引用箇所を方法の項で明記する。
Google Colabで作成した本記事のコードは、こちらにあります。

各種インポート

各種インポート

import random
import numpy as np
import matplotlib.pyplot as plt

実行時の各種バージョン

Name	Version
numpy	1.21.6
matplotlib	3.2.2

方法

ルービックキューブの状態を表すクラス

ルービックキューブの状態を表すクラスを用意する。[1]の記事から、 apply_moveまでのcp, coの操作箇所を引用している。状態の説明は[1]の記事を参照されたい。

• apply_scramble: 同じ面を連続して回さないように適当に回すメソッドで、同じ面を連続で回さないようにした条件下で回す面、回す角度90°, 180°, -90°が等確率で発生する乱数により崩す。
• is_solved: 完成状態を判定するメソッドである。完成状態は、cp = [0, 1, 2, 3, 4, 5, 6, 7], co = [0, 0, 0, 0, 0, 0, 0, 0]で定義される。

ルービックキューブの状態を表すクラス

class State:
    """
    ルービックキューブの状態を表すクラス
    """

    def __init__(self, cp, co):
        self.cp = cp
        self.co = co

    def apply_move(self, move):
        """
        操作を適用し、新しい状態を返す
        """
        new_cp = [self.cp[p] for p in move.cp]
        new_co = [(self.co[p] + move.co[i]) % 3 for i, p in enumerate(move.cp)]
        return State(new_cp, new_co)

    def apply_scramble(self, moves, moves_num):
        """
        適当に崩した状態を返す
        moves: 9種類の1手操作
        moves_num: 適当に回す回数
        """
        faces = sorted(list(moves.keys()))
        list_noF = [3, 4, 5, 6, 7, 8]
        list_noL = [6, 7, 8, 0, 1, 2]
        list_noU = [0, 1, 2, 3, 4, 5]
        choice_list = [list_noF, list_noL, list_noU]
        choice_next = random.randrange(3)
        for _ in range(moves_num):
            rand = random.choice(choice_list[choice_next])
            move_state = moves[faces[rand]]
            self = self.apply_move(move_state)
            choice_next = rand // 3
        return State(self.cp, self.co)
    
    def is_solved(self):
        """
        完成状態か判定
        """
        return (self.co == [0, 0, 0, 0, 0, 0, 0, 0] and self.cp == [0, 1, 2, 3, 4, 5, 6, 7])

ルービックキューブの1手操作を定義

ルービックキューブの1手操作を定義する。[1]の記事から、1手操作U, L, Fのcp, coの状態を引用している。
普通の3x3x3ルービックキューブでは各面に90°, 180°, -90°の3種類の手があるので、全部で（正六面体）6×（回す角度）3=18通りの1手操作がある。一方、本記事の2x2x2のルービックキューブでは中心のキューブがないため、軸となる1つのキューブを動かさないような1手操作のみで、全ての手を表すことができる。また、そうすることで完成状態の状態を一意に表現できる。例えば、上面を90°回すのは、下面を-90°回すことで表現できるので、3x3x3ルービックキューブの半分の手である9種類のみ定義すると良い。

ルービックキューブの1手操作を定義

moves = {
    'U': State([3, 0, 1, 2, 4, 5, 6, 7],
               [0, 0, 0, 0, 0, 0, 0, 0]),
    'L': State([4, 1, 2, 0, 7, 5, 6, 3],
               [2, 0, 0, 1, 1, 0, 0, 2]),
    'F': State([0, 1, 3, 7, 4, 5, 2, 6],
               [0, 0, 1, 2, 0, 0, 2, 1]
               )}
faces = list(moves.keys())
for face_name in faces:
    moves[face_name + '2'] = moves[face_name].apply_move(moves[face_name])
    moves[face_name + '\''] = moves[face_name].apply_move(moves[face_name]).apply_move(moves[face_name])
print(moves.keys())  # -> dict_keys(['U', 'L', 'F', 'U2', "U'", 'L2', "L'", 'F2', "F'"])

適当に回して解く

手順を順を追って説明する。

1. 完成状態を用意する。完成状態は、cp = [0, 1, 2, 3, 4, 5, 6, 7], co = [0, 0, 0, 0, 0, 0, 0, 0]である。
2. 完成状態から100回適当に回して崩す。solved_state.apply_scramble(moves=moves, moves_num=100)は、完成状態(solved_state)から、同じ面を連続で回さないようにした条件下で回す面、回す角度90°, 180°, -90°が等確率で発生する乱数により100回崩す。
3. 2と同様に、同じ面を連続で回さないようにした条件下で回す面、回す角度90°, 180°, -90°が等確率で発生する乱数により適当に回す。イメージは以下のようである。
   1. まず、初期値として0~3のいずれかを出す乱数を用意して、choice_nextに代入する。（例として、choice_next = 1のとき、choice_list[choice_next]により、list_noLが選ばれる。list_noLはL面以外の手のindexが入っている。）
   2. rand = random.choice(choice_list[choice_next])で、list_noLから乱数により1つの数字を選び、randに代入する。（例として、rand = 2が出たとする。）
   3. move_state = moves[faces[rand]]で、その名前の1手をmove_stateに代入する。（例では、rand = 2なので、move_state = 'F2'である。）
   4. scrambled_state = scrambled_state.apply_move(move_state)で、今の状態をその名前の手で回す。
   5. choice_next = rand // 3でrandの3で割った商をchoice_nextに代入する。（例では、rand = 2で、move_state = 'F2'だったので、choice_next = 0である。そして、3.1の手順に戻るが、次はF面を連続して回してはならないので、choice_list[choice_next]により、list_noFが選ばれる。こうして、同じ面を連続して回さないようなアルゴリズムが構成されている。）
      
4. もし、適当に100000000回回しても揃わなかった場合は、for文のelse文によりf'It could not be solved in {num_move} moves.'が出力される。

適当に回して解く

# 完成状態を用意
solved_state = State(
    [0, 1, 2, 3, 4, 5, 6, 7],
    [0, 0, 0, 0, 0, 0, 0, 0]
)

# 完成状態から100回適当に回して崩す
scrambled_state = solved_state.apply_scramble(moves=moves, moves_num=100)

# 同じ面を続けて回さないようにするためのリスト
faces = sorted(list(moves.keys()))
list_noF = [3, 4, 5, 6, 7, 8]
list_noL = [6, 7, 8, 0, 1, 2]
list_noU = [0, 1, 2, 3, 4, 5]
choice_list = [list_noF, list_noL, list_noU]

choice_next = random.randrange(3)
for num_move in range(100000000):
    if scrambled_state.is_solved():
        print(f'It was solved in {num_move} moves.')
        break
    rand = random.choice(choice_list[choice_next])
    move_state = moves[faces[rand]]
    scrambled_state = scrambled_state.apply_move(move_state)
    choice_next = rand // 3
else:
    print(f'It could not be solved in {num_move} moves.')

出力結果（乱数の種を指定してないので毎回変わる）

It was solved in 719531 moves.

1万個の異なるルービックキューブを適当に回して解くコード

適当に回して解くの項のコードを1万回試行する。つまり、1万個の異なるルービックキューブを適当に回して解く。下記のコードは1.5日くらいと結構時間がかかる。事前に1万回試行した結果は、nums_solve_move.npyより、ダウンロードできる。以降、こちらのファイルを使用する。

1万個のルービックキューブを適当に回して解く

import random
import numpy as np


def main():
    # 9種類の1手操作を全部定義する
    moves = {
        'U': State([3, 0, 1, 2, 4, 5, 6, 7],
                [0, 0, 0, 0, 0, 0, 0, 0]),
        'L': State([4, 1, 2, 0, 7, 5, 6, 3],
                [2, 0, 0, 1, 1, 0, 0, 2]),
        'F': State([0, 1, 3, 7, 4, 5, 2, 6],
                [0, 0, 1, 2, 0, 0, 2, 1]
                )}
    faces = list(moves.keys())
    for face_name in faces:
        moves[face_name + '2'] = moves[face_name].apply_move(moves[face_name])
        moves[face_name + '\''] = moves[face_name].apply_move(moves[face_name]).apply_move(moves[face_name])

    random.seed(2022)  # 再現性のため乱数の種を指定

    # 完成状態を用意
    solved_state = State(
        [0, 1, 2, 3, 4, 5, 6, 7],
        [0, 0, 0, 0, 0, 0, 0, 0]
    )

    trials = 10000  # 試行回数
    nums_solve_move= np.zeros(trials)  # 保存先の配列

    # 同じ面を連続で回さないようにするためのリスト
    faces = sorted(list(moves.keys()))
    list_noF = [3, 4, 5, 6, 7, 8]
    list_noL = [6, 7, 8, 0, 1, 2]
    list_noU = [0, 1, 2, 3, 4, 5]
    choice_list = [list_noF, list_noL, list_noU]

    for trial in range(trials):
        print(f'[trial num {trial}]')
        # 完成状態から100回適当に回して崩す
        scrambled_state = solved_state.apply_scramble(moves=moves, moves_num=100)

        choice_next = random.randrange(3)
        for num_move in range(100000000):
            if scrambled_state.is_solved():
                print(f'It was solved in {num_move} moves.')
                nums_solve_move[trial] = num_move
                break
            rand = random.choice(choice_list[choice_next])
            move_state = moves[faces[rand]]
            scrambled_state = scrambled_state.apply_move(move_state)
            choice_next = rand // 3
        else:
            print(f'It could not be solved in {num_move} moves.')
            nums_solve_move[trial] = np.nan
    np.save('nums_solve_move.npy', nums_solve_move)


class State:
    """
    ルービックキューブの状態を表すクラス
    """

    def __init__(self, cp, co):
        self.cp = cp
        self.co = co

    def apply_move(self, move):
        """
        操作を適用し、新しい状態を返す
        """
        new_cp = [self.cp[p] for p in move.cp]
        new_co = [(self.co[p] + move.co[i]) % 3 for i, p in enumerate(move.cp)]
        return State(new_cp, new_co)
    
    def apply_scramble(self, moves, moves_num):
        """
        適当に崩した状態を返す
        moves: 9種類の1手操作
        moves_num: 適当に回す回数
        """
        faces = sorted(list(moves.keys()))
        list_noF = [3, 4, 5, 6, 7, 8]
        list_noL = [6, 7, 8, 0, 1, 2]
        list_noU = [0, 1, 2, 3, 4, 5]
        choice_list = [list_noF, list_noL, list_noU]
        choice_next = random.randrange(3)
        for _ in range(moves_num):
            rand = random.choice(choice_list[choice_next])
            move_state = moves[faces[rand]]
            self = self.apply_move(move_state)
            choice_next = rand // 3
        return State(self.cp, self.co)
    
    def is_solved(self):
        """
        完成状態か判定
        """
        return (self.co == [0, 0, 0, 0, 0, 0, 0, 0] and self.cp == [0, 1, 2, 3, 4, 5, 6, 7])

if __name__ == '__main__':
    main()

結果

事前に1万回試行した結果は、nums_solve_move.npyより、ダウンロードできる。

1万個のルービックキューブを適当に回して解いた結果

nums_solve_move = np.load('nums_solve_move.npy')

# 統計情報
print('試行回数:', len(nums_solve_move))
print('平均値:', round(np.mean(nums_solve_move)))
print('中央値:', int(np.median(nums_solve_move)))
print('不偏標準偏差:', round(np.std(nums_solve_move,ddof=1)))
print('最小値:', int(np.min(nums_solve_move)))
print('最大値:', int(np.max(nums_solve_move)))

# ヒストグラムに表示
fig, axs = plt.subplots(2, 2, figsize=(14, 14))
axs[0, 0].hist(nums_solve_move)
axs[0, 0].set_xlabel('moves')
axs[0, 0].set_ylabel('counts')
axs[0, 1].hist(nums_solve_move, density=True)
axs[0, 1].set_xlabel('moves')
axs[0, 1].set_ylabel('probability density')
axs[1, 0].hist(nums_solve_move, bins=1000, cumulative=True)
axs[1, 0].set_xlabel('moves')
axs[1, 0].set_ylabel('counts')
axs[1, 1].hist(nums_solve_move, bins=1000, density=True, cumulative=True)
axs[1, 1].set_xlabel('moves')
axs[1, 1].set_ylabel('cumulative probability')

plt.show()

出力結果

試行回数: 10000
平均値: 3608165
中央値: 2507114
不偏標準偏差: 3591931
最小値: 15
最大値: 32444318

平均が約360万回となりました。最小値はなんと15回でとても幸運な猿ですね。最大値は約3200万回と平均から9倍もかかってしまいました。中央値は、平均値よりも小さく約250万回となりました。
ちなみに、360万回で解けた場合、どれくらい時間がかかるでしょうか。1秒に1回回す猿に1日1時間ルービックキューブで遊ばせた場合、

\frac{3600000\ 回}{3600\ 回/日} = 1000\ 日

となる。2×2×2のルービックキューブは、1000日（約3年）あれば解けるので結構現実的に解けそうですね。ただ、10000万匹の猿に1つずつルービックキューブを与えると、最大値を引いた猿は、約3200万回かかるとして約24年ですね（笑）。2世代で解くかもしれないですね。

考察

ルービックキューブを適当に解くときの結果のヒストグラムの分布は数学的に何に従うのでしょうか。ここで、ルービックキューブの取り得る組み合わせは、wikipediaによると、以下の通りである。

コーナーキューブの取り得る位置の順列が8!通り、コーナーキューブが独立に取り得る向きは37通り、空間に対するキューブの向きは24通りある。
よって考えられる組み合わせは：

{\frac  {8!\times 3^{7}}{24}}=7!\times 3^{6}=3674160

(wikipedia ポケットキュープより引用)

ここで、以下の3つの方針が成り立っている仮定して、結果のヒストグラムの分布を考察する。

1. 全てのルービックキューブの状態は、同様に確からしい。
2. 本記事で定義した適当に回すは、実際に適当に回せている。
3. 本記事の崩した状態は本当に崩した状態である。

1については、例えばある状態から完成状態に揃えるのと、ある状態から任意の別の状態に揃えるのは同じくらいの難しさで解けることを仮定する。人にとって完成状態は特別であるが、数学的には特別な意味を持たないと考えられる。
2については、同じ面を連続で回さないようにした条件下で回す面、回す角度90°, 180°, -90°が等確率で回すことを適当に回すと呼んでいるが、実際に適当に回せていると仮定する。少し引っかかるところとして、回す角度90°, 180°, -90°が等確率で良いのか、180°と90°を同じ価値として扱って良いのかということであるが問題ないと仮定する。
3については、完成状態から100回適当に回して崩した状態としているが、実際に適当に崩した状態となっていると仮定する。もちろん、よりたくさん回した方が崩せると考えられるが、処理速度とトレードオフになる。

2, 3の仮定より、例えば360万回目で完成状態になった場合、
0回目は、完成状態から適当に100回回した状態 崩した状態
1回目は、完成状態から適当に101回回した状態 崩した状態
2回目は、完成状態から適当に102回回した状態 崩した状態
・・・
3599999回は、完成状態から適当に3600099回回した状態 崩した状態
3600000回は、完成状態から適当に3600100回回した状態 崩した状態 = 完成状態!
のようなイメージである。つまり、100回以上であれば全て同じように崩した状態に変わりない。完成状態から適当に100回回した状態と完成状態から適当に3600099回回した状態は崩した状態という意味でその状態は同じ価値である。
そして、1の仮定より、崩した状態は全ての状態のいずれかになるので、崩した状態が完成状態になる確率は、状態の組み合わせが3674160通りより、1/3674160である。つまり、2×2×2のルービックキューブを適当に解くまでに回した回数は、「当たりが1本ある3674160本のくじ引きで当たりが出るまで引いた回数は？ただし、引いたくじは元に戻す。」という問題設定と同値である。

当たる確率が1/3674160のくじ引きの当たるまでの分布

当たる確率$p$が$p=1/3674160$のくじ引きを引いたくじを元に戻すとき当たるまでに引く回数$moves$は、$moves$回外れを引いて$moves$回目で当たりを引くので、

P(moves)=p(1-p)^{moves}

となる。この式をグラフで表す。

1 / 3674160で当たりが出るくじの確率分布

moves = 30000000
p = 1 / 3674160
probability = [p*(1-p)**move for move in range(moves)]
cumsum_probability = np.cumsum(probability)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 6))
# 確率分布を表示
ax1.plot(np.arange(moves)[::100], probability[::100])
ax1.set_xlabel('moves')
ax1.set_ylabel('probability')
# 累積確率分布を表示
ax2.plot(np.arange(moves)[::100], cumsum_probability[::100])
ax2.set_xlabel('moves')
ax2.set_ylabel('cumulative probability')

plt.show()

出力結果

*[::100]としている理由は、プロット点が多すぎるとグラフに表示するまでの時間がかるためである。3000万点も表示する必要もないので、100点置きに30万点プロットしている。

1 / 3674160で当たりが出るくじの確率分布と結果を重ねて表示

当たる確率が1/3674160のくじ引きの当たるまでの分布と1万個の異なるルービックキューブを適当に回して解いた結果を重ねて表示する。
事前に1万回試行した結果は、nums_solve_move.npyより、ダウンロードできる。

1 / 3674160で当たりが出るくじの確率分布と結果を重ねて表示

moves = 30000000
nums_solve_move = np.load('nums_solve_move.npy')

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
# 確率分布を表示
ax1.hist(nums_solve_move, bins=np.linspace(0, moves, 100), density=True, 
         label="Results of 10,000 Rubik's Cubes")
ax1.plot(np.arange(moves)[::100], probability[::100], color='red', label='1 / 3674160 distribution')
ax1.set_xlabel('moves')
ax1.set_ylabel('probability density')
ax1.legend()
# 累積確率分布を表示
ax2.hist(nums_solve_move, bins=np.linspace(0, moves, 1000), density=True, cumulative=True, 
          label="Results of 10,000 Rubik's Cubes")
ax2.plot(np.arange(moves)[::100], cumsum_probability[::100], color='red', label='1 / 3674160 distribution')
ax2.set_xlabel('moves')
ax2.set_ylabel('cumulative probability')
ax2.legend()

plt.show()

出力結果

確かに見た目が一致していることがわかる。

おまけ(3×3×3のルービックキューブを適当に回して解くまでに回した回数の分布の予想)

3×3×3のルービックキューブの取り得る組み合わせは、wikipediaによると、以下の通りである。

ルービックキューブをいったん分解して組み立てなおしたときに考えられる色の配置の総数を求めると、まずコーナーキューブの位置が$8!$通り、向きが$38$通り、エッジキューブの位置が$12!$通り、向きが$2^{12}$通り、これらを全てかけあわせて$(8!×38)×(12!×2^{12})$通りとなる。しかし、実際には完全に揃った状態のキューブに回転操作を施すだけではこれだけの組み合わせは実現できない。

1. コーナーキューブとエッジキューブの順列の偶奇は一致する
2. 全てのエッジキューブの位置が揃っている場合、向きが異なっているエッジキューブの個数は偶数個である
3. 全てのコーナーキューブの位置が揃っている場合、時計回りに向きがずれているコーナーキューブの個数と反時計回りに向きがずれているコーナーキューブの個数は3を法として合同である
   以上の3つの条件から、完全に揃った状態のキューブに回転操作を施してできる組み合わせの総数は前述の値を$(2×2×3)$で割ったものとなる。逆に、上記の3つの条件を満たしていれば6面がそろった初期状態に戻すことができる。すなわち、このパズルで考えられる配置は

\frac{(8!×38)×(12!×212)}{2×2×3}=4325京2003兆2744億8985万6000通り

である。
(wikipedia ルービックキューブより引用)

考察の項から、3×3×3のルービックキューブを適当に回して解くまでに回した回数の分布は、「当たりが1本ある43252003274489856000本のくじ引きで当たりが出るまで引いた回数は？ただし、引いたくじは元に戻す。」という問題設定と同値であると考えられる。
もし、4325京回で解けた場合、どれくらい時間がかかるでしょうか。1秒に1回回す猿に24時間フル稼働体制で回してもらった場合は、

\frac{4325 \times 10^{16}\ 回}{86400\ 回/日} \sim 5.0 \times 10^{14}\ 日 \sim 約1.4兆年

途方に暮れますね。では、世界人口80億人全ての人に3×3×3のルービックキューブを与え、1秒に1回1日8時間適当に回してもらった場合は、

\frac{4325 \times 10^{16}\ 回}{28800\ 回/日 \times 80\times 10^8\ 人 } \sim 190000 \ 日/人 \sim 520\ 年/人

となる。つまり、520年に誰か1人は適当に回して3×3×3のルービックキューブを解くことができる。

まとめ

以上より、2×2×2のルービックキューブを適当に回して解くまでに回した回数は、「当たりが1本ある3674160（2×2×2のルービックキューブの組み合わせの総数）本のくじ引きで当たりが出るまで引いた回数は？ただし、引いたくじは元に戻す。」という問題設定と同値であると考えられる。1秒に1回回す猿に1日1時間2×2×2のルービックキューブで遊ばせた場合、360万回で解けるとすると約3年かかることになる。
この考え方に基づいて、3×3×3のルービックキューブでは約4325京通りの組み合わせがあるので、世界人口80億人全ての人に3×3×3のルービックキューブを与え、1秒に1回1日8時間適当に回してもらった場合は約520年かかることになる。
やはり、3×3×3のルービックキューブを適当に揃えることはとんでもない確率であるが、世界が立ち向かえば5世紀程度で解けてしまうという面白さが、無限の猿定理にはある。ただ、毎日8時間も適当に回していたら、もはや猿になっているかもわかりませんが（笑）。

関連記事

参考資料