蟻コロニー最適化でモンティ・ホール問題を解く

本記事について

勉強のため、私は様々な課題を「人工知能」の手法で解いてきました。

#	課題	手法	記事
1	三角関数	ニューラルネットワーク	機械学習とは何か（ついでにsin(θ)を人工知能で計算）
2	FizzBuzz問題	再帰型ニューラルネットワーク	KerasのRNN(LSTM)でFizzBuzz問題を解いてみた
3	素数判定	ニューラルネットワーク	素数を求めて（ディープラーニング編）
4	モンティ・ホール問題	強化学習	モンティ・ホール問題を機械学習（強化学習）で解く
5	モンティ・ホール問題	deep Q-network	モンティ・ホール問題をゲーム化してDQNにプレーさせてみた
6	じゃんけんグリコ	deep Q-network	じゃんけんグリコでDQNに挑む
7	秘書問題	deep Q-network	DQNに秘書を面接させたら、美しい結果が得られた話
8	ダブルアップ	deep Q-network	強化学習の報酬に関する一研究（ダブルアップを題材に）
9	Hello, world	遺伝的アルゴリズム	ムダにエボリューショナルな"hello, world"の出力方法(Python)
10	モンティ・ホール問題	蟻コロニー最適化	本記事

#10に示すとおり、本記事では、モンティ・ホール問題を、蟻コロニー最適化の手法で解きます。

本記事の流れは、以下のとおりです。

• モンティ・ホール問題についての説明
• 蟻コロニー最適化についての説明
• プログラムについての説明
• プログラムの実行結果
• まとめ

参考文献は以下です。

• 小高知宏(著)『機械学習と深層学習《Pythonによるシミュレーション》』オーム社

それでは参りましょう。

モンティ・ホール問題についての説明

本記事では、Wikipedia『モンティ・ホール問題』を少しアレンジします。

モンティ・ホール問題

あなたは愛車の鍵を盗まれました。あなたの前にモンティと名乗る屈強な男が現れ、3つの箱を出してみせました。モンティが言うには、鍵はどれか1つの箱に入っているが、箱を1つ開けるのに1万円必要とのことです。力づくで箱を奪うのは無理そうですし、この場から逃げることもできそうにありません。あなたは最悪3万円の出費を覚悟し、まずは1万円を払い、1つの箱に手を伸ばしますが、そのときモンティが「出血大サービスだ。残りの箱の内、はずれの箱を1つ開けてやろう」と言い、箱を開けてくれました。その箱は、確かにはずれの箱でした。

——さて、以上の状況下で、あなたがとるべき行動は、A)手を伸ばしたその箱を開ける、B)選択を変え、残っているもう1つの箱を開ける、のどちらでしょうか。ただし、行動基準は、出費の期待値が小さいほうを選ぶものとします。また、心の広いモンティは、Bの行動を許すものとします。

問題の正解

正解はB。選択を変え、残っているもう1つの箱を開けてください。Aと比べ、出費の期待値が小さくて済みます。ホントの話。

説明

選択を変えて当たりを引くのは以下1のケース。選択を変えずに当たりを引くのは以下2のケース。

1. 最初に選んだ箱がはずれの場合、モンティがもう一つのはずれを開けてくれるため、残った箱が当たり
2. 最初に選んだ箱が当たりの場合、モンティの行動に関わらず、残った箱ははずれ

1の確率は2/3、2の確率は1/3。ここから、出費の期待値を計算すると、箱の選択を変更する場合は1.33..万円、箱の選択を変更しない場合は1.66..万円となります。

さて、蟻はどんな答えを出してくれるでしょうか。

蟻コロニー最適化についての説明

蟻コロニー最適化

蟻コロニー最適化とは、「グラフを使ってよい経路を探すことで単純化できるような計算問題」の確率的解法です。蟻がコロニー（＝群れ）から餌までの経路を見つける際の挙動からヒントを得たものであり、巡回セールスマン問題に良好な成績を示すことで知られています。

モンティ・ホール問題も、以下のようにグラフ化できるので、蟻コロニー最適化で解くことができるでしょう。

経路A→B→Cと、経路A→C→Bの勝負となります。

ポイント

本記事では、以下をポイントに、蟻コロニー最適化を実装しました。

