「続・わかりやすいパターン認識」11章「ノンパラメトリックベイズ」の「クラスタリングの事前確率」について

昨日の「続・わかりやすいパターン認識」読書会にて、「ホップの壺や中華料理店過程のシミュレーションをみると、これを使うと均等にクラスタリングされるのではなく、クラスタサイズが大きいものから順に小さくなっていくようなクラスタリングがされるように見えるのだが、その認識で正しいのか」といった感じの質疑があった。
いい質問。


実は「続・わかりやすいパターン認識」(以降「ぞくパタ」)では、 p225 の「クラスタリングの事前確率の考え方」のところに、ダイレクトにではないがその質問の答えにつながることが書いてあったりする。coffee break というコラムの形になっているので、つい読み飛ばしちゃった人も多いかもしれないが、結構大事なことが書いてあるので一度じっくり読んでみるといい。


そのあたりも含めて読書会でフォローした内容をここにメモしておく。


まずそもそもの話として。
ベイズにおいて、事前確率(事前分布)の選び方に論理的な裏付けがありうるのか、というと、実のところ全く何もない。ぶっちゃけ事前分布はなんでもいい。
手前味噌だが、そのあたりのことは gihyo.jp での機械学習連載で書いたので、興味があれば一度。


理論的にはなんでもいいが、そこに人間様の都合が絡むと、こういう事前分布は困る、とか、こういうのは嬉しい、とか出てくる。
望ましい事前分布の条件として、どんな問題でも共通して求められるとすればこの 2つだろう。

  • ありうる事象(今回の場合は分割方法)の確率が非ゼロ
  • なんらかの方法で事後分布が計算できる

共役事前分布は後者の「計算できる」を最優先にした選び方になる。共役でなくても、今回のディリクレ過程のように「ギブスサンプリングできる」でもいい。そこらへんの詳細は 12章でやる話なので、ここではとりあげない。


ベイズ公式を見ればわかるが、事前確率がゼロである事象は、たとえどんなデータをもってきても、事後確率もゼロである。よって「ありうる事象」、つまり「答えとして考えられるすべての可能性」には非ゼロの確率を与えておかないといけない。
今やりたいことは「クラスタ数をあらかじめ決めない」なので、クラスタ数が3個でも4個でも5個でも、何個であってもちゃんと非ゼロの確率を持つ事前分布が必要なわけである。


さて、これで冒頭の質問に答える準備ができた。
ここで大事なことは、事前確率は非ゼロでさえあればどんなに小さくても(理論的には)いいということ。
ホップの壺(中華料理店過程)は、色(テーブル)ごとの玉(客)数の偏りが大きい分布だが、均等な分割、たとえば 500個の玉について 100, 100, 100, 100, 100 という割り方にもちゃんと非ゼロの確率が割り振られる。現実には到底観測されることなんかあり得ないほど低い確率だろうけど、ゼロでさえなければ、データが何とかしてくれる期待が持てるわけだ。


ただし、いくら何でもいいといっても、事前分布が「変」であるほど、「正しい答え」にたどり着くには膨大なデータ数が必要になる。逆に言えば、データ数が少ないと事前分布に引きずられる。
例えば 3 クラス× 50 データの iris データセットを 12章のディリクレ過程混合モデルに食わせると、100, 50 の2個のクラスタにまでは分かれるのだが、そこから分割が進まなかったりする*1。k-means (もちろん3クラス決め打ち)などなら、50, 50, 50 に近いところまで行くのだが。
これはおそらく事前分布であるディリクレ過程に引きずられているのだろう。


というわけで冒頭の疑問の答えは、均等な分割もちゃんと行われうるから大丈夫! という点では No だが、事前分布に負けて、希望よりも偏ったクラスタになることもあるという点では Yes でもある。
それなら、中華料理店過程のような偏った分割モデルではなく、もっと一様な分割を考えたらいいのでは? という発想が当然出てくるだろうが、実は可算無限な台を持つ一様分布は作れない(ちょっと考えれば雰囲気はわかる、はず)のと同じ理由で、どうしたって傾斜を持たせるしかなかったりする。


つまり、ディリクレ過程(やホップの壺や中華料理店過程)は、上の「任意クラスタ数で非ゼロ&計算できる」という条件を満たす中で一番易しいものなのだ。
とはいえ、他にもノンパラメトリックベイズのための無限分割モデルは解きたい問題にあわせて考えられており、ぞくパタ 11.4 でちらっと紹介されているピットマン・ヨー過程*2はその代表的な一つである。


読書会では、ピットマン・ヨー過程がなぜかついでに紹介されているっぽい刺身のツマのような扱いを受けていて、あわてて言語処理ではめっちゃ重要ですから! と突っ込んだり(笑)。
ある分野でどんなに有名であっても、その隣接分野ですら知られているとは限らないというよくある現象なんだけどね。

*1:ぞくパタに書かれている「全パターンが 1 つのクラスタに所属」という初期化だと、実はそれすらもままならず、150個が 1 クラスタのままから動かなかったりする

*2:ぞくパタで「ピットマン・ヨー過程」だといって紹介されているのは、正確には「中華料理店過程をピットマン・ヨー過程に対応するように拡張したもの」ではあるが

中華料理店過程のテーブル数の分布が見てみたい・実験編 #ぞくパタ

今日のぞくパタ読書会の予習に「続・わかりやすいパターン認識」(以降「ぞくパタ」)の 11章をつらつら読む。
p227 に「CRP における使用テーブル数の変化」というグラフがあり、αが 2 と 10 のそれぞれにおいて、来客数が 1000 になるまでシミュレーションしたときのテーブル数の推移が示されているのだが、当然ながら一回の試行しか描かれていない。
中華料理店過程(CRP)ないしホップの壺 をちょっこり試すだけなら 10行ちょっとのスクリプトでできる話であり、実際やってみると、実行するたびにテーブルの増え方が結構変わる。つまり一回の試行ではぶっちゃけよくわからん。


というわけで、α=2 について 1000 回シミュレーションを行い、来客数がそれぞれ 100, 200, 400, 600, 800, 1000 になったときのテーブル数の分布を箱ひげで描いてみた。



ほとんどは中心の周辺にきれいにまとまるけど、散らばりもそれなりにある。
また、テーブル数がある程度増えた以降は、分散はほとんど増えなくなっていく(散らばり具合は n が少ないうちにだいたい決まっちゃう)ということもわかる。テーブルの増える確率がどんどん小さくなることから想像できる通りではあるが。


