梅村長幸、考え中。

蛇の瞳に恋してる

奥様の妊娠

誰も知り合い見てないから書くけど、奥様が妊娠した。まだ三ヶ月ぐらいで安定してない。結婚して6年目、妊娠とか都市伝説ではないかと思い始めた矢先である。 夫婦生活もいい意味でマンネリ化してきて、いい意味で退屈ではあった。不安も多いが 概ねポジティブに考えている。ごめん、というか超嬉しい。 奥様は仕事していて、仕事が好き。まあ、色々大変だろうし、制限も増えてしまうだろうが出来る限りの協力をして楽しんでいけたらと思う。

統計検定2級 結果

合格でした。まじでほっとした。PDFのならびが上から下に並んでいるので落ちたかと思ったよ。。。

異常検知と変化検知 - 第3章 k近傍法

モチベーション

この本読み始めました。専門知識とか、数学の基礎知識が多く私が理解するにはちょっとハードルが高めだったのですが、Pythonで可視化して、どのようなアルゴリズムなのかを理解しながらゆっくり進めたいと思います。このエントリーでは 3章の「近傍法による異常検知」から k近傍法を可視化したうえで方法を説明できればと思っています。 アルゴリズムの詳細をかききるのは難しいかと思いますので、気になった方は本を参照いただければと思います。

k近傍法のサマリ

簡単な例として、異常と正常のラベルがついた (x, y) の二次元のサンプルがあるとします。サンプルそれぞれの点から k 個の近い点をチェックしてその正常値、異常値が含まれる割合を求めて、異常値と判断する割合の閾値を計算します。その後、新しい値が来たときに同じような割合の計算をして設定した閾値を超える場合はAlertを出すようにします。

詳細は本書と、外部ですが以下のスライドがわかりやすいかと思います。 異常検知と変化検知 第4章 近傍法による異常検知

今回は正規分布をランダムにだす numpy の関数 randn() を使ってサンプル値を作成します。

num = 100 # 正常値の数
inum = 5 # 異常値の数
# randn()で正常値を 100個、異常値を中心をずらして5個 設定
smpl = [[randn(), randn(), 0] for i in range(num)]
smpl += [[randn()*0.4 + 2.0, randn()*0.4 + 2.0, 1] for i in range(inum)]

図示すると以下のようになります。 f:id:bython-chogo:20170624200213p:plain

参照数 k と 閾値 F値 をサンプルより計算する

少し驚いたのですが、サンプルから最適な参照する近傍数kと閾値であるF値が高くなるa_thの値を変更しながら一番最適な閾値を見つけるのが事前にやる作業になります。これが経験分布ってことで正しいのでしょうか。そのあたりわかっていません。

kの候補の値を1~6まで、a_thの候補の値を0~5まで0.1まで細かく確認して、F値を計算します。

        klist = [1, 2, 3, 4, 5, 6]
        athlist = [0.1*(1+i) for i in range(50)]
        # 現在までの最大値を記録 [k値、a_thの値、最終的な閾値F値]
        maxset = [0, 0, 0.0]
        for kv in klist:
            for av in athlist:
                self.k = kv
                self.ath = av
                tmp_set = [kv,av,self.calc_ath()]
                #print "k, ath, f", tmp_set
                if tmp_set[2] >= maxset[2]:
                    maxset = list(tmp_set)

閾値を比較するためにのF値の計算は calc_ath()の関数で計算します。F値に関しては一章に書かれてますので、ご参照いただければと思います。

k 値と a_thの値から導き出された F値をプロットしてみます。 f:id:bython-chogo:20170624190238p:plain ちょっとわかりずらいからと思いますが、k=4の時のF値が 0.3636…と一番高い値になります。このF値が基準として、異常値の判定をおこないます。

F値が超える範囲にいろをつけてみた