1. 餌を持たない蟻は、餌を得るまで巣に帰らない。分岐点では、フェロモンの匂いの比に応じ、どちらに行くかを確率的に決定する。
2. 餌を持つ蟻は、来た道を、フェロモンを放出しつつ帰る。巣に着いたら餌を保管する。
3. 蟻の持つフェロモンは有限で、巣までの距離を計算し、ちょうど使い切るように調整しつつ、経路上に均等に放出される。巣に着いたらフェロモンは充電される。
4. 経路上のフェロモンは時間と共に蒸発する。

以上が、アルゴリズムの要点です。詳細は、プログラムを見ていきましょう。

プログラムについての説明

主要クラスは、経路、モンティ、蟻、蟻の管理者の4つです。以下に、プログラムのソースを示します。ソースの全体は、GitHubにあります。

ライブラリのインポート

以下、インポートします。

from enum import IntEnum, auto
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from operator import add

定数クラス

主要クラスの説明に入る前に、2つの定数クラスについて説明しておきます。各クラスが表すものは、コメントに示すとおりです。

class Position(IntEnum):
    '''蟻の位置を表す定数クラス'''
    change = 0
    stay = auto()
    home = auto()


class Select(IntEnum):
    '''蟻の選択を表す定数クラス'''
    change = 0
    stay = auto()

配列のインデクスとして使うことを想定しているため、Intと比較可能なIntEnumとし、かつ、0から開始するよう指定しています。

経路

位置と位置を結ぶ経路について知っており、また、経路上のフェロモンについても管理するクラスです。

class Path():
    '''経路'''
    def __init__(self, evaporation_rate):
        # 蒸発率
        self.evaporation_rate = evaporation_rate

        # フェロモン2サイクル分
        self.lp = len(Position)
        self.p_pheromone = np.full((self.lp, self.lp), 0.001)
        self.pheromone = np.full((self.lp, self.lp), 0.001)

フェロモンを保持する器は、from × toの、2次元の配列とします。numpy.full()に形状と初期値を与えると、全要素を引数に従い初期化した配列を返してくれます。フェロモンを微小な値で初期化しているのは、のちの計算の都合上です。フェロモンを2サイクル分保持するのは、何匹もの蟻が同時に行動している様子を表現するためです。蟻1匹が行動するたびに経路上のフェロモンを更新していたら、蟻がお行儀よく順番に行動していることになってしまいますので。

    def get_dests(self, position):
        '''目的地の候補と、1サイクル前のフェロモンを返す'''
        if position == Position.home:
            dests = [Position.change, Position.stay]
            pheromones = [self.p_pheromone[Position.home][Position.change],
                          self.p_pheromone[Position.home][Position.stay]]

        elif position == Position.change:
            dests = [Position.stay]
            pheromones = [self.p_pheromone[Position.change][Position.stay]]

        elif position == Position.stay:
            dests = [Position.change]
            pheromones = [self.p_pheromone[Position.stay][Position.change]]

        return dests, pheromones

    def append_pheromone(self, from_, to, pheromone):
        '''フェロモンを上塗りする'''
        self.pheromone[from_][to] += pheromone

以上のメソッドは、いずれも、蟻から呼ばれるものです。前者は蟻の往路で呼ばれ、その行動基準となります。後者は蟻の復路で呼ばれ、蟻の放出したフェロモンをバッファ的に蓄積します。

    def update(self):
        '''フェロモンを蒸発、上塗りする'''
        self.p_pheromone *= self.evaporation_rate
        self.p_pheromone += self.pheromone
        self.pheromone = np.full((self.lp, self.lp), 0.001)

updateメソッドは、蟻の管理者から呼ばれます。1サイクル完了するたびに呼ばれ、バッファ的に蓄積したフェロモンを反映します。この更新式は、蟻コロニー最適化の特徴と言えるでしょう。代替として、例えば以下の更新式でもよさそうですが、そうではなく、蒸発と上塗りで更新します。

    def update(self):
        '''フェロモンを更新する'''
        self.p_pheromone += 0.1 * (self.pheromone - self.p_pheromone)
        self.pheromone = np.full((self.lp, self.lp), 0.001)

モンティ

諸悪の根源であるモンティです。モンティ・ホール問題について知っており、箱や鍵、選択について管理します。やりとりの相手は蟻のみです。

class Monty:
    '''モンティ'''
    def __init__(self):
        self.reset()

    def reset(self):
        '''状況を初期化する'''
        self.boxes = np.zeros(3)
        self.choices = np.zeros(3)
        self.key = np.zeros(3)
        self.key[np.random.randint(3)] = 1.