ちなみに、それぞれの平均値と標準偏差はこんな感じだった。

100: mean=8.31(8.39) sd=2.47(2.42)
200: mean=9.69(9.77) sd=2.78(2.68)
400: mean=11.03(11.14) sd=3.07(2.93)
600: mean=11.86(11.95) sd=3.20(3.06)
800: mean=12.41(12.53) sd=3.25(3.15)
1000: mean=12.88(12.97) sd=3.27(3.22)

カッコ内はそれぞれの理論値。1000 サンプルにしては悪くない?
平均値(期待値)は前回の記事「ディリクレ過程(中華料理店過程)のトピック数(テーブル数)の期待値を導出してみる #ぞくパタ - Mi manca qualche giovedi`?」で求めていたが、分散も「独立な確率変数の和の分散は、分散の和」*1であることに注意すればまったく同様に求めることができる。


ここでは厳密さが欲しかったわけじゃないのでシミュレーションでお茶を濁したけど、実は二手間くらいかければ確率質量関数も計算できる。もし知りたかったら「ポアソン二項分布」で検索。


大したものではないが、シミュレーションのスクリプトも一応貼っておく。
中華料理店過程は、事後分布とか絡まなければ、実装はとても易しく楽しいので、興味あればこんなん見ないで自分で実装してみる方がいい。

import numpy

alpha = 2.0

labels = [100, 200, 400, 600, 800, 1000]
target = dict()
data = []
for i, x in enumerate(labels):
    target[x] = i
    data.append([])

for i in xrange(1000):
    n = []
    for N in xrange(1001):
        z = N + alpha
        p = [ni / z for ni in n]
        p.append(alpha / z)
        k = numpy.random.choice(len(p), p=p)
        if k<len(n):
            n[k] += 1
        else:
            n.append(1)
        if N in target:
            data[target[N]].append(len(n))

for i, N in enumerate(labels):
    p = alpha / (numpy.arange(N) + alpha)
    m = p.sum()
    v = ((1 - p) * p).sum()
    print "%d: mean=%.2f(%.2f) sd=%.2f(%.2f)" % (N, numpy.mean(data[i]), m, numpy.std(data[i]), numpy.sqrt(v))

import matplotlib.pyplot as plt
plt.boxplot(data, labels=labels)
plt.show()

*1:正確には無相関で十分

ディリクレ過程(中華料理店過程)のトピック数(テーブル数)の期待値を導出してみる #ぞくパタ

「続・わかりやすいパターン認識」(以降「ぞくパタ」)の11章「ノンパラメトリックベイズモデル」を読んでいる。


ぞくパタはこの11章・12章のための本、と言ったらさすがに言いすぎなのかもしれないけど、そのつもりで読んでる人はけっこう多そう。
ノンパラベイズをしたかったら英語の論文読むしかない状態だったのが、日本語で読めるようになったというだけでも十分価値があるわけだが、ディリクレ過程と、無限次元ディリクレ分布や中華料理店過程との関係などもきちんとフォローされているし、ディリクレ過程とはどのようなふるまいをするものなのか、あれこれ実験して示してもくれていて、これは力作だなあと頭が下がる。


といっても、そもそも題材が本質的に易しくないので、これほどていねいに解説されてあっても、さらっと読んだだけではまったく意味はわからない。紙と鉛筆無しで理解できたら、ウン千人かに一人の天才以上だろう。
読書会もまもなくノンパラベイズに突入するが、それなりの覚悟を持って予習するのがおすすめ。


とまあ前振りはこれくらいで。
ぞくパタ 11章の p227 にて

CRP(中華料理店過程)では、客に使用されているテーブル数は、来客数 n の増加とともに増大し、その使用テーブル数の期待値 E[c] は、理論的には E[c] = O(αlog n) となることが知られている。

とあった。
テーブル数(つまりトピック数)の期待値が分かれば、集中度 α を決める目安になるので、どうやって求めるのかちょっと興味を持って探してみたら、英語 wikipedia の Chinese Restaurant Process のページに載っていた。

(ぞくパタとは記号の使い方が異なるので注意)


が、式見た瞬間に導出方法がわかるくらい簡単なものだったので(答え見る前に自力でわかればよかったのだけどw)、せっかくだから紹介してみよう。


n 人の客が着席したときの、テーブル数を表す確率変数 c = c_n を考える。
また確率変数 X_k を「 k 人目の客が来たときに新しいテーブルに着いたら 1、すでに客のいるテーブルに着いたら 0 」(つまり k 人目の客のときのテーブル数の増分)とすると、 c は X_k の和に等しい。

  • c = \sum_{k=1}^n X_k


k 人目の客が新しいテーブルに着く確率は \frac{\alpha}{k-1+\alpha} なので、

  • E[X_k]=1\cdot\frac{\alpha}{k-1+\alpha}+0\cdot\frac{k-1}{k-1+\alpha}=\frac{\alpha}{k-1+\alpha}


ここで、確率変数の和の期待値は、期待値の和になる(独立でなくても成立するが、今回は独立なのでどっちでもいい)ので、

  • E[c] = \sum_{k=1}^n E[X_k] = \sum_{k=1}^n \frac{\alpha}{k-1+\alpha}


とわかる。
さらに digamma 関数Ψの特性 \Psi(z+1)=\Psi(z)+\frac1z を使えば E[c]=\alpha\{\Psi(\alpha+n)-\Psi(\alpha)\} と表せる。実装や見積もりではこちらが速くて便利だろう。
クラスタ数(トピック数)はわからないとはいっても目星くらいはあるだろうから、この期待値がその目安となる数字を少し超えるくらいにαを選ぶと良さそう。


あとは E[c] = O(αlog n) であることを大雑把に示してみる。
といっても \alpha\approx 1 として、さらに和を積分で近似してしまえばあっさり出る*1

  • E[c]=\sum_{k=1}^n \frac{\alpha}{k-1+\alpha}\approx \sum_{k=1}^n \frac{\alpha}{k} \approx \int_1^n \frac{\alpha}{k}dk = \alpha\log n


ちなみに、alpha = 1.0, n = 500 のとき(ぞくパタ 12章の実験の設定)、E[c] = 6.8 なのだが、αlog n = 6.2 とそれなりの近さ。
まあ元の式も十分簡単なので近似の意味はたいしてなさそうだけど、計算機を使わずに頭の中で目安を立てるくらいになら使えるかも。

*1:真面目にやるなら、ちゃんと上からおさえなきゃなんだろうけど、ごちゃごちゃするだけなのでネグる