参照する近傍数K=4で設定して、0.3636を超える場所を薄赤で表示してみます。

    def plot_abn(self):
        k, ath = self.check_better_comb()
        print k, ath
        xl = []
        yl = []
        
        # 縦横-3.0 から 3.0まで0.1でF値を計算してプロット
        for x in range(60):
            xd = x*0.1 - 3.0
            for y in range(60):
                yd = y*0.1 - 3.0
                dltList = []
                # 各地点とサンプル間の距離を計算
                for j, smp in enumerate(self.smpl):
                    dltList.append([self.delta(smp[0:2], [xd, yd]), j])

                nml, abnml = 0, 0
                n_cor, ab_cor = 0, 0
                # 近い順に k 個のサンプルを抽出してF値を算出。
                # 閾値を超えた値をxl, ylに記録
                for j, dl in enumerate( sorted(dltList)):
                    if j == k:
                        break

                    if self.smpl[dl[1]][2] == 1:
                        abnml += 1
                    else:
                        nml += 1

                if np.log(self.pi0*abnml)-np.log(self.pi1*nml) > ath:
                    xl.append(xd)
                    yl.append(yd)

        plt.xlim(-3.0, 3.0)
        plt.ylim(-3.0, 3.0)
        plt.scatter(xl, yl, color='r', alpha=0.1)
        self.ploting()
        plt.show()

表示すると以下のようになります。異常値である赤の周りが異常値として判定されるいるのがわかります。多少正常値のサンプルも入ってきますが、アラートを上げる範囲としては妥当といえるかと思います。 f:id:bython-chogo:20170624190618p:plain

ソースコード

以下、今回利用したソースコードになります。ご確認を

import numpy as np
import matplotlib.pyplot as plt
import random
from numpy.random import *
%matplotlib inline

num = 100
inum = 5
smpl = [[randn(), randn(), 0] for i in range(num)]
smpl += [[randn()*0.4 + 2.0, randn()*0.4 + 2.0, 1] for i in range(inum)]

class KNearest():
    def __init__(self):

        self.num = num
        self.inum = inum

        self.pi0 = self.num*1.0/(self.num+self.inum)
        self.pi1 = self.inum*1.0/(self.num+self.inum)
        
        self.k = 5
        self.ath = 0.1
        self.cand = []
        
        self.smpl = smpl

    # サンプルのPlot
    def ploting(self):
        plt.xlim(-3.0, 3.0)
        plt.ylim(-3.0, 3.0)
        for s in self.smpl:

            if s[2] == 0:
                plt.plot(s[0], s[1], 'bo')
            else:
                plt.plot(s[0], s[1], 'ro')

    # 2地点の距離の計算
    def delta(self, a, b):
        return (a[0]-b[0])**2 + (a[1]-b[1])**2

    # a_thからF値計算。。。データ消してしまった orz..後ほど
    def calc_ath(self):

    # 最適なF値の計算
    def check_better_comb(self):
        klist = [1, 2, 3, 4, 5, 6]
        athlist = [0.1*(1+i) for i in range(50)]
        
        maxset = [0, 0, 0.0]
        for kv in klist:
            for av in athlist:
                self.k = kv
                self.ath = av
                tmp_set = [kv,av,self.calc_ath()]
                #print "k, ath, f", tmp_set
                if tmp_set[2] >= maxset[2]:
                    maxset = list(tmp_set)

        return maxset[0], maxset[2]

    # F値のプロット
    def plot_better_comb(self):
        klist = [1, 2, 3, 4, 5, 6]
        athlist = [0.1*(1+i) for i in range(50)]
        
        result = {}
        maxset = [0, 0, 0.0]
        for kv in klist:
            result[kv] = []
            for av in athlist:
                self.k = kv
                self.ath = av
                result[kv].append(self.calc_ath())

        for kv in klist:
            plt.plot(athlist, result[kv], label='k =' + str(kv))
        plt.legend()
    
    # F値を超越して異常と判定される範囲をプロット
    def plot_abn(self):
        k, ath = self.check_better_comb()
        print k, ath
        xl = []
        yl = []
        
        
        for x in range(60):
            
            xd = x*0.1 - 3.0
            for y in range(60):
                yd = y*0.1 - 3.0
                dltList = []
                for j, smp in enumerate(self.smpl):
                    dltList.append([self.delta(smp[0:2], [xd, yd]), j])

                nml, abnml = 0, 0
                n_cor, ab_cor = 0, 0
                for j, dl in enumerate( sorted(dltList)):
                    if j == k:
                        break

                    if self.smpl[dl[1]][2] == 1:
                        abnml += 1
                    else:
                        nml += 1

                if np.log(self.pi0*abnml)-np.log(self.pi1*nml) > ath:
                    #print xd, yd
                    xl.append(xd)
                    yl.append(yd)

        plt.xlim(-3.0, 3.0)
        plt.ylim(-3.0, 3.0)
        plt.scatter(xl, yl, color='r', alpha=0.1)
        self.ploting()
        plt.show()