まずは、初期化関連。箱は3つあるので、箱、鍵、選択のすべてについて、要素数が3の配列で保持するのが自然でしょう。いずれも0で初期化し、鍵はランダムに選んだ1つの要素に1を立てます。

    def open(self):
        '''箱を開ける'''
        # 最初に選択する箱はランダム
        self.choices[np.random.randint(3)] = 1.

        # はずれの箱を一つ開ける
        boxes = self.key + self.choices
        choices = np.where(boxes == 0.)[0]
        self.boxes[np.random.choice(choices)] = 1.

        # stayとchangeのインデックスを保持
        self.stay = self.choices.argmax()
        self.change = (self.boxes + self.choices).argmin()

openメソッドは、「あなた」が箱に手を伸ばすところから、モンティがはずれの箱を開けるところまでを実装します。最初に選択する箱は、どれを選んでも等価と考えられるので、ランダムに選択としています。はずれの箱を一つ開ける処理は、numpyの力を借りて、3行のコードで実装しています。配列の要素同士の足し算ができたり、検索条件に一致する要素のインデクスを取得できたりするので、実装が楽です。最後に、箱の選択を変える場合と変えない場合の、対応するインデクスをそれぞれ保持しておきます。numpy.ndarray.argmax()は値が最大の要素のインデクスを、numpy.ndarray.argmin()は値が最小の要素のインデクスを返してくれます。

    def judge(self, select):
        '''蟻の選択を判定'''
        # 蟻の選択を反映
        if select == Select.change:
            self.choices = np.zeros(3)
            self.choices[self.change] = 1.

        elif select == Select.stay:
            self.choices = np.zeros(3)
            self.choices[self.stay] = 1.

        # 正解判定
        return self.choices.argmax() == self.key.argmax()

judgeメソッドは、蟻の位置＝選択を受け、正解かどうかを判定します。選択と鍵のそれぞれについて、0でなく1が立っている要素が1つあり、そのインデクスが一致していれば正解です。numpy.where()でも可能ですが、よりシンプルにnumpy.ndarray.argmax()でインデクスを取得しています。

蟻

蟻です。蟻として行動します。往路では、フェロモンに導かれ経路を歩き、餌を探します。復路では、フェロモンを放出しつつ、巣に帰ります。帰巣本能を備えており、来た道を帰ることができます。

class Ant():
    '''蟻'''
    def __init__(self, path, action_rate):
        self.path = path
        self.action_rate = action_rate
        self.monty = Monty()
        self.pheromone = 1.
        self.position = Position.home
        self.hasFood = False
        self.memory = []

蟻は管理者によって生成されます。その際、経路と行動確率（後述）を教えられます。蟻はモンティを生成します。——なんだか奇妙な表現ですが。そして初期値として、フェロモンは充電満タンの1、位置は巣、餌は持っておらず、まだ歩いていないので歩いた経路についての「記憶」は白紙の状態です。

    def action(self):
        '''行動を起こす'''
        if np.random.random() < self.action_rate:
            self.walk()
        else:
            self.sleep()

    def sleep(self):
        '''何もしない'''
        pass

管理者から指示を受けた蟻は、所定の行動確率に従い、歩きます。これは、一般の蟻コロニー最適化には無い概念です。この概念を加えたのは、働きアリの法則を取り入れようという遊び心——ではなく、サイクルの初めのほうで、蟻のクロックが揃うのを嫌った（行動を適当にバラけさせたかった）ためです。蟻のクロックが揃うと、経路上のフェロモンをグラフに描いたとき、ギザギザした形状になってしまいます。より滑らかで美しく、傾向がわかりやすいグラフにしたかったので、行動確率という概念を加えてみました。

    def walk(self):
        '''歩く'''
        # 往路
        if not self.hasFood:
            if self.position == Position.home:
                # モンティに箱を開けさせる
                self.monty.reset()
                self.monty.open()

            # 目的地の候補とフェロモンを取得し、目的地（選択）を決定
            dests, pheromones = self.path.get_dests(self.position)
            sp = sum(pheromones)
            probability = [p / sp for p in pheromones]
            self.select = np.random.choice(dests, p=probability)

            # 選択した位置に移動
            self.memory.append(self.position)
            self.position = self.select

            # 餌を得られたか判定
            self.hasFood = self.monty.judge(self.select)