Python Lasagne でニューラルネットするチュートリアル その 2

昨日の記事の続き。

まずは、Lasagne のハマリポイントを紹介しながらコードの説明。
そのあと、ディープラーニング? ディープニューラルネットワーク? どちらの呼び名が正しいのかビミョウに自信ないが、要は畳み込み& max-pooling &ドロップアウトを織り交ぜたモデルを学習させようとするとハマリポイントが増えるので、そのあたりの注意点とサンプルコード。

def digits_dataset(test_N = 400):
    import sklearn.datasets
    data = sklearn.datasets.load_digits()

    numpy.random.seed(0)
    z = numpy.arange(data.data.shape[0])
    numpy.random.shuffle(z)
    X = data.data[z>=test_N, :]
    y = numpy.array(data.target[z>=test_N], dtype=numpy.int32)
    test_X = data.data[z<test_N, :]
    test_y = numpy.array(data.target[z<test_N], dtype=numpy.int32)
    return X, y, test_X, test_y

データセットづくり。
ありきたりなコードだが、ここにすでに一つハマリポイントが隠されている。

特徴量 X の方は特筆することはなく、レコード数×次元のシンプルな ndarray でいいのだが。
Lasagne の出力層の非線形関数に softmax を指定した場合、出力層に与える正解ラベル y の型は暗黙に int32 が要求される。
うっかりここに普通の int の ndarray とか渡したりすると、