2017年6月統計検定2級 反省会会場

日曜に受けた統計検定2級の公式解答がでてきました。
http://www.toukei-kentei.jp/wp-content/uploads/ans2017j_grade2.pdf

bython-chogo.hatenablog.com

結論から言うと、24/35 、率としては68%となりました。うーん微妙。よくよく見てみると基本的なところが間違っていて全然理解できないことがよくわかりました。。というわけで反省会!赤のところが間違っているところで、コメントをつけたしておきました。ただ見直してもよくわかってないようで、いろいろまずい部分が浮き彫りになっています。


[1] ④
Ⅰ ○ TV保有率はずば抜けていてPC保有率と被らない故
Ⅱ ○ PC保有率の第3四分位数12都道府県はDVD/BD保有率と被らない
Ⅲ × SP保有率とMP保有率はオーバーラップするため高いかどうかの判断がつかない
[2] ③ MPの最高値は75~80でb, SPの最高値は80~85でa,
[3] ③ 茨城の中央値は28841で足していくとD内
[4] ② 同じく長野の第1、第3四分位数の値まで足していくと、BとD
[5] ⑤ 北海道のA~Fまで足しても20%に満たない(エ)、秋田のA~Bを足して10%に満たない(ウ)
[6] ③ 1668/995 = 1.68, 68%増
[7] ②
Ⅰ × 季節性の傾向ならば前年同月比でもよくね?
Ⅱ ○ 図より拡大している
Ⅲ × 伸び率は訪問数の差分に影響しない
[8] ③ 自信なし!自己相関係数とコレログラムをお勉強
=> やはりまちがっていた!でも正解の絞り方がわからない。。。①なのかー①なんだなー。

[9] ④ 目分量、強い負の相関だけど、-0.9まではないか
[10] ② 目分量、ちょうど間とおてる
<span style="color: #ff0000">[11] ③ F値が、推定値/標準誤差で -6.27, 自由度はサンプル数25-1</span>
=> いやー、どうどうと間違った解答をw 自由度は 25-2 みたいです。。。おいうことでたぶん⑤

[12] ④
① × 制度は悪くなる
② × 迷った。各層内の散らばりを小さくする必要があるのか?
③ × 渋谷でやったら渋谷ピーポーの傾向が強くなる。無作為にはならない
④ ○ 非がない
⑤ × 渋谷ピーポーのお友達紹介したら渋谷ピーポーの傾向が強くなる。
[13] ⑤ 自信なし。④だったか。
=> ③でした。wで恥の上塗り。
  b = (X-Y)/2 - (ε1 - ε2)/2 だから、σ**2/2, なのかなぁ。。

[14] ④ 超パニックになっていたが、0.1 * 0.4 だわね
[15] ③ 規格外率 0.1*0.4 + 0.05*0.3 + 0.02*0.3
[16] ② 条件付き確率 (0.1*0.4) / [15]
[17] ③ 5C4 (0.7)**4 (0.3)**1
=> ああ、間違えてもうた。そうか、4戦目までに名人が3勝1敗の場合、4C3 (0.7)**3 (0.3)**1 で5戦目に名人が一勝するために、0.7をさらにかけるのね。すると、0.28812で①でした!
[18] ② 多分間違えた⑤が正解か。6戦目までで引き分けている確率 6C3(0.7)**3 (0.3)**3 = 0.185<
=> やはり⑤でした。稼げた問題だけに超もったいない2問です。
[19] ④ 自信ない。複数の標本が合わさるやつ苦手ね。
[20] ④
Ⅰ ○ 足しても引いても平均0は0
Ⅱ ○ [19]が正しいこと前提だと真だが、果たして
Ⅲ × わからない。。。
=> Ⅲは○だったようです!これはどういうことだ。
[21] ③ 自信ない。正規分布で、0.05の上側は1.64となるので、カイ2乗だから、1.64の2乗つーことで選んでみました。
=> 完全に導き方ミスりましたね。カイ2乗分布表から0.5 自由値1の値をみつければ。3.84ですので④
[22] ⓶ [21]から、Wn/2 = 1.355, カイ2乗分布でそれを超えるのが、自由度7のとき、じゃあ nは1足して8! それらしいあてずっぽだし、[21]の解答に依存する
=> ラッキーパンチ!一定確率であてずっぽが当たるのは皆平等に与えられた権利です。これはやっぱり分布表からでしょうね。2nが0.05を超える直前の値は7ですね。するとnは8になるようです。
[23] できなかった。とりあえず③と埋めてみた