walkメソッドの前半部分、往路です。巣にいる場合のみ、モンティに初期化とはずれの箱を開けるよう指示します。

蟻は経路クラスから、目的地の候補、および、行動基準となるフェロモンを得ます。目的地を選択するアルゴリズムとしては、イプシロン・グリーディ法も考えられますが、今回は「フェロモンの匂いの比に応じ、どちらに行くかを確率的に決定する」とのことですので、ルーレット選択としました。実装はnumpy.random.choice()を利用。ここで重要なのは、フェロモンに多少の差があっても、選択される余地を残すことです。そうしないと、初期の偶然性に左右されるプログラムとなってしまいます。

目的地を選択できたら、移動し、餌を得られたかどうかを、モンティに判定させます。なお、移動の際には、来た道を帰ることができるように、リストself.memoryをスタックとして使います。往路ではappend()しながら歩き、逆に、帰路ではpop()しながら帰ることになります。

        # 復路
        else:
            if self.position == Position.home:
                # 餌を巣に保管し、フェロモン充電
                self.hasFood = False
                self.pheromone = 1.

            else:
                # フェロモン放出量を計算
                lm = len(self.memory)
                pheromone = self.pheromone / lm

                # フェロモンを放出しつつ、帰る
                dest = self.memory.pop()
                self.path.append_pheromone(dest, self.position, pheromone)
                self.pheromone -= pheromone
                self.position = dest

walkメソッドの後半部分、復路です。位置が巣であれば、餌を保管し、フェロモンを充電します。

位置が巣でない場合は、巣に帰るわけですが、来た道を、フェロモンを放出しながら帰るのがポイントです。来た道については、リストself.memoryに記憶しているので、pop()しながら帰ればOKです。フェロモンについては、巣に帰るまで持たせる必要があるので、self.memoryの長さ、すなわち、記憶している経路の数で割ることにより、放出量を決定します。経路を餌場から巣までの一気通貫で見たとき、長い経路ほど単位距離あたりのフェロモンが薄くなるというのが、蟻コロニー最適化の特徴です。この特徴により、長い経路はやがて選択されないようになります。

蟻の管理者

蟻の管理者は、全体を統括します。

class AntController():
    '''管理者'''
    def start(self, cycles, ants_num, action_rate, evaporation_rate):
        '''シミュレーションを開始する'''
        path = Path(evaporation_rate)
        ants = [Ant(path, action_rate) for _ in range(ants_num)]
        history = [[] for _ in range(4)]

        for _ in range(cycles):
            for ant in ants:
                ant.action()
            path.update()
            history[0].append(path.p_pheromone[Position.home][Position.change])
            history[1].append(path.p_pheromone[Position.change][Position.stay])
            history[2].append(path.p_pheromone[Position.home][Position.stay])
            history[3].append(path.p_pheromone[Position.stay][Position.change])

        title = (f'ants_num = {ants_num}, '
                 f'action_rate = {action_rate:.1f}, '
                 f'evaporation_rate = {evaporation_rate:.1f}')

        print('1/2')
        visualize1(title, history)
        print('2/2')
        visualize2(title, history)

経路、蟻を生成し、さらに、経路上のフェロモンの履歴を保持するリストを生成します。その後、サイクルの数だけ、蟻への行動指示と、経路への更新指示、および、履歴の追記を繰り返します。全サイクルが終了したら、履歴を元にグラフを描きます。

プログラムの実行結果

さて、蟻10,000匹を200回転させた結果、経路上のフェロモンは以下グラフのように推移しました。

1. p_home2change：箱の選択を変更する経路のフェロモン
2. p_change2stay：1で開けた箱がはずれだったので、残った1つの箱を開ける経路のフェロモン
3. p_home2stay：箱の選択を変更しない経路のフェロモン
4. p_stay2change：3で開けた箱がはずれだったので、残った1つの箱を開ける経路のフェロモン

まず、12サイクルあたりまで見ましょう。フェロモンは、当たった場合1、外れた場合0.5放出されるため、経路のフェロモンは当たる確率に依存します。箱の選択を変更したほうが当たる確率が2倍高くなるので、p_home2change > p_home2stay となります。箱の選択を変更しないほうが、はずれる（その先に進む）確率が高いので、p_stay2change > p_change2stayとなります。