TypeError: ('Bad input argument to theano function with name "****.py:**" at index 1(0-based)', 'TensorType(int32, vector) cannot store a value of dtype int64 without risking loss of precision. If you do not mind this loss, you can: 1) explicitly cast your data to int32, or 2) set "allow_input_downcast=True" when calling "function".', array([...

と怒られて、どこのことを指しているのか見当つかず悩むハメになる。
このエラーが出たときは、正解ラベル y を numpy.array(*, dtype=numpy.int32) でくるむといい。

#### model
n_classes = 10
batch_size=100

l_in = lasagne.layers.InputLayer(
    shape=(batch_size, input_dim),
)
l_hidden1 = lasagne.layers.DenseLayer(
    l_in,
    num_units=512,
    nonlinearity=lasagne.nonlinearities.rectify,
)
l_hidden2 = lasagne.layers.DenseLayer(
    l_hidden1,
    num_units=64,
    nonlinearity=lasagne.nonlinearities.rectify,
)
model = lasagne.layers.DenseLayer(
    l_hidden2,
    num_units=n_classes,
    nonlinearity=lasagne.nonlinearities.softmax,
)

モデルの定義は入力層 lasagne.layers.InputLayer から始めて、lasagne.layers.* の第一引数に前の層を指定しながらつないでいき、最後に適当な出力層を宣言すると、それがそのままモデルの参照となる。


入力層は shape で入力層の次元を指定するのは自然だが、一度に入力するデータのレコード数もこのモデルを定義する段階で指定する必要がある。これがハマリポイントその……いくつ目だっけ?
おそらく実装上の都合だろうが、適当な batch_size を決めて、入力データを batch_size 件ずつに分けて学習を回すのが Lasagne 流になる。ここで shape の値に、入力データの全データ件数を指定すると、学習のところでデータを分割して回す必要がなくなってコードの一部がすっきりするのだが、(おそらく)最急降下法相当になり、足がめちゃめちゃ遅くなる。
というわけで適当な batch_size を指定する必要がある。あまり小さいと速度が落ちるし、学習が落ち着かずに loss や accuracy が跳ねまわる。大きいと足が遅くなるし、余りのデータが無駄になる。
十分大きいデータなら 500 前後がほどよい印象( mnist.py は batch_size=600 になっている)。このサンプルコードでは、データが 2000件にも満たないので batch_size=100 にしている。


出力層は、ニューラルネットワークで何が作りたいかによって変わるだろうが、いちばん一般的な他クラス分類器なら、サンプルにある通り lasagne.layers.DenseLayer を使って、 num_units にクラス数を、 nonlinearity に lasagne.nonlinearities.softmax を指定すればいい。
クラス数は y.max()+1 とかしてもよかったけど、わかりきってるのでリテラル書いちゃった。

#### loss function
objective = lasagne.objectives.Objective(model,
    loss_function=lasagne.objectives.categorical_crossentropy)

X_batch = T.matrix('x')
y_batch = T.ivector('y')
loss_train = objective.get_loss(X_batch, target=y_batch)

定義したニューラルネットワークモデルから目的関数を生成してくれるところが Lasagne の真骨頂である。
ああ、あとおそらくこの Objective をインスタンス化するまでのどこかのタイミングに、パラメータを格納する SharedVariable も用意してくれている。


loss_train は Theano の expression になっているので、theano.function に食わせれば実行コードにコンパイル済みの関数が得られる。Theano すげー。
あとは、この loss_train を使って、必要な関数や処理を作っていくことになる。

#### update function
learning_rate = 0.01
momentum = 0.9
all_params = lasagne.layers.get_all_params(model)
updates = lasagne.updates.nesterov_momentum(
    loss_train, all_params, learning_rate, momentum)

#### training
train = theano.function(
    [X_batch, y_batch], loss_train,
    updates=updates
)

Lasagne における学習は、Theano の updates の仕組みを使ってパラメータを更新する。
updates に渡す更新関数も、適切な lasagne.updates.* を呼べば Lasagne が作ってくれる。
Lasagne には単純な SGD も用意されているが、lasagne.updates.nesterov_momentum は Nesterov Moment を使った SGD になる。これは、「パラメータを前回動かした方向にしばらく動かす & 勾配を取る点も前回動かした方向に少しずらす」というもの。DNN 向けに SGD を改良したもので、収束が速くなる……のかな? mnist.py がこれを使っていたので、ここのサンプルもそれに倣った。


lasagne.layers.get_all_params は、前回の最後にもチラッと出ていたが、パラメータを格納した SharedVariable のリストを返すものである。
更新関数はもちろん対象となるそのパラメータたちを知らなければならないので必要性はわかるのだが、それをユーザが書かなければならないところには納得いかない(苦笑)。


学習率や Nesterov Moment の moment パラメータが固定で与えられているのは違和感があるかもしれない。
これを T.dscalar にして、theano.function の引数で与えるようにすれば可変にもできることは確認した。が、生 SGD より制御が難しく、結果をうまく改善できなかったので、mnist.py と同じように固定で与えている。


mnist.py では、この train 関数はデータを SharedVariable に置いて givens で渡す形で書かれている。そして引数はそのデータの中でその呼出で対象となる範囲を指すインデックスのみを指定する。特に GPU を用いる場合はデータの転送コストが一番高いわけで、そうすると mnist.py と同じ方式のほうが確実に効率いいだろう。
このサンプルでは簡素なコードを優先したことと、Lasagne を試している環境が GPGPU を利用できないものだったので(笑)、データを引数で渡すシンプルな形にした。

#### prediction
loss_eval = objective.get_loss(X_batch, target=y_batch,
                               deterministic=True)
pred = T.argmax(
    lasagne.layers.get_output(model, X_batch, deterministic=True),
    axis=1)
accuracy = T.mean(T.eq(pred, y_batch), dtype=theano.config.floatX)
test = theano.function([X_batch, y_batch], [loss_eval, accuracy])

予測のための関数を定義しているところ。
deterministic というパラメータは、おそらく、Dropout などのノイズ層にスルーで渡されて、True のときにはランダムに捨てるのをやめる(学習時のみドロップアウトする)という制御のためだと思われる。
それ以外には特に疑問はないだろう。

#### inference
numpy.random.seed()
nlist = numpy.arange(N)
for i in xrange(100):
    numpy.random.shuffle(nlist)
    for j in xrange(N / batch_size):
        ns = nlist[batch_size*j:batch_size*(j+1)]
        train_loss = train(X[ns], y[ns])
    loss, acc = test(test_X, test_y)
    print("%d: train_loss=%.4f, test_loss=%.4f, test_accuracy=%.4f" % (i+1, train_loss, loss, acc))

ようやく準備が全て終わって推論である。が、batch_size ずつしか訓練に与えることができないので、ここでも自前でコードをちょいと書く必要がある。といっても、読めばわかるレベルなので大丈夫だろう。
せっかく引数でデータを渡すのだから、渡すデータの順序がランダムになるようにした。batch_size に対して余る分はそのイテレーションでは使われないわけだが、ランダムにしておくことで使われないデータの偏りをなくすことも期待している。
テストデータの評価の方は batch_size に分ける必要がないので、一発呼び出しで済んでいる(このときは……)。


以上、これが一番シンプルな Lasagne の使い方。
でも、Lasagne でもっとディープラーニングっぽい画像処理したいという場合にはもう二手間くらい必要になる。
サンプルで使っているデータセット digits は 8x8 の画像なので、これを入力とした畳み込み& max-pooling &ドロップアウトを織り交ぜたモデルのサンプルコードがこちら。

import numpy
import lasagne
import theano
import theano.tensor as T

# dataset
def digits_dataset(input_width, input_height, test_N = 400):
    import sklearn.datasets
    data = sklearn.datasets.load_digits()
    N = data.data.shape[0]
    X = data.data.reshape(N, 1, input_width, input_height)
    y = numpy.array(data.target, dtype=numpy.int32)

    numpy.random.seed(0)
    z = numpy.arange(data.data.shape[0])
    numpy.random.shuffle(z)
    test_X = X[z<test_N]
    test_y = y[z<test_N]
    X = X[z>=test_N]
    y = y[z>=test_N]
    return X, y, test_X, test_y

n_classes = 10
input_width = input_height = 8
X, y, test_X, test_y = digits_dataset(input_width, input_height)
N = X.shape[0]
test_N = test_X.shape[0]
print(X.shape, test_X.shape)


#### model
batch_size=100

l_in = lasagne.layers.InputLayer(
    shape=(batch_size, 1, input_width, input_height),
)
l_conv1 = lasagne.layers.Conv2DLayer(
    l_in,
    num_filters=8,
    filter_size=(3, 3),
    nonlinearity=lasagne.nonlinearities.rectify,
    W=lasagne.init.GlorotUniform(),
    )
l_pool1 = lasagne.layers.MaxPool2DLayer(l_conv1, pool_size=(2, 2))

l_hidden1 = lasagne.layers.DenseLayer(
    l_pool1,
    num_units=256,
    nonlinearity=lasagne.nonlinearities.rectify,
)
l_hidden1_dropout = lasagne.layers.DropoutLayer(l_hidden1, p=0.2)
l_hidden2 = lasagne.layers.DenseLayer(
    l_hidden1_dropout,
    num_units=64,
    nonlinearity=lasagne.nonlinearities.rectify,
)
model = lasagne.layers.DenseLayer(
    l_hidden2,
    num_units=n_classes,
    nonlinearity=lasagne.nonlinearities.softmax,
)

#### loss function
objective = lasagne.objectives.Objective(model,
    loss_function=lasagne.objectives.categorical_crossentropy)

X_batch = T.tensor4('x')
y_batch = T.ivector('y')
loss_train = objective.get_loss(X_batch, target=y_batch)

#### update function
learning_rate = 0.01
momentum = 0.9
all_params = lasagne.layers.get_all_params(model)
updates = lasagne.updates.nesterov_momentum(
    loss_train, all_params, learning_rate, momentum)

#### training
train = theano.function(
    [X_batch, y_batch], loss_train,
    updates=updates
)

#### prediction
loss_eval = objective.get_loss(X_batch, target=y_batch,
                               deterministic=True)
pred = T.argmax(
    lasagne.layers.get_output(model, X_batch, deterministic=True),
    axis=1)
accuracy = T.mean(T.eq(pred, y_batch), dtype=theano.config.floatX)
test = theano.function([X_batch, y_batch], [loss_eval, accuracy])


#### inference
numpy.random.seed()
nlist = numpy.arange(N)
for i in xrange(100):
    numpy.random.shuffle(nlist)
    for j in xrange(N / batch_size):
        ns = nlist[batch_size*j:batch_size*(j+1)]
        train_loss = train(X[ns], y[ns])
    result = []
    for j in xrange(test_N / batch_size):
        j1, j2 = batch_size*j, batch_size*(j+1)
        result.append(test(test_X[j1:j2], test_y[j1:j2]))
    loss, acc = numpy.mean(result, axis=0)
    print("%d: train_loss=%.4f, test_loss=%.4f, test_accuracy=%.4f" % (i+1, train_loss, loss, acc))

元のサンプルコードとよく似ているが、細かいところが結構違うので、そのあたりを中心に説明して Lasagne チュートリアルを終わろう。

  • 2次元データを入力に使うときは、4次元テンソルで渡す必要がある。 digits_dataset() で X を reshape している行を見てもらえば手っ取り早いが、その形も (データ件数, 1, 横次元, 縦次元) と、2次元目がなぜか 1 でないといけない(理由は調べてない。 Theano の制限?)
  • モデルの定義で Conv2DLayer や MaxPool2DLayer や DropoutLayer で畳み込みや max-pooling やドロップアウトを記述できるが、これはコードを見ればわかると思うので説明略。
  • loss function の定義で、入力データを表す変数 X_batch が、1次元のときは T.matrix だったが、2次元では T.tensor4 とする。
  • 1次元のときは updates 関数だけが batch_size のしばりがあったのだが、2次元では loss function そのものにも batch_size のしばりが及ぶようになる。つまり test も一発呼び出しができなくなるので、こちらでも分割してループして結果を平均、といった処理を行う必要がある。テストは訓練と違ってランダムサンプリングするわけにはいかないので、テストデータのサイズは batch_size の整数倍であることが望ましい。

Python Lasagne でニューラルネットするチュートリアル その 1

@nishio さんに教えてもらったのだが、Lasagne というニューラルネットワークPython ライブラリが Kaggle でけっこうよく使われているらしい。
イタリア語読みすると「ラザーニェ」、Lasagna(ラザニア) の複数形なので、まあ日本人が呼ぶ分には「ラザニア」でいい気がする。


2015年6月現在でバージョンが 0.1.dev と、今手を出すのは人柱感満載。
実際、自分の思ったとおりのモデルを動かすのはなかなかに大変だったので、そのメモを残しておく。


インストールは別に難しいところはない。
ただ Theano 前提なので、Python 2.7 でないと動かないし、Windows で動かすのは茨の道だろう。


また、ドキュメントには "Install from PyPI" とあるくせに、pip ではインストールできない(ワナその1)。
ぐぐると、

Lasagne が PyPI からインストールできないんですけど
git clone でインストールできるよ
そりゃそうだろうけど、ドキュメントには "Install from PyPI" って書いてあるよ?
そんなこと言ってもできないもんはできないんだから、ガタガタぬかさず git から入れとけ

みたいなやりとりが引っかかって、ウケる。
というわけで、おとなしく git clone & python setup.py しよう。


インストール後、git clone した場所に examples というディレクトリがあって、かの MNIST を使ったサンプルコードが置いてある。
GPGPU が叩けない環境でも mnist.py と mnist_conv.py というサンプルは問題なく動くので、まずはそれで遊んでみるといい。
標準のサンプルなのにいきなり

The uniform initializer no longer uses Glorot et al.'s approach to determine the bounds, but defaults to the range (-0.01, 0.01) instead. Please use the new GlorotUniform initializer to get the old behavior. GlorotUniform is now the default for all layers.

みたいなワーニングが出るのだが、多分気にしたら負け。


mnist.py は 512 個ずつのユニットを持つ2段の隠れ層からなる古き良きニューラルネットワークで、環境にもよるだろうが2時間半くらい学習して 98.5% くらいの精度が出る。
mnist_conv.py は 5x5 の畳込みと 2x2 の max-pooling を2回重ねたあと、256 ユニットの隠れ層、そしてドロップアウトという今風のディープなニューラルネット。さすがに重く、それでも 27時間ほどで学習を終えて、99.4% の精度を叩き出す。
これが Python のコードをちょちょっと書くだけで動く(ウソ)んだから、楽しそうでしょう?


mnist.py のコードを見るとモデルを定義するのは簡単そうなので、簡単に使えるのかと思って、mnist.py を改造して自前のデータを僕の考えた最強のモデルに食わせようとしたら、図ったように動かない。
まず mnist.py のコードが無駄に複雑で、汎用化しているつもりなんだろうけど、明示していない仕様があれこれあるようで、謎の型エラーがバンバン出る。


よし、改造はあきらめて一からコードを書こう。ドキュメントにはちゃんと TUTORIAL の文字がある(ワナその2)。
開くと、

Understand the MNIST example
TODO:

良かった、紙のマニュアルだったら壁に叩きつけているところだった。電子化バンザイ。


しかたない、MNIST サンプルコードを理解してやろうじゃないか。
と、勢い込んで読み始めるが、学習や予測のためのコードが 100行以上あって、わずか数行で機械学習できる scikit-learn(ぬるま湯) に慣らされた ゆとり には大層ツライ。


ともあれ、そうして一応理解したつもりで、必要最小限にしぼった Lasagne のスモールサンプルコードがこちら。

import numpy
import lasagne
import theano
import theano.tensor as T

#### dataset
def digits_dataset(test_N = 400):
    import sklearn.datasets
    data = sklearn.datasets.load_digits()

    numpy.random.seed(0)
    z = numpy.arange(data.data.shape[0])
    numpy.random.shuffle(z)
    X = data.data[z>=test_N, :]
    y = numpy.array(data.target[z>=test_N], dtype=numpy.int32)
    test_X = data.data[z<test_N, :]
    test_y = numpy.array(data.target[z<test_N], dtype=numpy.int32)
    return X, y, test_X, test_y

X, y, test_X, test_y = digits_dataset()
N, input_dim = X.shape
n_classes = 10
print(X.shape, test_X.shape)


#### model
batch_size=100

l_in = lasagne.layers.InputLayer(
    shape=(batch_size, input_dim),
)
l_hidden1 = lasagne.layers.DenseLayer(
    l_in,
    num_units=512,
    nonlinearity=lasagne.nonlinearities.rectify,
)
l_hidden2 = lasagne.layers.DenseLayer(
    l_hidden1,
    num_units=64,
    nonlinearity=lasagne.nonlinearities.rectify,
)
model = lasagne.layers.DenseLayer(
    l_hidden2,
    num_units=n_classes,
    nonlinearity=lasagne.nonlinearities.softmax,
)

#### loss function
objective = lasagne.objectives.Objective(model,
    loss_function=lasagne.objectives.categorical_crossentropy)

X_batch = T.matrix('x')
y_batch = T.ivector('y')
loss_train = objective.get_loss(X_batch, target=y_batch)

#### update function
learning_rate = 0.01
momentum = 0.9
all_params = lasagne.layers.get_all_params(model)
updates = lasagne.updates.nesterov_momentum(
    loss_train, all_params, learning_rate, momentum)

#### training
train = theano.function(
    [X_batch, y_batch], loss_train,
    updates=updates
)

#### prediction
loss_eval = objective.get_loss(X_batch, target=y_batch,
                               deterministic=True)
pred = T.argmax(
    lasagne.layers.get_output(model, X_batch, deterministic=True),
    axis=1)
accuracy = T.mean(T.eq(pred, y_batch), dtype=theano.config.floatX)
test = theano.function([X_batch, y_batch], [loss_eval, accuracy])


#### inference
numpy.random.seed()
nlist = numpy.arange(N)
for i in xrange(100):
    numpy.random.shuffle(nlist)
    for j in xrange(N / batch_size):
        ns = nlist[batch_size*j:batch_size*(j+1)]
        train_loss = train(X[ns], y[ns])
    loss, acc = test(test_X, test_y)
    print("%d: train_loss=%.4f, test_loss=%.4f, test_accuracy=%.4f" % (i+1, train_loss, loss, acc))


このコードは何をやっているか。

  • データは scikit-learn の datasets に含まれる digits 。0 から 9 までの数字画像(16階調 8x8 ピクセル)が 1797 件。今回 scikit-learn はこのためだけw*1
  • 400 件をテストデータに、残り 1397 件を訓練データに回している。テストデータを切りの良い数字にしているのは次回への振り
  • モデルは隠れ層2層(1層目 512ユニット、2層目 64ユニット)。100周の学習で 97% くらいの精度になる。


細かい解説は次回に回すが、とりあえず Lasagne の守備範囲は、内部 DSL 的に記述されたモデルから、目的関数を生成するところだけということを念頭に置けば、このコードは特に苦もなく読めると思う。
学習におけるパラメータ更新とか、テストデータの評価とかはほぼ Theano 頼みで、現状はそこのつなぎを利用者が書く必要がある(だから書かないといけないコードが多い)。まあ 0.1.dev なんで。


またこのコードでは学習後のモデルを保存していないが(このサンプルデータの規模なら保存する必要もないだろうし)、まじめにやるなら当然その要望は出てくるだろう。
そのときは lasagne.layers.get_all_params(model) がパラメータを格納した Theano の SharedVariable のリストを返すので、こいつらを何らかの方法で永続化するといい。

続き。

*1:なので、わざわざ def してスコープを分けたところで sklearn を import している

「続・わかりやすいパターン認識」の8章「隠れマルコフモデル」の問題点 2つ #ぞくパタ

【追記】
本記事の内容は公式の正誤表ですでに修正済みです。第1版第4刷以降が出ることがあれば、そちらに反映されていることが期待されます。

【/追記】


昨日は ぞくパタ読書会 にのこのこ行ってきた。主催者、発表者、参加者の皆さん、会場を提供してくださったドワンゴさんに感謝。

続・わかりやすいパターン認識―教師なし学習入門―

続・わかりやすいパターン認識―教師なし学習入門―

「続・わかりやすいパターン認識」(以降「ぞくパタ」)の8章「隠れマルコフモデル」を読んだわけだが、この章には理解のさまたげになりうる問題点が大きく2つあると感じた。

  • 自明ではない条件付き独立性を、言及なく使っている
  • ビタービアルゴリズムで求める ψ_t(j) の定義が明記されていない


前者は読書会でも突っ込んだのだが、残念ながらピンときている人はいない様子だった? 式変形を追えばそこで詰まると思うのだが……いや、全員まったく問題なく式変形できた可能性もまだ残っている。

まあそれはともかく。簡単に問題点を整理・フォローしてみる。

自明ではない条件付き独立性を、言及なく使っている

書籍の中で対象となるのは次の番号の式。


見ればわかるが「隠れマルコフで重要な式のほぼ全て」である。
それもそのはず、そもそも隠れマルコフの仮定はまさにその条件付き独立性を得るため導入されており、その条件付き独立性のおかげでこのモデルは計算できるようになっているのだから。


ぞくパタはグラフィカルモデルをやってなく、グラフィカルモデルなしでモデルに含まれる条件付き独立性を導出するのは面倒なので、ネグったのだろう。そこまでは理解できるのだが、条件付き独立性を使っていることを言いもしないのはさすがにまずい。

「この式変形では、隠れマルコフの仮定から導かれる 〜 という条件付き独立性を用いている」と一言あるだけでも、無条件で成立するものではないと知れるし、ちゃんと式の導出を追っかけたいと思った人は条件付き独立性について調べ考える機会が与えられる。
今の状態では、隠れマルコフの仮定がいかにして強力な枠組みを生み出しているのか、その様子を目の当たりにすることもできない。


もったいない。


確かに何もかも理解するのは難しいし、そんな必要なんてないという意見もあるだろう。
が、捨てる順番的には、確率モデルで一番大切な条件付き独立性は最後ちゃうかなあ、と個人的には強く思う。


というわけで、使われている条件付き独立性を紹介し、それを使って実際に式変形する。
式変形は全部ではなくあえて 1つだけにしておくので、興味があれば残りはぜひ自力で。


まずは条件付き独立性。
つか一応「条件付き独立性」って何ってところからやっとくか。
ぞくパタ 7ページに書いてあるのをそのまま引くと、「 S が与えられた下で事象 X, Y は条件付き独立である」とは次のいずれかが成り立つこと。

  • P(X|Y,S)=P(X|S)
  • P(Y|X,S)=P(Y|S)
  • P(X,Y|S)=P(X|S)P(Y|S)


ぞくパタ では「事象」になっているが「確率変数」に対しても同様に定義する。


そして「S が与えられた下で X, Y は条件付き独立である」ことを X\perp Y | S と書くことにする。
本来なら条件付き独立性の ⊥ の縦棒は2本なのだが、はてダの環境では出ないので1本にさせてもらってる。手抜きですまん。


さて、隠れマルコフモデルでは以下のような条件付き独立性が成立する。記号はぞくパタのままなので、notation は省略。

  • (1) x_1,\cdots,x_t\;\perp\;x_{t+1},\cdots,x_n\;|\;s_t
  • (2) x_1,\cdots,x_t \;\perp\;s_{t+1}\;|\;s_t
  • (3) x_1,\cdots,x_{t-1},s_{t-1} \;\perp\; x_t \;|\;s_t


他にもあるのだが、とりあえず。
導出は……省略(苦笑)。上にも書いた通り、隠れマルコフの仮定から導出できなくないのだが、めんどくさすぎる。
ここでの主眼は「条件付き独立性を示す」ことではなく「条件付き独立性を使っていることを言う」なので、いいのだ(開き直り)。


ちなみにグラフィカルモデルは、そのめんどくさすぎる条件付き独立性の導出を「見ただけでわかる」ようにしてくれる強力なツールである。
興味があれば PRML 8章などをどうぞ。


さて条件付き独立性を使って上の式の 1つをやっつけよう。前向きアルゴリズム更新式 (8.10) あたりでいいか。

  • \alpha_t(j)\;=\;\left\{\sum_{i=1}^c\alpha_{t-1}(i)\;a_{ij}\right\}\;b(\omega_j,x_t)   (8.10)

a とかαとかをもとの確率になおす。

  • P(x_1,x_2,\cdots,x_t,s_t=\omega_j)\\ \;=\;\left\{\sum_{i=1}^c P(x_1,x_2,\cdots,x_{t-1},s_{t-1}=\omega_i) \;P(s_t=\omega_j|s_{t-1}=\omega_i) \right\}\;P(x_t|s_t=\omega_j)

長い。


x_1,x_2,\cdots,x_{t-1} はセットでしか動かさないので、それを {\bf x}_1^{t-1} と書くことにする。
また s_t=\omega_js_{t-1}=\omega_i は常に値が指定されているものとして、それぞれ単に s_t, s_{t-1} と記す。
それだけだと右辺の sum のインデックス i の指すものが行方不明になるので、sum は s_{t-1} についてとる形で表す。

  • P({\bf x}_1^{t-1},x_t,s_t)\;=\;\left\{\sum_{s_{t-1}} P({\bf x}_1^{t-1},s_{t-1}) \;P(s_t|s_{t-1}) \right\}\;P(x_t|s_t)

見通しが良くなった。これを示す。
確率の加法定理から

  • P({\bf x}_1^{t-1},x_t,s_t)\;=\;\sum_{s_{t-1}} P({\bf x}_1^{t-1},x_t,s_t,s_{t-1})

である。s_{t-1} を増やして消しただけ。
右辺の sum の中身を、乗法定理を2回使って変形する。

  • P({\bf x}_1^{t-1},x_t,s_t,s_{t-1})
  • =P(x_t|{\bf x}_1^{t-1},s_t,s_{t-1})\;P({\bf x}_1^{t-1},s_t,s_{t-1})
  • =P(x_t|{\bf x}_1^{t-1},s_t,s_{t-1})\;P(s_t|{\bf x}_1^{t-1},s_{t-1})\;P({\bf x}_1^{t-1},s_{t-1})

\alpha_{t-1}(i)=P({\bf x}_1^{t-1},s_{t-1}) なので、分解の3項目は片付いた。
(ここまでは乗法・加法定理しか使っていないので、隠れマルコフでなくても成立する)


1項目は上の条件付き独立性の (3) x_1,\cdots,x_{t-1},s_{t-1}\;\perp\;x_t\;|\;s_t を使うと、

  • P(x_t|{\bf x}_1^{t-1},s_t,s_{t-1})=P(x_t|,s_t)

となり、これは b(\omega_j,x_t) である。

2項目は同じく条件付き独立性 (2) x_1,\cdots,x_t\;\perp\;s_{t+1}\;|\;s_t から

  • P(s_t|{\bf x}_1^{t-1},s_{t-1})=P(s_t|s_{t-1})

となり、a_{ij} である。これで (8.10) が示された。


ここで一番大事なことは、「前向きアルゴリズムやビタービを導出するには隠れマルコフの仮定がちゃんと必要だった」とわかること。
これらのアルゴリズムは隠れマルコフ以外でも出てくるので、そのときになんで使えるかもこの辺りを理解していれば納得しやすい。

ビタービアルゴリズムで求める ψ_t(j) の定義が明記されていない

ビタービは、ψ_t(j) という値を再帰的に求めることで状態の系列を推定するアルゴリズムだが、その ψ_t(j) について、ぞくパタでは次のように説明している。


「ある時点 t で状態 ω_j に到達し、かつ x_t が観測される確率を考える。その確率は、1時点前の状態が ω_1〜ω_c のいずれであるかによって異なり、c 種存在する。その中で最大となる確率を ψ_t(j) で表す」


正直、この説明では ψ_t(j) がなんなのかわからなかった。
ψ_t(j) が何かわからなくても、(8.24) 式では漸化式が与えられるので、それを認めれば計算はできる。
が、その ψ_t(j) で最適な状態の系列が推定できることはなぜわかるのだろう?


ψ_t(j) が定義されていれば、漸化式はそこから導ける(ここでも条件付き独立性を使う)し、それを求められれば最適系列を与えることもわかる。問題解決。
というわけで定義をしよう。*1

  • \psi_t(j)\;:=\;\max_{s_1,\cdots,s_{t-1}}\;P(x_1,\cdots,x_t,s_1,\cdots,s_{t-1},s_t=\omega_j)


ちなみに Ψ_t(j) の方は、この最大値を与えるときの s_{t-1}=\omega_i のインデックス i である。このことをふまえると、Ψ を逆向きにたどると最適系列を得られることも間違いなく理解できる。


条件付き独立性についてはグラフィカルモデルをやってないという同情点があったが、こちらは単なる定義の明記漏れである。
次版があるならぜひ改善してほしい……と思うが、話を聞くとそういったフィードバックは受け付けられにくそうな雰囲気なので、期待はしない。


次回読書会参加は……ツッコミの反応イマイチだったし、11章のディリクレ過程と、あともしかして行くとしたら12章かなあ。
teramonagi さんと sfchaos さんはガチツッコミしてもいいと聞いているのでw

*1:この定義を見ると、本の「説明」が実はまずいということもわかる

プチコン3号 ショートサンプル&ドリル / 十字キー編

プチコン3号の短いサンプルプログラムと、かんたんな演習問題。ドリルを解くと、だんだんプチコンプログラミングを覚える、みたいな。
「新しい命令」は プチコンで 3DS のゲームを作ろう #petitcom - 木曜不足 に出てきてないもの。

■サンプル 1. 十字キーで移動(斜め移動できない)

X=200:Y=120
SPSET 0,600  ' 600 は好きなキャラクタの番号に変えていいよ
WHILE 1
  B=BUTTON()
  IF B==1 THEN Y=Y-1
  IF B==2 THEN Y=Y+1
  IF B==4 THEN X=X-1
  IF B==8 THEN X=X+1
  SPOFS 0,X,Y
  VSYNC
WEND

キャラクタ(スプライト)の番号はスマイルツールか プチコンの標準 BG とスプライトの一覧を作ってみた #petitcom - 木曜不足 を見てね。


▼新しい命令:

BUTTON(ボタン)
どのボタンが押されているかを返す関数。引数の意味は SmileBasic – 『プチコン3号 SmileBASIC』公式サイト 参照。押されているボタンの数値の足し算を返す。A と B の両方が押されていたら(それ以外が押されてなかったら) 16+32=48 が返る。
十字ボタン上 1
十字ボタン下 2
十字ボタン左 4
十字ボタン右 8
A 16
B 32
X 64
Y 128
L 256
R 512
ZR(拡張パッド) 2048
ZL(拡張パッド) 4096


▼問題:

  • [問題 1-1] 作ったけど、移動がものすごく遅い。速さを3倍にするには?
    • (ヒント) IF B==1 THEN Y=Y-1 は「上を押されたら、上に 1 動く」。動くのを3倍にするには? 新しい命令: なし
  • [問題 1-2] ずっと左を押してると画面から出て行ってしまう。画面のはしっこにきたらそれ以上左に行かないようにするには?
    • (ヒント) 今は「左が押されたら X=X-1」だけなので画面の外に行けてしまう。「左が押されていて、まだ左に行っても大丈夫なら X=X-1」に書き換える。「まだ左に行っても大丈夫」を X の式で書くと? 新しい命令: &&
&&(アンド)
IF 文で「〜と〜の両方が成り立ったら」という、2つ(以上)の条件を書きたいときに使う。「 IF A>0 && B>0 THEN 」と書くと「 A>0 かつ B>0 なら〜」という意味になる
    • (さらにヒント) X や Y がどうなったら「画面のはしっこ」なのかわからないってときは、 X と Y の値を画面に表示してみるといい。 VSYNC の前に次の1行を入れてみよう。
  CLS:PRINT X,Y
PRINT
画面に文字や数を表示する
CLS
画面の文字を消す。 ACLS はスプライトなども全部(オール)消すけど、 CLS が消すのは文字だけ。

■サンプル 2. 十字キーで移動(斜め移動できる)


サンプル 1 では、上ボタンと右ボタンを同時に押したら動かない。
斜めに動いて欲しいとか、「右に移動しながら A ボタンで撃ちたい」とか、そういうのは AND を使って、 IF の条件 B==1 を (B AND 1)==1 のように書き直すと思った通りの動きになる。

X=200:Y=120
SPSET 0,600
WHILE 1
  B=BUTTON()
  IF (B AND 1)==1 THEN Y=Y-1
  IF (B AND 2)==2 THEN Y=Y+1
  IF (B AND 4)==4 THEN X=X-1
  IF (B AND 8)==8 THEN X=X+1
  SPOFS 0,X,Y
  VSYNC
WEND


▼新しい命令:

AND(アンド)
くわしく説明するには「二進数」が必要。学研のまんが「算数頭をつくるひみつ」は「二進数」をやさしくおもしろく説明してくれている。ここではかんたんな説明をする。
たとえば「下ボタンだけ」を押したら B=BUTTON() は 2 だけど、左を押しながらだと 2+4=6 、B と同時押しだと 2+32=34 になる。つまり IF B==2 THEN は「ほかのボタンはどれも押さず、下ボタンだけが押されたら〜」ということ。だから斜めに行けないし、「移動しながら弾をうつ」とかもできない。
IF B==6 THEN, IF B==34 THEN って押されるかもしれない全部のパターンを並べる……なんて無理だから、ここはさっきのような足し算の中にどんな数が入っているのかわかる AND を使う(なんでもわかるわけじゃなくて、この場合はわかるってこと。ここが「二進数」のポイント)。2 AND 2 も、 6 AND 2 も、34 AND 2 も全部、下ボタンが押されていれば B AND 2 は 2 になるけど、下ボタンが押されてなかったら、ほかにどんなボタンが押されていても B AND 2 は 0 になる。
だから、IF (B AND 2)==2 THEN って書けば、「ほかに押されたボタンがあってもなくても、下ボタンが押されたら〜」って命令になる。(B AND 2) にカッコが付くのは、AND より == の方が優先順位が高く、カッコを忘れたら B AND (2==2) になってしまうから。
AND も && も「アンド」なのはたまたまじゃあなくて関係があるんだけど、「二進数」について知らないあいだは別のもの( && は複数の条件を並べるとき用、AND は BUTTON 用)と思っておいた方が間違いがない。プチコン3号のよく使う命令で、「二進数」のことを知ってた方がうれしいのは BUTTON くらいしかないから、急がなくても大丈夫。*1


▼問題:

  • [問題 2-1] 移動がやっぱり遅い。速さを3倍にするには?
    • (ヒント) 問題 1-1 といっしょ。
  • [問題 2-2] いつでも速いと操作しにくい。 B ダッシュ、つまり B ボタンを押しているときだけ速さを3倍にするには?
    • (ヒント) 移動の大きさを V=1 のように変数を使って表せば、B ボタンが押されたときに V に別の値を入れるだけでスピードを変えることできる! でも IF B==32 THEN だと「 B ボタンだけが押されたら」になってしまうので……。新しい命令: なし
  • [問題 2-3] どっちに動かしてもキャラクタの向きが変わらないのはつまらない。上に動かしたときは上に、右に動かしたときは右に向かせたい。
    • (ヒント) お姫様のスプライトの番号は、右=596、下=600、左=604、上=608 である。今は SPSET 0,600 ってなっているからずっと下(正面)向きなわけ。移動の大きさを変えたいときは V=1 と新しい変数を作った。向きに合わせてスプライトの番号を変えたいときは……。新しい命令:IF THEN のあとの :(コロン)
IF(イフ) 〜 THEN(ゼン)
THEN の後ろの命令は条件が成立したら実行されるけど、そこに2つ以上の命令を書きたいとき、命令の後に :(コロン) で区切れば、2つ目・3つ目の命令を続けることができる。
  • [問題 2-4](ちょい難) 向きは変わるようになったけど、動きがないからスケートで滑ってるみたい。歩いてるみたいに足を動かしたい。
    • (ヒント) スプライトの番号を 596,597,598,599 の順番に変えて、599 の次は 596 に戻せば右に歩いているように見える。同じように下は 600,601,602,603、左は 604,605,606,607、上は 608,609,610,611。向きごとにプログラムを書いたら大変そうだけど、よく見ると全部、(最初の数), (最初の数)+1, (最初の数)+2, (最初の数)+3 ってなっている。つまり、[2-3] で用意した変数と、ループを回るたびに 0, 1, 2, 3 の順番に変わる ( 3 のあとは 0 に戻る) 変数があれば……。

*1:次の機会は SPCOL のマスクを使いたいときだろう。