=> やっぱりミスってた。んでまだ答えがわからんぞ。がんばろ。
[24] ⓸ 標準誤差 = 標準偏差/root(n)、正直パニクった
[25] ⑤ 計算して出てきたのがこれだけど、ミスってないか不安。1.64で算出
[26] ①
Ⅰ ○ 不偏・・・だよな?
Ⅱ × 許容範囲を広めるので、広くなる
Ⅲ × 1/root(3)倍になるはず。。。
[27] ② 自由度はそれぞれ1引く
[28] ① 27*49/100
[29] ② (2-1)*(2-1)
[30] ④ あれ、足せばいいんですよね・・・?

=> あちゃー、不偏分散は(n-1)で割らなきゃダメです!割ったら ①ですね。
[31] ⑤ 地域は4つなので、自由度は3、F値は平均平方の比率
[32] ② Prが0.0405で5%以下、なので棄却できない

=> これはやばい。検定の基本的な部分を間違えて覚えちまっている。5%以下なので、帰無仮説を棄却!しなきゃだめ。①

[33] ② 問15は理解するのに時間がかかりパニクった。時間なくなるし。31.9 + 0.3 * 28 - 4.43
[34] ① Prが全部5%未満なので採用

=> 有意に「正」となるパラメータね。。。ひっかかった。β2は負のパラメータだわさ。
[35] ⑤ 
① × 博士号と理論研究の関係性はわからない
② × PhDAgeが少ないが最盛期も少なくなる
③ × 理論研究者のほうが4年半ほど早く最盛期に到達する
④ × t値の判定から、4年半の差は統計的に有意である
⑤ ○ だから、理論研究やれば4年半早く最盛期迎えられるの!

 

 なんどもいうようですが、やばい。合格できるのかはなぞですが、基本的な部分がことごとく間違った理解をしていることが明らかになりました。[8][13][20][23]についてはきちんと解答が導けないありさま。受かるかまだわかりませんが、まだまだ精進が必要!

統計検定試験その後の散歩

前回のエントリの通り、本日統計検定2級をうけてきました。実は午後に準1級も申し込んでいたのですが、ほとんど勉強できなかったうえ、2級でぎりぎりの状態だったのでやめることにしました。 会場は文京区の中央大学キャンパス。そこの近くに文京シビックセンターによってから帰ることにしました。文京区の区役所なんですが実はここの上が無料の展望台なんです。東京23区のど真ん中にあるため、23区のランドマークがたくさん見渡せるのでおすすめ。 f:id:bython-chogo:20170618231923j:plain

統計検定会場となった中央大学のキャンパスです。いい眺めでしょ。

シビックセンターをあとにしたあと、ふらふらと歩いていたら神社が。そこには金毘羅宮とかかれていました。そういえばブラタモリでやっていたことを思い出してお参りに。

f:id:bython-chogo:20170618232129j:plain

そこから変なスイッチが入ってしまいました。よっしゃこのまま歩いて神田明神に行ってみようじゃないの!って感じで、水道橋駅から神田川をのぼって御茶ノ水神田明神に到着。

f:id:bython-chogo:20170618232251j:plain

なんか巫女さんとか宮司さんとかがいそいそと歩いていて、行事をやっている様子。境内の真ん中に藁のワッカができていて説明書きを読むとその中を通りながら8の字歩くとご利益があるようだけど、さすがに一人では恥ずかしいので却下。お参りだけして後にしました。

もうここまで来たら湯島もいっちゃうでしょ。一度はいったスイッチぶっ壊れちゃってるようです。

f:id:bython-chogo:20170618232602j:plain

境内に入るとなにやら儀式をしている音が。そこでは神前の結婚式がとりおこなわれているようでした。音に合わせて2人の巫女さんが舞を踊っているようで、思わず前前前世を唄ってしまいそうなのをこらえながらお参りを済ませました。