サイクルを重ねるごとに「箱の選択を変更しない」のは、餌を得られる確率が低い経路ということで、徐々に選択されなくなり、p_home2stayが減少。連れ立ってp_stay2changeも減少します。200サイクルあたりでは、p_home2stay、p_stay2change、共にほぼ0となりました。

以下の図で言うところの、経路A→C→Bは、フェロモンがほぼ0となったわけです。すなわち、この蟻達は、経路A→B→Cが「最短」経路だと学習しました。

以下に、巣から餌場までのフェロモンを、一気通貫で合計したものの推移を示します。

1. p_change：箱の選択を変更する経路A→B→Cのフェロモン（p_home2change + p_change2stay）
2. p_stay：箱の選択を変更しない経路A→C→Bのフェロモン（p_home2stay + p_stay2change）

これが蟻達の出した答えです。

先述のとおり、出費の期待値は、箱の選択を変更する場合は1.33..万円、箱の選択を変更しない場合は1.66..万円です。確かに、箱の選択を変更する経路A→B→Cは、出費が少なくて済む「最短」経路であることが分かります。

蟻達は、見事に正解を導き出しました。

まとめ

勉強のため、私は様々な課題を「人工知能」の手法で解いてきました。

#	課題	手法	記事
1	三角関数	ニューラルネットワーク	機械学習とは何か（ついでにsin(θ)を人工知能で計算）
2	FizzBuzz問題	再帰型ニューラルネットワーク	KerasのRNN(LSTM)でFizzBuzz問題を解いてみた
3	素数判定	ニューラルネットワーク	素数を求めて（ディープラーニング編）
4	モンティ・ホール問題	強化学習	モンティ・ホール問題を機械学習（強化学習）で解く
5	モンティ・ホール問題	deep Q-network	モンティ・ホール問題をゲーム化してDQNにプレーさせてみた
6	じゃんけんグリコ	deep Q-network	じゃんけんグリコでDQNに挑む
7	秘書問題	deep Q-network	DQNに秘書を面接させたら、美しい結果が得られた話
8	ダブルアップ	deep Q-network	強化学習の報酬に関する一研究（ダブルアップを題材に）
9	Hello, world	遺伝的アルゴリズム	ムダにエボリューショナルな"hello, world"の出力方法(Python)
10	モンティ・ホール問題	蟻コロニー最適化	本記事

今回の試み（#10）は、同じ課題を取り上げた#4、5と比較すると、蟻という具体的なイメージを伴うので#4よりも面白く、また、ゲームを実装しなくて済むので#5よりは易しかったです。

また、#5〜8（DQNはどんな攻略法を見つけてくれる？）や、#9（ランダムな文字列は「hello, world」に進化してくれる？）と同じワクワクを感じることができました。

ただし、#3（データの濃度を自然数から実数に変更することで近似に成功）や、#8（DQNに与える報酬を対数に変更することで対極的な攻略法を導出）のような「発見」があったかというと、ありませんでした。それでも、蟻達が美しいグラフを描き出してくれたので、私は満足です。

さて、モンティ・ホール問題や秘書問題のような課題に対し、理論どおりの答えを出してくれる「人工知能」の手法はすごいなと思います。ビジネスの現場で理論がわかっていない課題に対して、統計的、経験的な角度から答えを出してくれるかもしれない、そんな可能性を感じます。

その一方で、素数という概念を学ばせるのに苦労する始末。#3では、ひらめきにより、データの濃度を自然数から実数に変更し、過学習させた結果、100以下の素数を素数判定させることに成功しましたが、当然、外挿結果は正視に耐えないものでした。いや、あえて正視しましょう。

Microsoft社のEric Nichols氏は、「次の素数を発見するようにモデルを訓練することは可能か」という問いに対し「現時点では、一般の数nに対する素数判定さえできない」と答えています。

Nichols氏は、数論的概念を学ぶためには、今現在「ディープ」と呼んでいるネットワークよりも、はるかによりディープなネットワークが必要ではないかと考えています。

これは、裏を返して言えば、既に多種多様な実用化が始まっているにも関わらず、まだ伸びしろがあるとポジティブに捉えることもできるわけで、人類史上でも指折りのフロンティアじゃないかと感じたりもします。

これからも、少しコンセプトは変えつつも、蟻のようにコツコツと、機械学習やディープラーニングに取り組んでいきたいと思います。