いや、ここまでいろんな神様にお願いしたから統計検定も大丈夫でしょう!とルンルンで途中で見つけた定食屋にはいって、遅めの昼飯でラーメン・チャーハンセットを幸せにいただきました。きっと神様にも私の祈りが届いたのでしょうか、もしくは普段の行いかよかったのか、店をでたら激しい雨が降り始めてきてしまい、駅までダッシュしたのですがびしょ濡れになるはめに。。。きっと神様が僕の浅はかな信仰心とあるく不謹慎ぶりを見抜いて罰をあたえたもうたのでしょう。

統計検定の結果も思いやられます。

2017年6月統計検定2級 自己チェック

2017年統計検定2級を受けてきました。とりあえずわからなかったところも含めて自分の解答をさらします。決して模範解答ではありません、ていうか、間違いだらけだと思いますので、つっこみ大歓迎です。落ちてたら超恥ずかしいですが、公式回答がくるまで落ち着かないのでやってみました。

[1] ④
Ⅰ ○ TV保有率はずば抜けていてPC保有率と被らない故

Ⅱ ○ PC保有率の第3四分位数12都道府県はDVD/BD保有率と被らない

Ⅲ × SP保有率とMP保有率はオーバーラップするため高いかどうかの判断がつかない
[2] ③ MPの最高値は75~80でb, SPの最高値は80~85でa,
[3] ③ 茨城の中央値は28841で足していくとD内
[4] ② 同じく長野の第1、第3四分位数の値まで足していくと、BとD
[5] ⑤ 北海道のA~Fまで足しても20%に満たない(エ)、秋田のA~Bを足して10%に満たない(ウ)
[6] ③ 1668/995 = 1.68, 68%増
[7] ②
Ⅰ × 季節性の傾向ならば前年同月比でもよくね?
Ⅱ ○ 図より拡大している
Ⅲ × 伸び率は訪問数の差分に影響しない
[8] ③ 自信なし!自己相関係数とコレログラムをお勉強
[9] ④ 目分量、強い負の相関だけど、-0.9まではないか
[10] ② 目分量、ちょうど間とおてる
[11] ③ F値が、推定値/標準誤差で -6.27, 自由度はサンプル数25-1
[12] ④
① × 制度は悪くなる
② × 迷った。各層内の散らばりを小さくする必要があるのか?
③ × 渋谷でやったら渋谷ピーポーの傾向が強くなる。無作為にはならない
④ ○ 非がない
⑤ × 渋谷ピーポーのお友達紹介したら渋谷ピーポーの傾向が強くなる。
[13] ⑤ 自信なし。④だったか。
[14] ④ 超パニックになっていたが、0.1 * 0.4 だわね
[15] ③ 規格外率 0.1*0.4 + 0.05*0.3 + 0.02*0.3
[16] ② 条件付き確率 (0.1*0.4) / [15]
[17] ③ 5C4 (0.7)**4 (0.3)**1
[18] ② 多分間違えた⑤が正解か。6戦目までで引き分けている確率 6C3(0.7)**3 (0.3)**3 = 0.185
[19] ④ 自信ない。複数の標本が合わさるやつ苦手ね。 
[20] ④
Ⅰ ○ 足しても引いても平均0は0
Ⅱ ○ [19]が正しいこと前提だと真だが、果たして
Ⅲ × わからない。。。
[21] ③ 自信ない。正規分布で、0.05の上側は1.64となるので、カイ2乗だから、1.64の2乗つーことで選んでみました。
[22] ⓶ [21]から、Wn/2 = 1.355, カイ2乗分布でそれを超えるのが、自由度7のとき、じゃあ nは1足して8! それらしいあてずっぽだし、[21]の解答に依存する
[23] できなかった。とりあえず③と埋めてみた
[24] ⓸ 標準誤差 = 標準偏差/root(n)、正直パニクった
[25] ⑤ 計算して出てきたのがこれだけど、ミスってないか不安。1.64で算出
[26] ①
Ⅰ ○ 不偏・・・だよな?
Ⅱ × 許容範囲を広めるので、広くなる
Ⅲ × 1/root(3)倍になるはず。。。
[27] ② 自由度はそれぞれ1引く
[28] ① 27*49/100
[29] ② (2-1)*(2-1)
[30] ④ あれ、足せばいいんですよね・・・?
[31] ⑤ 地域は4つなので、自由度は3、F値は平均平方の比率
[32] ② Prが0.0405で5%以下、なので棄却できない
[33] ② 問15は理解するのに時間がかかりパニクった。時間なくなるし。31.9 + 0.3 * 28 - 4.43
[34] ① Prが全部5%未満なので採用
[35] ⑤ 
① × 博士号と理論研究の関係性はわからない
② × PhDAgeが少ないが最盛期も少なくなる
③ × 理論研究者のほうが4年半ほど早く最盛期に到達する
④ × t値の判定から、4年半の差は統計的に有意である
⑤ ○ だから、理論研究やれば4年半早く最盛期迎えられるの!

あらためて見直してみると自信ないし基本的なところ抑えられてないし、めっちゃ恥ずかしい。自信ないとこ全部だめだったら落ちてるなぁ、こりゃ。答え出てきたら改めて採点してみます。

ランダム迷路をコーディング

モチベーション

昔よく遊んだゲームで「トルネコの大冒険」というのがありました。いわゆるローグ系の走りで毎回変わるダンジョンの配置に驚かされたものです。今回はそのダンジョンをランダムに作る・・・と思ったのですが、頭で考えた限りではちょっと設定がめんどくさそうかなぁと。そこで、トルネコのダンジョンのかなり下の階にでてくる、ランダム迷路をつくってみようかと思いました。いまの私の実力にちょうどいい。よし、それをPythonでやってみようぞ!

まずは結果から

こんなんを出力しました。にょきにょき伸びていく感じがなんとも気持ち悪くて俺好みですよ。 f:id:bython-chogo:20170608232524g:plain

きちんと迷路になっています。定義としてはループをつくらない、マスを全部埋めていく、ということを念頭においてました。まあ、よくみると埋まっていないますもあるのですが、まあ、おおめに見てください。

全部のコード

import numpy as np
import random
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.patches as patches
%matplotlib inline
from IPython.display import HTML

class MakeMaize():
    def __init__(self):
        self.size = [60,60]
        self.holes = [[[1,1,1,1] for i in range(self.size[1])] for j in range(self.size[0])]
        self.holes[0][0] = [0, 0, 1, 1]
        self.driller = []
        self.died = []
        self.dlog = []

    #上下左右のマスに動けるかチェック。[下、左、上、右]でいけたら1
    def check_available(self, plr):
        direction = [1,1,1,1]
        
        # 今来た道は戻れない
        direction[(plr["f"]-2)%4] = 0
        # 壁にいるときは壁側に移動できない
        if plr["p"][0] == 0: direction[1] =0
        if plr["p"][1] == 0: direction[0] =0
        if plr["p"][0] == self.size[0] - 1: direction[3] = 0
        if plr["p"][1] == self.size[1] - 1: direction[2] = 0

        # それ以外で、すでに埋まっている場所を確認(holesに記録)
        x = plr["p"][0]; y = plr["p"][1];
        for i, d in enumerate([[0,-1],[-1,0],[0,1],[1,0]]):
            dx = x + d[0]; dy = y + d[1];
            if direction[i] ==1:
                #print "DX:DY:", dx, dy, self.holes[dx][dy]
                ix = 1
                for j in self.holes[dx][dy]:
                    ix *= j
                direction[i] = ix
        # 結果を返す
        return direction
    # 3から9までの値をランダムに返す。分岐、折れ曲がりのタイミングを設定
    def ran_gene(self):

        return random.choice(range(3,10))
        
    # ルーティン drillerが場所を定義
    def step_forward(self, c):
        def direc(f, p, g):
            x, y = 0, 0
            if f == 0: y = - 1
            if f == 1: x = - 1
            if f == 2: y = 1
            if f == 3: x = 1
            
            while len(self.dlog) != c+1:
                self.dlog.append([])
                

            self.dlog[c].append([[p[0], p[0]+x], [p[1], p[1]+y]])

            return [p[0]+x, p[1]+y]
        
        # Drillerの1ステップ
        for d, drl in enumerate(self.driller):

            # 現在の動ける場所を確認
            av_dire = self.check_available(self.driller[d])
            # 動ける場所がない場合は動ける場所にワープ
            if av_dire == [0, 0, 0, 0]:
                cands = []
                # 過去ログ中からランダムに動ける場所に移動
                pick = random.choice(self.dlog)
                random.shuffle(pick)
                for p in pick:
                    cand = self.holes[p[0][0]][p[1][0]]
                    if sum(cand) > 0:
                        ordr = [0, 1, 2, 3]
                        random.shuffle(ordr)
                        for o in ordr:
                            if cand[o] == 1:
                                self.driller[d]["f"] = o
                                self.driller[d]["p"] = [p[0][0], p[1][0]]
                                break
                        if self.driller[d]["p"] == [p[0][0], p[1][0]]: break
                        
                    
            # ターン、分岐のカウントを一つずつ減らす
            self.driller[d]["tc"] -= 1
            self.driller[d]["sc"] -= 1
            # 現在の進む向き
            f = self.driller[d]["f"]
            av_dire[(f-2)%4] = 0

            # 現在の向きに進める場合、ターンするカウントが0出ない場合
            # 現在の方向に一つ進む
            if av_dire[f] == 1 and self.driller[d]["tc"] != 0:
                self.driller[d]["p"] = direc(f, self.driller[d]["p"], d)
                av_dire[f] = 0
                
            # 進める方向からランダムに一つ方向をえらぶ
            # 前に進めない場合の対処
            way = []           
            for i, v in enumerate( av_dire):
                if v == 1:
                    way.append(i)

            if len(way) > 0:
                af = random.choice(way)

                # ターンカウントが0のとき
                # 進める方向にターンする
                if self.driller[d]["tc"] <= 0:
                    self.driller[d]["p"] = direc(af, self.driller[d]["p"], d)
                    self.driller[d]["f"] = af

                # 分岐カウントが0のとき
                # 新しいdrillerを作成してリスト追加、進める方向に進ませる
                elif self.driller[d]["sc"] <= 0:
                            
                    #print "\tAF", af, av_dire, 
                    
                    self.driller.append({
                            #"p": direc(f, self.driller[d]["p"]),
                            "p": self.driller[d]["p"],
                            "f": af,
                            "sc": self.ran_gene(),
                            "tc": self.ran_gene()
                        })

                av_dire[af] = 0

            else:
                if d not in self.died:
                    self.died.append(d)

            # 各カウントが0になったときにランダムに値を設定
            if self.driller[d]["tc"] == 0:
                self.driller[d]["tc"] = self.ran_gene()
            if self.driller[d]["sc"] == 0:
                self.driller[d]["sc"] = self.ran_gene()
                    
            self.holes[self.driller[d]["p"][0]][self.driller[d]["p"][1]] = list(av_dire)
        return 0

    # ログの表示
    def print_now(self, c):
        print c
        log = []
        for i, drl in enumerate(self.driller):
            
            print "   ", "PLY", i, "P:", drl["p"], "F:", drl["f"], "SC:", drl["sc"], "TC", drl["tc"]
            log.append(drl["p"])
        self.logs.append(log)
            
    def start(self):
        # 初期のdrillerの配置
        self.driller.append({"p": [0, 0], "f": 3, "tc": self.ran_gene(), "sc": self.ran_gene()})
        #self.driller.append({"p": [self.size[0]-1, self.size[1]-1], "f": 1, "tc": self.ran_gene(), "sc": self.ran_gene()})

        # 各drillerを動かす。
        c = 0
        while(True):
            
            #self.print_now(c)
            code = self.step_forward(c)
            if code < 0:
                break

            if len(self.died) == len(self.driller):
                break
            c += 1
            
            if c == 500:
                break
        
        dlog = []
        p_old = []

        print len(mm.dlog) 

これを使って可視化します。

mm = MakeMaize()
mm.start()

fig = plt.figure()
ax = plt.gca(axisbg='black')
ax.set_xlim(-1, mm.size[0]+1)
ax.set_ylim(-1, mm.size[1]+1)
fig.patch.set_facecolor('black')

def init():
    #line.set_data([0,0], [0,0])
    #return line,
    return []

def animated(i):
    patches = []
    for xy in mm.dlog[i]:
        x, y = xy
        patches.append(ax.add_patch(
                plt.Rectangle((x[0], y[0]), x[1]-x[0], y[1]-y[0], 
                color="white", edgecolor="white", linewidth="2.5")))

    return patches



ani = animation.FuncAnimation(fig, animated, frames=len(mm.dlog), init_func=init, interval=100, blit=True)
ani.save("makemaize.gif",writer='imagemagick')
HTML(ani.to_html5_video())